pyhiscat.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhiscat.f

    
C*********************************************************************  
    
      SUBROUTINE pyhiscat 
    
C...Finds outgoing flavours and event type; sets up the kinematics  
C...and colour flow of the hard scattering. 
      common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
      SAVE /lujets/ 
      common/ludat1/mstu(200),paru(200),mstj(200),parj(200) 
      SAVE /ludat1/ 
      common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)    
      SAVE /ludat2/ 
      common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)    
      SAVE /ludat3/ 
      common/pyhisubs/msel,msub(200),kfin(2,-40:40),ckin(200) 
      SAVE /pyhisubs/ 
      common/pyhipars/mstp(200),parp(200),msti(200),pari(200) 
      SAVE /pyhipars/ 
      common/pyhiint1/mint(400),vint(400) 
      SAVE /pyhiint1/ 
      common/pyhiint2/iset(200),kfpr(200,2),coef(200,20),icol(40,4,2) 
      SAVE /pyhiint2/ 
      common/pyhiint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)  
      SAVE /pyhiint3/ 
      common/pyhiint4/widp(21:40,0:40),wide(21:40,0:40),wids(21:40,3) 
      SAVE /pyhiint4/ 
      common/pyhiint5/ngen(0:200,3),xsec(0:200,3) 
      SAVE /pyhiint5/ 
      dimension wdtp(0:40),wdte(0:40,0:5),pmq(2),z(2),cthe(2),phi(2)    
    
C...Choice of subprocess, number of documentation lines.    
      isub=mint(1)  
      idoc=6+iset(isub) 
      IF(isub.EQ.95) idoc=8 
      mint(3)=idoc-6    
      IF(idoc.GE.9) idoc=idoc+2 
      mint(4)=idoc  
      ipu1=mint(84)+1   
      ipu2=mint(84)+2   
      ipu3=mint(84)+3   
      ipu4=mint(84)+4   
      ipu5=mint(84)+5   
      ipu6=mint(84)+6   
    
C...Reset K, P and V vectors. Store incoming particles. 
      DO 100 jt=1,mstp(126)+10  
      i=mint(83)+jt 
      DO 100 j=1,5  
      k(i,j)=0  
      p(i,j)=0. 
  100 v(i,j)=0. 
      DO 110 jt=1,2 
      i=mint(83)+jt 
      k(i,1)=21 
      k(i,2)=mint(10+jt)    
      p(i,1)=0. 
      p(i,2)=0. 
      p(i,5)=vint(2+jt) 
      p(i,3)=vint(5)*(-1)**(jt+1)   
  110 p(i,4)=sqrt(p(i,3)**2+p(i,5)**2)  
      mint(6)=2 
      kfres=0   
    
C...Store incoming partons in their CM-frame.   
      sh=vint(44)   
      shr=sqrt(sh)  
      shp=vint(26)*vint(2)  
      shpr=sqrt(shp)    
      shuser=shr    
      IF(iset(isub).GE.3) shuser=shpr   
      DO 120 jt=1,2 
      i=mint(84)+jt 
      k(i,1)=14 
      k(i,2)=mint(14+jt)    
      k(i,3)=mint(83)+2+jt  
  120 p(i,5)=ulmass(k(i,2)) 
      IF(p(ipu1,5)+p(ipu2,5).GE.shuser) THEN    
        p(ipu1,5)=0.    
        p(ipu2,5)=0.    
      ENDIF 
      p(ipu1,4)=0.5*(shuser+(p(ipu1,5)**2-p(ipu2,5)**2)/shuser) 
      p(ipu1,3)=sqrt(max(0.,p(ipu1,4)**2-p(ipu1,5)**2)) 
      p(ipu2,4)=shuser-p(ipu1,4)    
      p(ipu2,3)=-p(ipu1,3)  
    
C...Copy incoming partons to documentation lines.   
      DO 130 jt=1,2 
      i1=mint(83)+4+jt  
      i2=mint(84)+jt    
      k(i1,1)=21    
      k(i1,2)=k(i2,2)   
      k(i1,3)=i1-2  
      DO 130 j=1,5  
  130 p(i1,j)=p(i2,j)   
    
C...Choose new quark flavour for relevant annihilation graphs.  
      IF(isub.EQ.12.OR.isub.EQ.53) THEN 
        CALL pyhiwidt(21,shr,wdtp,wdte)   
        rkfl=(wdte(0,1)+wdte(0,2)+wdte(0,4))*rlu(0) 
        DO 140 i=1,2*mstp(1)    
        kflq=i  
        rkfl=rkfl-(wdte(i,1)+wdte(i,2)+wdte(i,4))   
        IF(rkfl.LE.0.) goto 150 
  140   CONTINUE    
  150   CONTINUE    
      ENDIF 
    
C...Final state flavours and colour flow: default values.   
      js=1  
      mint(21)=mint(15) 
      mint(22)=mint(16) 
      mint(23)=0    
      mint(24)=0    
      kcc=20    
      kcs=isign(1,mint(15)) 
    
      IF(isub.LE.10) THEN   
      IF(isub.EQ.1) THEN    
C...f + fb -> gamma*/Z0.    
        kfres=23    
    
      ELSEIF(isub.EQ.2) THEN    
C...f + fb' -> W+/- .   
        kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))   
        kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))   
        kfres=isign(24,kch1+kch2)   
    
      ELSEIF(isub.EQ.3) THEN    
C...f + fb -> H0.   
        kfres=25    
    
      ELSEIF(isub.EQ.4) THEN    
C...gamma + W+/- -> W+/-.   
    
      ELSEIF(isub.EQ.5) THEN    
C...Z0 + Z0 -> H0.  
        xh=sh/shp   
        mint(21)=mint(15)   
        mint(22)=mint(16)   
        pmq(1)=ulmass(mint(21)) 
        pmq(2)=ulmass(mint(22)) 
  240   jt=int(1.5+rlu(0))  
        zmin=2.*pmq(jt)/shpr    
        zmax=1.-pmq(3-jt)/shpr-(sh-pmq(jt)**2)/(shpr*(shpr-pmq(3-jt)))  
        zmax=min(1.-xh,zmax)    
        z(jt)=zmin+(zmax-zmin)*rlu(0)   
        IF(-1.+(1.+xh)/(1.-z(jt))-xh/(1.-z(jt))**2.LT.  
     &  (1.-xh)**2/(4.*xh)*rlu(0)) goto 240 
        sqc1=1.-4.*pmq(jt)**2/(z(jt)**2*shp)    
        IF(sqc1.LT.1.e-8) goto 240  
        c1=sqrt(sqc1)   
        c2=1.+2.*(pmas(23,1)**2-pmq(jt)**2)/(z(jt)*shp) 
        cthe(jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1   
        cthe(jt)=min(1.,max(-1.,cthe(jt)))  
        z(3-jt)=1.-xh/(1.-z(jt))    
        sqc1=1.-4.*pmq(3-jt)**2/(z(3-jt)**2*shp)    
        IF(sqc1.LT.1.e-8) goto 240  
        c1=sqrt(sqc1)   
        c2=1.+2.*(pmas(23,1)**2-pmq(3-jt)**2)/(z(3-jt)*shp) 
        cthe(3-jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1 
        cthe(3-jt)=min(1.,max(-1.,cthe(3-jt)))  
        phir=paru(2)*rlu(0) 
        cphi=cos(phir)  
        ang=cthe(1)*cthe(2)-sqrt(1.-cthe(1)**2)*sqrt(1.-cthe(2)**2)*cphi    
        z1=2.-z(jt) 
        z2=ang*sqrt(z(jt)**2-4.*pmq(jt)**2/shp) 
        z3=1.-z(jt)-xh+(pmq(1)**2+pmq(2)**2)/shp    
        z(3-jt)=2./(z1**2-z2**2)*(z1*z3+z2*sqrt(z3**2-(z1**2-z2**2)*    
     &  pmq(3-jt)**2/shp))  
        zmin=2.*pmq(3-jt)/shpr  
        zmax=1.-pmq(jt)/shpr-(sh-pmq(3-jt)**2)/(shpr*(shpr-pmq(jt)))    
        zmax=min(1.-xh,zmax)    
        IF(z(3-jt).LT.zmin.OR.z(3-jt).GT.zmax) goto 240 
        kcc=22  
        kfres=25    
    
      ELSEIF(isub.EQ.6) THEN    
C...Z0 + W+/- -> W+/-.  
    
      ELSEIF(isub.EQ.7) THEN    
C...W+ + W- -> Z0.  
    
      ELSEIF(isub.EQ.8) THEN    
C...W+ + W- -> H0.  
        xh=sh/shp   
  250   DO 280 jt=1,2   
        i=mint(14+jt)   
        ia=iabs(i)  
        IF(ia.LE.10) THEN   
          rvckm=vint(180+i)*rlu(0)  
          DO 270 j=1,mstp(1)    
          ib=2*j-1+mod(ia,2)    
          ipm=(5-isign(1,i))/2  
          idc=j+mdcy(ia,2)+2    
          IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 270  
          mint(20+jt)=isign(ib,i)   
          rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2)   
          IF(rvckm.LE.0.) goto 280  
  270     CONTINUE  
        ELSE    
          ib=2*((ia+1)/2)-1+mod(ia,2)   
          mint(20+jt)=isign(ib,i)   
        ENDIF   
  280   pmq(jt)=ulmass(mint(20+jt)) 
        jt=int(1.5+rlu(0))  
        zmin=2.*pmq(jt)/shpr    
        zmax=1.-pmq(3-jt)/shpr-(sh-pmq(jt)**2)/(shpr*(shpr-pmq(3-jt)))  
        zmax=min(1.-xh,zmax)    
        z(jt)=zmin+(zmax-zmin)*rlu(0)   
        IF(-1.+(1.+xh)/(1.-z(jt))-xh/(1.-z(jt))**2.LT.  
     &  (1.-xh)**2/(4.*xh)*rlu(0)) goto 250 
        sqc1=1.-4.*pmq(jt)**2/(z(jt)**2*shp)    
        IF(sqc1.LT.1.e-8) goto 250  
        c1=sqrt(sqc1)   
        c2=1.+2.*(pmas(24,1)**2-pmq(jt)**2)/(z(jt)*shp) 
        cthe(jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1   
        cthe(jt)=min(1.,max(-1.,cthe(jt)))  
        z(3-jt)=1.-xh/(1.-z(jt))    
        sqc1=1.-4.*pmq(3-jt)**2/(z(3-jt)**2*shp)    
        IF(sqc1.LT.1.e-8) goto 250  
        c1=sqrt(sqc1)   
        c2=1.+2.*(pmas(24,1)**2-pmq(3-jt)**2)/(z(3-jt)*shp) 
        cthe(3-jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1 
        cthe(3-jt)=min(1.,max(-1.,cthe(3-jt)))  
        phir=paru(2)*rlu(0) 
        cphi=cos(phir)  
        ang=cthe(1)*cthe(2)-sqrt(1.-cthe(1)**2)*sqrt(1.-cthe(2)**2)*cphi    
        z1=2.-z(jt) 
        z2=ang*sqrt(z(jt)**2-4.*pmq(jt)**2/shp) 
        z3=1.-z(jt)-xh+(pmq(1)**2+pmq(2)**2)/shp    
        z(3-jt)=2./(z1**2-z2**2)*(z1*z3+z2*sqrt(z3**2-(z1**2-z2**2)*    
     &  pmq(3-jt)**2/shp))  
        zmin=2.*pmq(3-jt)/shpr  
        zmax=1.-pmq(jt)/shpr-(sh-pmq(3-jt)**2)/(shpr*(shpr-pmq(jt)))    
        zmax=min(1.-xh,zmax)    
        IF(z(3-jt).LT.zmin.OR.z(3-jt).GT.zmax) goto 250 
        kcc=22  
        kfres=25    
      ENDIF 
    
      ELSEIF(isub.LE.20) THEN   
      IF(isub.EQ.11) THEN   
C...f + f' -> f + f'; th = (p(f)-p(f))**2.  
        kcc=mint(2) 
        IF(mint(15)*mint(16).LT.0) kcc=kcc+2    
    
      ELSEIF(isub.EQ.12) THEN   
C...f + fb -> f' + fb'; th = (p(f)-p(f'))**2.   
        mint(21)=isign(kflq,mint(15))   
        mint(22)=-mint(21)  
        kcc=4   
    
      ELSEIF(isub.EQ.13) THEN   
C...f + fb -> g + g; th arbitrary.  
        mint(21)=21 
        mint(22)=21 
        kcc=mint(2)+4   
    
      ELSEIF(isub.EQ.14) THEN   
C...f + fb -> g + gam; th arbitrary.    
        IF(rlu(0).GT.0.5) js=2  
        mint(20+js)=21  
        mint(23-js)=22  
        kcc=17+js   
    
      ELSEIF(isub.EQ.15) THEN   
C...f + fb -> g + Z0; th arbitrary. 
        IF(rlu(0).GT.0.5) js=2  
        mint(20+js)=21  
        mint(23-js)=23  
        kcc=17+js   
    
      ELSEIF(isub.EQ.16) THEN   
C...f + fb' -> g + W+/-; th = (p(f)-p(W-))**2 or (p(fb')-p(W+))**2. 
        kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))   
        kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))   
        IF(mint(15)*(kch1+kch2).LT.0) js=2  
        mint(20+js)=21  
        mint(23-js)=isign(24,kch1+kch2) 
        kcc=17+js   
    
      ELSEIF(isub.EQ.17) THEN   
C...f + fb -> g + H0; th arbitrary. 
        IF(rlu(0).GT.0.5) js=2  
        mint(20+js)=21  
        mint(23-js)=25  
        kcc=17+js   
    
      ELSEIF(isub.EQ.18) THEN   
C...f + fb -> gamma + gamma; th arbitrary.  
        mint(21)=22 
        mint(22)=22 
    
      ELSEIF(isub.EQ.19) THEN   
C...f + fb -> gamma + Z0; th arbitrary. 
        IF(rlu(0).GT.0.5) js=2  
        mint(20+js)=22  
        mint(23-js)=23  
    
      ELSEIF(isub.EQ.20) THEN   
C...f + fb' -> gamma + W+/-; th = (p(f)-p(W-))**2 or (p(fb')-p(W+))**2. 
        kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))   
        kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))   
        IF(mint(15)*(kch1+kch2).LT.0) js=2  
        mint(20+js)=22  
        mint(23-js)=isign(24,kch1+kch2) 
      ENDIF 
    
      ELSEIF(isub.LE.30) THEN   
      IF(isub.EQ.21) THEN   
C...f + fb -> gamma + H0; th arbitrary. 
        IF(rlu(0).GT.0.5) js=2  
        mint(20+js)=22  
        mint(23-js)=25  
    
      ELSEIF(isub.EQ.22) THEN   
C...f + fb -> Z0 + Z0; th arbitrary.    
        mint(21)=23 
        mint(22)=23 
    
      ELSEIF(isub.EQ.23) THEN   
C...f + fb' -> Z0 + W+/-; th = (p(f)-p(W-))**2 or (p(fb')-p(W+))**2.    
        kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))   
        kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))   
        IF(mint(15)*(kch1+kch2).LT.0) js=2  
        mint(20+js)=23  
        mint(23-js)=isign(24,kch1+kch2) 
    
      ELSEIF(isub.EQ.24) THEN   
C...f + fb -> Z0 + H0; th arbitrary.    
        IF(rlu(0).GT.0.5) js=2  
        mint(20+js)=23  
        mint(23-js)=25  
    
      ELSEIF(isub.EQ.25) THEN   
C...f + fb -> W+ + W-; th = (p(f)-p(W-))**2.    
        mint(21)=-isign(24,mint(15))    
        mint(22)=-mint(21)  
    
      ELSEIF(isub.EQ.26) THEN   
C...f + fb' -> W+/- + H0; th = (p(f)-p(W-))**2 or (p(fb')-p(W+))**2.    
        kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))   
        kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))   
        IF(mint(15)*(kch1+kch2).GT.0) js=2  
        mint(20+js)=isign(24,kch1+kch2) 
        mint(23-js)=25  
    
      ELSEIF(isub.EQ.27) THEN   
C...f + fb -> H0 + H0.  
    
      ELSEIF(isub.EQ.28) THEN   
C...f + g -> f + g; th = (p(f)-p(f))**2.    
        kcc=mint(2)+6   
        IF(mint(15).EQ.21) kcc=kcc+2    
        IF(mint(15).NE.21) kcs=isign(1,mint(15))    
        IF(mint(16).NE.21) kcs=isign(1,mint(16))    
    
      ELSEIF(isub.EQ.29) THEN   
C...f + g -> f + gamma; th = (p(f)-p(f))**2.    
        IF(mint(15).EQ.21) js=2 
        mint(23-js)=22  
        kcc=15+js   
        kcs=isign(1,mint(14+js))    
    
      ELSEIF(isub.EQ.30) THEN   
C...f + g -> f + Z0; th = (p(f)-p(f))**2.   
        IF(mint(15).EQ.21) js=2 
        mint(23-js)=23  
        kcc=15+js   
        kcs=isign(1,mint(14+js))    
      ENDIF 
    
      ELSEIF(isub.LE.40) THEN   
      IF(isub.EQ.31) THEN   
C...f + g -> f' + W+/-; th = (p(f)-p(f'))**2; choose flavour f'.    
        IF(mint(15).EQ.21) js=2 
        i=mint(14+js)   
        ia=iabs(i)  
        mint(23-js)=isign(24,kchg(ia,1)*i)  
        rvckm=vint(180+i)*rlu(0)    
        DO 220 j=1,mstp(1)  
        ib=2*j-1+mod(ia,2)  
        ipm=(5-isign(1,i))/2    
        idc=j+mdcy(ia,2)+2  
        IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 220    
        mint(20+js)=isign(ib,i) 
        rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2) 
        IF(rvckm.LE.0.) goto 230    
  220   CONTINUE    
  230   kcc=15+js   
        kcs=isign(1,mint(14+js))    
    
      ELSEIF(isub.EQ.32) THEN   
C...f + g -> f + H0; th = (p(f)-p(f))**2.   
        IF(mint(15).EQ.21) js=2 
        mint(23-js)=25  
        kcc=15+js   
        kcs=isign(1,mint(14+js))    
    
      ELSEIF(isub.EQ.33) THEN   
C...f + gamma -> f + g. 
    
      ELSEIF(isub.EQ.34) THEN   
C...f + gamma -> f + gamma. 
    
      ELSEIF(isub.EQ.35) THEN   
C...f + gamma -> f + Z0.    
    
      ELSEIF(isub.EQ.36) THEN   
C...f + gamma -> f' + W+/-. 
    
      ELSEIF(isub.EQ.37) THEN   
C...f + gamma -> f + H0.    
    
      ELSEIF(isub.EQ.38) THEN   
C...f + Z0 -> f + g.    
    
      ELSEIF(isub.EQ.39) THEN   
C...f + Z0 -> f + gamma.    
    
      ELSEIF(isub.EQ.40) THEN   
C...f + Z0 -> f + Z0.   
      ENDIF 
    
      ELSEIF(isub.LE.50) THEN   
      IF(isub.EQ.41) THEN   
C...f + Z0 -> f' + W+/-.    
    
      ELSEIF(isub.EQ.42) THEN   
C...f + Z0 -> f + H0.   
    
      ELSEIF(isub.EQ.43) THEN   
C...f + W+/- -> f' + g. 
    
      ELSEIF(isub.EQ.44) THEN   
C...f + W+/- -> f' + gamma. 
    
      ELSEIF(isub.EQ.45) THEN   
C...f + W+/- -> f' + Z0.    
    
      ELSEIF(isub.EQ.46) THEN   
C...f + W+/- -> f' + W+/-.  
    
      ELSEIF(isub.EQ.47) THEN   
C...f + W+/- -> f' + H0.    
    
      ELSEIF(isub.EQ.48) THEN   
C...f + H0 -> f + g.    
    
      ELSEIF(isub.EQ.49) THEN   
C...f + H0 -> f + gamma.    
    
      ELSEIF(isub.EQ.50) THEN   
C...f + H0 -> f + Z0.   
      ENDIF 
    
      ELSEIF(isub.LE.60) THEN   
      IF(isub.EQ.51) THEN   
C...f + H0 -> f' + W+/-.    
    
      ELSEIF(isub.EQ.52) THEN   
C...f + H0 -> f + H0.   
    
      ELSEIF(isub.EQ.53) THEN   
C...g + g -> f + fb; th arbitrary.  
        kcs=(-1)**int(1.5+rlu(0))   
        mint(21)=isign(kflq,kcs)    
        mint(22)=-mint(21)  
        kcc=mint(2)+10  
    
      ELSEIF(isub.EQ.54) THEN   
C...g + gamma -> f + fb.    
    
      ELSEIF(isub.EQ.55) THEN   
C...g + Z0 -> f + fb.   
    
      ELSEIF(isub.EQ.56) THEN   
C...g + W+/- -> f + fb'.    
    
      ELSEIF(isub.EQ.57) THEN   
C...g + H0 -> f + fb.   
    
      ELSEIF(isub.EQ.58) THEN   
C...gamma + gamma -> f + fb.    
    
      ELSEIF(isub.EQ.59) THEN   
C...gamma + Z0 -> f + fb.   
    
      ELSEIF(isub.EQ.60) THEN   
C...gamma + W+/- -> f + fb'.    
      ENDIF 
    
      ELSEIF(isub.LE.70) THEN   
      IF(isub.EQ.61) THEN   
C...gamma + H0 -> f + fb.   
    
      ELSEIF(isub.EQ.62) THEN   
C...Z0 + Z0 -> f + fb.  
    
      ELSEIF(isub.EQ.63) THEN   
C...Z0 + W+/- -> f + fb'.   
    
      ELSEIF(isub.EQ.64) THEN   
C...Z0 + H0 -> f + fb.  
    
      ELSEIF(isub.EQ.65) THEN   
C...W+ + W- -> f + fb.  
    
      ELSEIF(isub.EQ.66) THEN   
C...W+/- + H0 -> f + fb'.   
    
      ELSEIF(isub.EQ.67) THEN   
C...H0 + H0 -> f + fb.  
    
      ELSEIF(isub.EQ.68) THEN   
C...g + g -> g + g; th arbitrary.   
        kcc=mint(2)+12  
        kcs=(-1)**int(1.5+rlu(0))   
    
      ELSEIF(isub.EQ.69) THEN   
C...gamma + gamma -> W+ + W-.   
    
      ELSEIF(isub.EQ.70) THEN   
C...gamma + W+/- -> gamma + W+/-    
      ENDIF 
    
      ELSEIF(isub.LE.80) THEN   
      IF(isub.EQ.71.OR.isub.EQ.72) THEN 
C...Z0 + Z0 -> Z0 + Z0; Z0 + Z0 -> W+ + W-. 
        xh=sh/shp   
        mint(21)=mint(15)   
        mint(22)=mint(16)   
        pmq(1)=ulmass(mint(21)) 
        pmq(2)=ulmass(mint(22)) 
  290   jt=int(1.5+rlu(0))  
        zmin=2.*pmq(jt)/shpr    
        zmax=1.-pmq(3-jt)/shpr-(sh-pmq(jt)**2)/(shpr*(shpr-pmq(3-jt)))  
        zmax=min(1.-xh,zmax)    
        z(jt)=zmin+(zmax-zmin)*rlu(0)   
        IF(-1.+(1.+xh)/(1.-z(jt))-xh/(1.-z(jt))**2.LT.  
     &  (1.-xh)**2/(4.*xh)*rlu(0)) goto 290 
        sqc1=1.-4.*pmq(jt)**2/(z(jt)**2*shp)    
        IF(sqc1.LT.1.e-8) goto 290  
        c1=sqrt(sqc1)   
        c2=1.+2.*(pmas(23,1)**2-pmq(jt)**2)/(z(jt)*shp) 
        cthe(jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1   
        cthe(jt)=min(1.,max(-1.,cthe(jt)))  
        z(3-jt)=1.-xh/(1.-z(jt))    
        sqc1=1.-4.*pmq(3-jt)**2/(z(3-jt)**2*shp)    
        IF(sqc1.LT.1.e-8) goto 290  
        c1=sqrt(sqc1)   
        c2=1.+2.*(pmas(23,1)**2-pmq(3-jt)**2)/(z(3-jt)*shp) 
        cthe(3-jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1 
        cthe(3-jt)=min(1.,max(-1.,cthe(3-jt)))  
        phir=paru(2)*rlu(0) 
        cphi=cos(phir)  
        ang=cthe(1)*cthe(2)-sqrt(1.-cthe(1)**2)*sqrt(1.-cthe(2)**2)*cphi    
        z1=2.-z(jt) 
        z2=ang*sqrt(z(jt)**2-4.*pmq(jt)**2/shp) 
        z3=1.-z(jt)-xh+(pmq(1)**2+pmq(2)**2)/shp    
        z(3-jt)=2./(z1**2-z2**2)*(z1*z3+z2*sqrt(z3**2-(z1**2-z2**2)*    
     &  pmq(3-jt)**2/shp))  
        zmin=2.*pmq(3-jt)/shpr  
        zmax=1.-pmq(jt)/shpr-(sh-pmq(3-jt)**2)/(shpr*(shpr-pmq(jt)))    
        zmax=min(1.-xh,zmax)    
        IF(z(3-jt).LT.zmin.OR.z(3-jt).GT.zmax) goto 290 
        kcc=22  
    
      ELSEIF(isub.EQ.73) THEN   
C...Z0 + W+/- -> Z0 + W+/-. 
        xh=sh/shp   
  300   jt=int(1.5+rlu(0))  
        i=mint(14+jt)   
        ia=iabs(i)  
        IF(ia.LE.10) THEN   
          rvckm=vint(180+i)*rlu(0)  
          DO 320 j=1,mstp(1)    
          ib=2*j-1+mod(ia,2)    
          ipm=(5-isign(1,i))/2  
          idc=j+mdcy(ia,2)+2    
          IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 320  
          mint(20+jt)=isign(ib,i)   
          rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2)   
          IF(rvckm.LE.0.) goto 330  
  320     CONTINUE  
        ELSE    
          ib=2*((ia+1)/2)-1+mod(ia,2)   
          mint(20+jt)=isign(ib,i)   
        ENDIF   
  330   pmq(jt)=ulmass(mint(20+jt)) 
        mint(23-jt)=mint(17-jt) 
        pmq(3-jt)=ulmass(mint(23-jt))   
        jt=int(1.5+rlu(0))  
        zmin=2.*pmq(jt)/shpr    
        zmax=1.-pmq(3-jt)/shpr-(sh-pmq(jt)**2)/(shpr*(shpr-pmq(3-jt)))  
        zmax=min(1.-xh,zmax)    
        z(jt)=zmin+(zmax-zmin)*rlu(0)   
        IF(-1.+(1.+xh)/(1.-z(jt))-xh/(1.-z(jt))**2.LT.  
     &  (1.-xh)**2/(4.*xh)*rlu(0)) goto 300 
        sqc1=1.-4.*pmq(jt)**2/(z(jt)**2*shp)    
        IF(sqc1.LT.1.e-8) goto 300  
        c1=sqrt(sqc1)   
        c2=1.+2.*(pmas(23,1)**2-pmq(jt)**2)/(z(jt)*shp) 
        cthe(jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1   
        cthe(jt)=min(1.,max(-1.,cthe(jt)))  
        z(3-jt)=1.-xh/(1.-z(jt))    
        sqc1=1.-4.*pmq(3-jt)**2/(z(3-jt)**2*shp)    
        IF(sqc1.LT.1.e-8) goto 300  
        c1=sqrt(sqc1)   
        c2=1.+2.*(pmas(23,1)**2-pmq(3-jt)**2)/(z(3-jt)*shp) 
        cthe(3-jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1 
        cthe(3-jt)=min(1.,max(-1.,cthe(3-jt)))  
        phir=paru(2)*rlu(0) 
        cphi=cos(phir)  
        ang=cthe(1)*cthe(2)-sqrt(1.-cthe(1)**2)*sqrt(1.-cthe(2)**2)*cphi    
        z1=2.-z(jt) 
        z2=ang*sqrt(z(jt)**2-4.*pmq(jt)**2/shp) 
        z3=1.-z(jt)-xh+(pmq(1)**2+pmq(2)**2)/shp    
        z(3-jt)=2./(z1**2-z2**2)*(z1*z3+z2*sqrt(z3**2-(z1**2-z2**2)*    
     &  pmq(3-jt)**2/shp))  
        zmin=2.*pmq(3-jt)/shpr  
        zmax=1.-pmq(jt)/shpr-(sh-pmq(3-jt)**2)/(shpr*(shpr-pmq(jt)))    
        zmax=min(1.-xh,zmax)    
        IF(z(3-jt).LT.zmin.OR.z(3-jt).GT.zmax) goto 300 
        kcc=22  
    
      ELSEIF(isub.EQ.74) THEN   
C...Z0 + H0 -> Z0 + H0. 
    
      ELSEIF(isub.EQ.75) THEN   
C...W+ + W- -> gamma + gamma.   
    
      ELSEIF(isub.EQ.76.OR.isub.EQ.77) THEN 
C...W+ + W- -> Z0 + Z0; W+ + W- -> W+ + W-. 
        xh=sh/shp   
  340   DO 370 jt=1,2   
        i=mint(14+jt)   
        ia=iabs(i)  
        IF(ia.LE.10) THEN   
          rvckm=vint(180+i)*rlu(0)  
          DO 360 j=1,mstp(1)    
          ib=2*j-1+mod(ia,2)    
          ipm=(5-isign(1,i))/2  
          idc=j+mdcy(ia,2)+2    
          IF(mdme(idc,1).NE.1.AND.mdme(idc,1).NE.ipm) goto 360  
          mint(20+jt)=isign(ib,i)   
          rvckm=rvckm-vckm((ia+1)/2,(ib+1)/2)   
          IF(rvckm.LE.0.) goto 370  
  360     CONTINUE  
        ELSE    
          ib=2*((ia+1)/2)-1+mod(ia,2)   
          mint(20+jt)=isign(ib,i)   
        ENDIF   
  370   pmq(jt)=ulmass(mint(20+jt)) 
        jt=int(1.5+rlu(0))  
        zmin=2.*pmq(jt)/shpr    
        zmax=1.-pmq(3-jt)/shpr-(sh-pmq(jt)**2)/(shpr*(shpr-pmq(3-jt)))  
        zmax=min(1.-xh,zmax)    
        z(jt)=zmin+(zmax-zmin)*rlu(0)   
        IF(-1.+(1.+xh)/(1.-z(jt))-xh/(1.-z(jt))**2.LT.  
     &  (1.-xh)**2/(4.*xh)*rlu(0)) goto 340 
        sqc1=1.-4.*pmq(jt)**2/(z(jt)**2*shp)    
        IF(sqc1.LT.1.e-8) goto 340  
        c1=sqrt(sqc1)   
        c2=1.+2.*(pmas(24,1)**2-pmq(jt)**2)/(z(jt)*shp) 
        cthe(jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1   
        cthe(jt)=min(1.,max(-1.,cthe(jt)))  
        z(3-jt)=1.-xh/(1.-z(jt))    
        sqc1=1.-4.*pmq(3-jt)**2/(z(3-jt)**2*shp)    
        IF(sqc1.LT.1.e-8) goto 340  
        c1=sqrt(sqc1)   
        c2=1.+2.*(pmas(24,1)**2-pmq(3-jt)**2)/(z(3-jt)*shp) 
        cthe(3-jt)=(c2-(c2**2-c1**2)/(c2+(2.*rlu(0)-1.)*c1))/c1 
        cthe(3-jt)=min(1.,max(-1.,cthe(3-jt)))  
        phir=paru(2)*rlu(0) 
        cphi=cos(phir)  
        ang=cthe(1)*cthe(2)-sqrt(1.-cthe(1)**2)*sqrt(1.-cthe(2)**2)*cphi    
        z1=2.-z(jt) 
        z2=ang*sqrt(z(jt)**2-4.*pmq(jt)**2/shp) 
        z3=1.-z(jt)-xh+(pmq(1)**2+pmq(2)**2)/shp    
        z(3-jt)=2./(z1**2-z2**2)*(z1*z3+z2*sqrt(z3**2-(z1**2-z2**2)*    
     &  pmq(3-jt)**2/shp))  
        zmin=2.*pmq(3-jt)/shpr  
        zmax=1.-pmq(jt)/shpr-(sh-pmq(3-jt)**2)/(shpr*(shpr-pmq(jt)))    
        zmax=min(1.-xh,zmax)    
        IF(z(3-jt).LT.zmin.OR.z(3-jt).GT.zmax) goto 340 
        kcc=22  
    
      ELSEIF(isub.EQ.78) THEN   
C...W+/- + H0 -> W+/- + H0. 
    
      ELSEIF(isub.EQ.79) THEN   
C...H0 + H0 -> H0 + H0. 
      ENDIF 
    
      ELSEIF(isub.LE.90) THEN   
      IF(isub.EQ.81) THEN   
C...q + qb -> Q' + Qb'; th = (p(q)-p(q'))**2.   
        mint(21)=isign(mint(46),mint(15))   
        mint(22)=-mint(21)  
        kcc=4   
    
      ELSEIF(isub.EQ.82) THEN   
C...g + g -> Q + Qb; th arbitrary.  
        kcs=(-1)**int(1.5+rlu(0))   
        mint(21)=isign(mint(46),kcs)    
        mint(22)=-mint(21)  
        kcc=mint(2)+10  
      ENDIF 
    
      ELSEIF(isub.LE.100) THEN  
      IF(isub.EQ.95) THEN   
C...Low-pT ( = energyless g + g -> g + g).  
        kcc=mint(2)+12  
        kcs=(-1)**int(1.5+rlu(0))   
    
      ELSEIF(isub.EQ.96) THEN   
C...Multiple interactions (should be reassigned to QCD process).    
      ENDIF 
    
      ELSEIF(isub.LE.110) THEN  
      IF(isub.EQ.101) THEN  
C...g + g -> gamma*/Z0. 
        kcc=21  
        kfres=22    
    
      ELSEIF(isub.EQ.102) THEN  
C...g + g -> H0.    
        kcc=21  
        kfres=25    
      ENDIF 
    
      ELSEIF(isub.LE.120) THEN  
      IF(isub.EQ.111) THEN  
C...f + fb -> g + H0; th arbitrary. 
        IF(rlu(0).GT.0.5) js=2  
        mint(20+js)=21  
        mint(23-js)=25  
        kcc=17+js   
    
      ELSEIF(isub.EQ.112) THEN  
C...f + g -> f + H0; th = (p(f) - p(f))**2. 
        IF(mint(15).EQ.21) js=2 
        mint(23-js)=25  
        kcc=15+js   
        kcs=isign(1,mint(14+js))    
    
      ELSEIF(isub.EQ.113) THEN  
C...g + g -> g + H0; th arbitrary.  
        IF(rlu(0).GT.0.5) js=2  
        mint(23-js)=25  
        kcc=22+js   
        kcs=(-1)**int(1.5+rlu(0))   
    
      ELSEIF(isub.EQ.114) THEN  
C...g + g -> gamma + gamma; th arbitrary.   
        IF(rlu(0).GT.0.5) js=2  
        mint(21)=22 
        mint(22)=22 
        kcc=21  
    
      ELSEIF(isub.EQ.115) THEN  
C...g + g -> gamma + Z0.    
    
      ELSEIF(isub.EQ.116) THEN  
C...g + g -> Z0 + Z0.   
    
      ELSEIF(isub.EQ.117) THEN  
C...g + g -> W+ + W-.   
      ENDIF 
    
      ELSEIF(isub.LE.140) THEN  
      IF(isub.EQ.121) THEN  
C...g + g -> f + fb + H0.   
      ENDIF 
    
      ELSEIF(isub.LE.160) THEN  
      IF(isub.EQ.141) THEN  
C...f + fb -> gamma*/Z0/Z'0.    
        kfres=32    
    
      ELSEIF(isub.EQ.142) THEN  
C...f + fb' -> H+/-.    
        kch1=kchg(iabs(mint(15)),1)*isign(1,mint(15))   
        kch2=kchg(iabs(mint(16)),1)*isign(1,mint(16))   
        kfres=isign(37,kch1+kch2)   
    
      ELSEIF(isub.EQ.143) THEN  
C...f + fb' -> R.   
        kfres=isign(40,mint(15)+mint(16))   
      ENDIF 
    
      ELSE  
      IF(isub.EQ.161) THEN  
C...g + f -> H+/- + f'; th = (p(f)-p(f))**2.    
        IF(mint(16).EQ.21) js=2 
        ia=iabs(mint(17-js))    
        mint(20+js)=isign(37,kchg(ia,1)*mint(17-js))    
        ja=ia+mod(ia,2)-mod(ia+1,2) 
        mint(23-js)=isign(ja,mint(17-js))   
        kcc=18-js   
        IF(mint(15).NE.21) kcs=isign(1,mint(15))    
        IF(mint(16).NE.21) kcs=isign(1,mint(16))    
      ENDIF 
      ENDIF 
    
      IF(idoc.EQ.7) THEN    
C...Resonance not decaying: store colour connection indices.    
        i=mint(83)+7    
        k(ipu3,1)=1 
        k(ipu3,2)=kfres 
        k(ipu3,3)=i 
        p(ipu3,4)=shuser    
        p(ipu3,5)=shuser    
        k(ipu1,4)=ipu2  
        k(ipu1,5)=ipu2  
        k(ipu2,4)=ipu1  
        k(ipu2,5)=ipu1  
        k(i,1)=21   
        k(i,2)=kfres    
        p(i,4)=shuser   
        p(i,5)=shuser   
        n=ipu3  
        mint(21)=kfres  
        mint(22)=0  
    
      ELSEIF(idoc.EQ.8) THEN    
C...2 -> 2 processes: store outgoing partons in their CM-frame. 
        DO 390 jt=1,2   
        i=mint(84)+2+jt 
        k(i,1)=1    
        IF(iabs(mint(20+jt)).LE.10.OR.mint(20+jt).EQ.21) k(i,1)=3   
        k(i,2)=mint(20+jt)  
        k(i,3)=mint(83)+idoc+jt-2   
        IF(iabs(k(i,2)).LE.10.OR.k(i,2).EQ.21) THEN 
          p(i,5)=ulmass(k(i,2)) 
        ELSE    
          p(i,5)=sqrt(vint(63+mod(js+jt,2)))    
        ENDIF   
  390   CONTINUE    
        IF(p(ipu3,5)+p(ipu4,5).GE.shr) THEN 
          kfa1=iabs(mint(21))   
          kfa2=iabs(mint(22))   
          IF((kfa1.GT.3.AND.kfa1.NE.21).OR.(kfa2.GT.3.AND.kfa2.NE.21))  
     &    THEN  
            mint(51)=1  
            RETURN  
          ENDIF 
          p(ipu3,5)=0.  
          p(ipu4,5)=0.  
        ENDIF   
        p(ipu3,4)=0.5*(shr+(p(ipu3,5)**2-p(ipu4,5)**2)/shr) 
        p(ipu3,3)=sqrt(max(0.,p(ipu3,4)**2-p(ipu3,5)**2))   
        p(ipu4,4)=shr-p(ipu3,4) 
        p(ipu4,3)=-p(ipu3,3)    
        n=ipu4  
        mint(7)=mint(83)+7  
        mint(8)=mint(83)+8  
    
C...Rotate outgoing partons using cos(theta)=(th-uh)/lam(sh,sqm3,sqm4). 
        CALL ludbrb(ipu3,ipu4,acos(vint(23)),vint(24),0d0,0d0,0d0)  
    
      ELSEIF(idoc.EQ.9) THEN    
C'''2 -> 3 processes:   
    
      ELSEIF(idoc.EQ.11) THEN   
C...Z0 + Z0 -> H0, W+ + W- -> H0: store Higgs and outgoing partons. 
        phi(1)=paru(2)*rlu(0)   
        phi(2)=phi(1)-phir  
        DO 400 jt=1,2   
        i=mint(84)+2+jt 
        k(i,1)=1    
        IF(iabs(mint(20+jt)).LE.10.OR.mint(20+jt).EQ.21) k(i,1)=3   
        k(i,2)=mint(20+jt)  
        k(i,3)=mint(83)+idoc+jt-2   
        p(i,5)=ulmass(k(i,2))   
        IF(0.5*shpr*z(jt).LE.p(i,5)) p(i,5)=0.  
        pabs=sqrt(max(0.,(0.5*shpr*z(jt))**2-p(i,5)**2))    
        ptabs=pabs*sqrt(max(0.,1.-cthe(jt)**2)) 
        p(i,1)=ptabs*cos(phi(jt))   
        p(i,2)=ptabs*sin(phi(jt))   
        p(i,3)=pabs*cthe(jt)*(-1)**(jt+1)   
        p(i,4)=0.5*shpr*z(jt)   
        izw=mint(83)+6+jt   
        k(izw,1)=21 
        k(izw,2)=23 
        IF(isub.EQ.8) k(izw,2)=isign(24,luchge(mint(14+jt)))    
        k(izw,3)=izw-2  
        p(izw,1)=-p(i,1)    
        p(izw,2)=-p(i,2)    
        p(izw,3)=(0.5*shpr-pabs*cthe(jt))*(-1)**(jt+1)  
        p(izw,4)=0.5*shpr*(1.-z(jt))    
  400   p(izw,5)=-sqrt(max(0.,p(izw,3)**2+ptabs**2-p(izw,4)**2))    
        i=mint(83)+9    
        k(ipu5,1)=1 
        k(ipu5,2)=kfres 
        k(ipu5,3)=i 
        p(ipu5,5)=shr   
        p(ipu5,1)=-p(ipu3,1)-p(ipu4,1)  
        p(ipu5,2)=-p(ipu3,2)-p(ipu4,2)  
        p(ipu5,3)=-p(ipu3,3)-p(ipu4,3)  
        p(ipu5,4)=shpr-p(ipu3,4)-p(ipu4,4)  
        k(i,1)=21   
        k(i,2)=kfres    
        DO 410 j=1,5    
  410   p(i,j)=p(ipu5,j)    
        n=ipu5  
        mint(23)=kfres  
    
      ELSEIF(idoc.EQ.12) THEN   
C...Z0 and W+/- scattering: store bosons and outgoing partons.  
        phi(1)=paru(2)*rlu(0)   
        phi(2)=phi(1)-phir  
        DO 420 jt=1,2   
        i=mint(84)+2+jt 
        k(i,1)=1    
        IF(iabs(mint(20+jt)).LE.10.OR.mint(20+jt).EQ.21) k(i,1)=3   
        k(i,2)=mint(20+jt)  
        k(i,3)=mint(83)+idoc+jt-2   
        p(i,5)=ulmass(k(i,2))   
        IF(0.5*shpr*z(jt).LE.p(i,5)) p(i,5)=0.  
        pabs=sqrt(max(0.,(0.5*shpr*z(jt))**2-p(i,5)**2))    
        ptabs=pabs*sqrt(max(0.,1.-cthe(jt)**2)) 
        p(i,1)=ptabs*cos(phi(jt))   
        p(i,2)=ptabs*sin(phi(jt))   
        p(i,3)=pabs*cthe(jt)*(-1)**(jt+1)   
        p(i,4)=0.5*shpr*z(jt)   
        izw=mint(83)+6+jt   
        k(izw,1)=21 
        IF(mint(14+jt).EQ.mint(20+jt)) THEN 
          k(izw,2)=23   
        ELSE    
          k(izw,2)=isign(24,luchge(mint(14+jt))-luchge(mint(20+jt)))    
        ENDIF   
        k(izw,3)=izw-2  
        p(izw,1)=-p(i,1)    
        p(izw,2)=-p(i,2)    
        p(izw,3)=(0.5*shpr-pabs*cthe(jt))*(-1)**(jt+1)  
        p(izw,4)=0.5*shpr*(1.-z(jt))    
        p(izw,5)=-sqrt(max(0.,p(izw,3)**2+ptabs**2-p(izw,4)**2))    
        ipu=mint(84)+4+jt   
        k(ipu,1)=3  
        k(ipu,2)=kfpr(isub,jt)  
        k(ipu,3)=mint(83)+8+jt  
        IF(iabs(k(ipu,2)).LE.10.OR.k(ipu,2).EQ.21) THEN 
          p(ipu,5)=ulmass(k(ipu,2)) 
        ELSE    
          p(ipu,5)=sqrt(vint(63+mod(js+jt,2)))  
        ENDIF   
        mint(22+jt)=k(izw,2)    
  420   CONTINUE    
        IF(isub.EQ.72) k(mint(84)+4+int(1.5+rlu(0)),2)=-24  
C...Find rotation and boost for hard scattering subsystem.  
        i1=mint(83)+7   
        i2=mint(83)+8   
        bexcm=(p(i1,1)+p(i2,1))/(p(i1,4)+p(i2,4))   
        beycm=(p(i1,2)+p(i2,2))/(p(i1,4)+p(i2,4))   
        bezcm=(p(i1,3)+p(i2,3))/(p(i1,4)+p(i2,4))   
        gamcm=(p(i1,4)+p(i2,4))/shr 
        bepcm=bexcm*p(i1,1)+beycm*p(i1,2)+bezcm*p(i1,3) 
        px=p(i1,1)+gamcm*(gamcm/(1.+gamcm)*bepcm-p(i1,4))*bexcm 
        py=p(i1,2)+gamcm*(gamcm/(1.+gamcm)*bepcm-p(i1,4))*beycm 
        pz=p(i1,3)+gamcm*(gamcm/(1.+gamcm)*bepcm-p(i1,4))*bezcm 
        thecm=ulangl(pz,sqrt(px**2+py**2))  
        phicm=ulangl(px,py) 
C...Store hard scattering subsystem. Rotate and boost it.   
        sqlam=(sh-p(ipu5,5)**2-p(ipu6,5)**2)**2-4.*p(ipu5,5)**2*    
     &  p(ipu6,5)**2    
        pabs=sqrt(max(0.,sqlam/(4.*sh)))    
        cthwz=vint(23)  
        sthwz=sqrt(max(0.,1.-cthwz**2)) 
        phiwz=vint(24)-phicm    
        p(ipu5,1)=pabs*sthwz*cos(phiwz) 
        p(ipu5,2)=pabs*sthwz*sin(phiwz) 
        p(ipu5,3)=pabs*cthwz    
        p(ipu5,4)=sqrt(pabs**2+p(ipu5,5)**2)    
        p(ipu6,1)=-p(ipu5,1)    
        p(ipu6,2)=-p(ipu5,2)    
        p(ipu6,3)=-p(ipu5,3)    
        p(ipu6,4)=sqrt(pabs**2+p(ipu6,5)**2)    
        CALL ludbrb(ipu5,ipu6,thecm,phicm,dble(bexcm),dble(beycm),  
     &  dble(bezcm))    
        DO 430 jt=1,2   
        i1=mint(83)+8+jt    
        i2=mint(84)+4+jt    
        k(i1,1)=21  
        k(i1,2)=k(i2,2) 
        DO 430 j=1,5    
  430   p(i1,j)=p(i2,j) 
        n=ipu6  
        mint(7)=mint(83)+9  
        mint(8)=mint(83)+10 
      ENDIF 
    
      IF(idoc.GE.8) THEN    
C...Store colour connection indices.    
        DO 440 j=1,2    
        jc=j    
        IF(kcs.EQ.-1) jc=3-j    
        IF(icol(kcc,1,jc).NE.0.AND.k(ipu1,1).EQ.14) k(ipu1,j+3)=    
     &  k(ipu1,j+3)+mint(84)+icol(kcc,1,jc) 
        IF(icol(kcc,2,jc).NE.0.AND.k(ipu2,1).EQ.14) k(ipu2,j+3)=    
     &  k(ipu2,j+3)+mint(84)+icol(kcc,2,jc) 
        IF(icol(kcc,3,jc).NE.0.AND.k(ipu3,1).EQ.3) k(ipu3,j+3)= 
     &  mstu(5)*(mint(84)+icol(kcc,3,jc))   
  440   IF(icol(kcc,4,jc).NE.0.AND.k(ipu4,1).EQ.3) k(ipu4,j+3)= 
     &  mstu(5)*(mint(84)+icol(kcc,4,jc))   
    
C...Copy outgoing partons to documentation lines.   
        DO 450 i=1,2    
        i1=mint(83)+idoc-2+i    
        i2=mint(84)+2+i 
        k(i1,1)=21  
        k(i1,2)=k(i2,2) 
        IF(idoc.LE.9) k(i1,3)=0 
        IF(idoc.GE.11) k(i1,3)=mint(83)+2+i 
        DO 450 j=1,5    
  450   p(i1,j)=p(i2,j) 
      ENDIF 
      mint(52)=n    
    
C...Low-pT events: remove gluons used for string drawing purposes.  
      IF(isub.EQ.95) THEN   
        k(ipu3,1)=k(ipu3,1)+10  
        k(ipu4,1)=k(ipu4,1)+10  
        DO 460 j=41,66  
  460   vint(j)=0.  
        DO 470 i=mint(83)+5,mint(83)+8  
        DO 470 j=1,5    
  470   p(i,j)=0.   
      ENDIF 
    
      RETURN    
      END