lushow.f

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

    
C*********************************************************************  
    
      SUBROUTINE lushow(IP1,IP2,QMAX)   
    
C...Purpose: to generate timelike parton showers from given partons.    
      IMPLICIT DOUBLE PRECISION(d)  
      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/ 
      dimension pmth(5,40),ps(5),pma(4),pmsd(4),iep(4),ipa(4),  
     &kfla(4),kfld(4),kfl(4),itry(4),isi(4),isl(4),dp(4),dpt(5,4)   
    
C...Initialization of cutoff masses etc.    
      IF(mstj(41).LE.0.OR.(mstj(41).EQ.1.AND.qmax.LE.parj(82)).OR.  
     &qmax.LE.min(parj(82),parj(83)).OR.mstj(41).GE.3) RETURN   
      pmth(1,21)=ulmass(21) 
      pmth(2,21)=sqrt(pmth(1,21)**2+0.25*parj(82)**2)   
      pmth(3,21)=2.*pmth(2,21)  
      pmth(4,21)=pmth(3,21) 
      pmth(5,21)=pmth(3,21) 
      pmth(1,22)=ulmass(22) 
      pmth(2,22)=sqrt(pmth(1,22)**2+0.25*parj(83)**2)   
      pmth(3,22)=2.*pmth(2,22)  
      pmth(4,22)=pmth(3,22) 
      pmth(5,22)=pmth(3,22) 
      pmqth1=parj(82)   
      IF(mstj(41).EQ.2) pmqth1=min(parj(82),parj(83))   
      pmqth2=pmth(2,21) 
      IF(mstj(41).EQ.2) pmqth2=min(pmth(2,21),pmth(2,22))   
      DO 100 if=1,8 
      pmth(1,if)=ulmass(if) 
      pmth(2,if)=sqrt(pmth(1,if)**2+0.25*pmqth1**2) 
      pmth(3,if)=pmth(2,if)+pmqth2  
      pmth(4,if)=sqrt(pmth(1,if)**2+0.25*parj(82)**2)+pmth(2,21)    
  100 pmth(5,if)=sqrt(pmth(1,if)**2+0.25*parj(83)**2)+pmth(2,22)    
      pt2min=max(0.5*parj(82),1.1*parj(81))**2  
      alams=parj(81)**2 
      alfm=log(pt2min/alams)    
    
C...Store positions of shower initiating partons.   
      m3jc=0    
      IF(ip1.GT.0.AND.ip1.LE.min(n,mstu(4)-mstu(32)).AND.ip2.EQ.0) THEN 
        npa=1   
        ipa(1)=ip1  
      ELSEIF(min(ip1,ip2).GT.0.AND.max(ip1,ip2).LE.min(n,mstu(4)-   
     &mstu(32))) THEN   
        npa=2   
        ipa(1)=ip1  
        ipa(2)=ip2  
      ELSEIF(ip1.GT.0.AND.ip1.LE.min(n,mstu(4)-mstu(32)).AND.ip2.LT.0.  
     &and.ip2.GE.-3) THEN   
        npa=iabs(ip2)   
        DO 110 i=1,npa  
  110   ipa(i)=ip1+i-1  
      ELSE  
        CALL luerrm(12, 
     &  '(LUSHOW:) failed to reconstruct showering system') 
        IF(mstu(21).GE.1) RETURN    
      ENDIF 
    
C...Check on phase space available for emission.    
      irej=0    
      DO 120 j=1,5  
  120 ps(j)=0.  
      pm=0. 
      DO 130 i=1,npa    
      kfla(i)=iabs(k(ipa(i),2)) 
      pma(i)=p(ipa(i),5)    
      IF(kfla(i).NE.0.AND.(kfla(i).LE.8.OR.kfla(i).EQ.21))  
     &pma(i)=pmth(3,kfla(i))    
      pm=pm+pma(i)  
      IF(kfla(i).EQ.0.OR.(kfla(i).GT.8.AND.kfla(i).NE.21).OR.   
     &pma(i).GT.qmax) irej=irej+1   
      DO 130 j=1,4  
  130 ps(j)=ps(j)+p(ipa(i),j)   
      IF(irej.EQ.npa) RETURN    
      ps(5)=sqrt(max(0.,ps(4)**2-ps(1)**2-ps(2)**2-ps(3)**2))   
      IF(npa.EQ.1) ps(5)=ps(4)  
      IF(ps(5).LE.pm+pmqth1) RETURN 
      IF(npa.EQ.2.AND.mstj(47).GE.1) THEN   
        IF(kfla(1).GE.1.AND.kfla(1).LE.8.AND.kfla(2).GE.1.AND.  
     &  kfla(2).LE.8) m3jc=1    
        IF(mstj(47).GE.2) m3jc=1    
      ENDIF 
    
C...Define imagined single initiator of shower for parton system.   
      ns=n  
      IF(n.GT.mstu(4)-mstu(32)-5) THEN  
        CALL luerrm(11,'(LUSHOW:) no more memory left in LUJETS')   
        IF(mstu(21).GE.1) RETURN    
      ENDIF 
      IF(npa.GE.2) THEN 
        k(n+1,1)=11 
        k(n+1,2)=21 
        k(n+1,3)=0  
        k(n+1,4)=0  
        k(n+1,5)=0  
        p(n+1,1)=0. 
        p(n+1,2)=0. 
        p(n+1,3)=0. 
        p(n+1,4)=ps(5)  
        p(n+1,5)=ps(5)  
        v(n+1,5)=ps(5)**2   
        n=n+1   
      ENDIF 
    
C...Loop over partons that may branch.  
      nep=npa   
      im=ns 
      IF(npa.EQ.1) im=ns-1  
  140 im=im+1   
      IF(n.GT.ns) THEN  
        IF(im.GT.n) goto 380    
        kflm=iabs(k(im,2))  
        IF(kflm.EQ.0.OR.(kflm.GT.8.AND.kflm.NE.21)) goto 140    
        IF(p(im,5).LT.pmth(2,kflm)) goto 140    
        igm=k(im,3) 
      ELSE  
        igm=-1  
      ENDIF 
      IF(n+nep.GT.mstu(4)-mstu(32)-5) THEN  
        CALL luerrm(11,'(LUSHOW:) no more memory left in LUJETS')   
        IF(mstu(21).GE.1) RETURN    
      ENDIF 
    
C...Position of aunt (sister to branching parton).  
C...Origin and flavour of daughters.    
      iau=0 
      IF(igm.GT.0) THEN 
        IF(k(im-1,3).EQ.igm) iau=im-1   
        IF(n.GE.im+1.AND.k(im+1,3).EQ.igm) iau=im+1 
      ENDIF 
      IF(igm.GE.0) THEN 
        k(im,4)=n+1 
        DO 150 i=1,nep  
  150   k(n+i,3)=im 
      ELSE  
        k(n+1,3)=ipa(1) 
      ENDIF 
      IF(igm.LE.0) THEN 
        DO 160 i=1,nep  
  160   k(n+i,2)=k(ipa(i),2)    
      ELSEIF(kflm.NE.21) THEN   
        k(n+1,2)=k(im,2)    
        k(n+2,2)=k(im,5)    
      ELSEIF(k(im,5).EQ.21) THEN    
        k(n+1,2)=21 
        k(n+2,2)=21 
      ELSE  
        k(n+1,2)=k(im,5)    
        k(n+2,2)=-k(im,5)   
      ENDIF 
    
C...Reset flags on daughers and tries made. 
      DO 170 ip=1,nep   
      k(n+ip,1)=3   
      k(n+ip,4)=0   
      k(n+ip,5)=0   
      kfld(ip)=iabs(k(n+ip,2))  
      itry(ip)=0    
      isl(ip)=0 
      isi(ip)=0 
  170 IF(kfld(ip).GT.0.AND.(kfld(ip).LE.8.OR.kfld(ip).EQ.21)) isi(ip)=1 
      islm=0    
    
C...Maximum virtuality of daughters.    
      IF(igm.LE.0) THEN 
        DO 180 i=1,npa  
        IF(npa.GE.3) p(n+i,4)=(ps(4)*p(ipa(i),4)-ps(1)*p(ipa(i),1)- 
     &  ps(2)*p(ipa(i),2)-ps(3)*p(ipa(i),3))/ps(5)  
        p(n+i,5)=min(qmax,ps(5))    
        IF(npa.GE.3) p(n+i,5)=min(p(n+i,5),p(n+i,4))    
  180   IF(isi(i).EQ.0) p(n+i,5)=p(ipa(i),5)    
      ELSE  
        IF(mstj(43).LE.2) pem=v(im,2)   
        IF(mstj(43).GE.3) pem=p(im,4)   
        p(n+1,5)=min(p(im,5),v(im,1)*pem)   
        p(n+2,5)=min(p(im,5),(1.-v(im,1))*pem)  
        IF(k(n+2,2).EQ.22) p(n+2,5)=pmth(1,22)  
      ENDIF 
      DO 190 i=1,nep    
      pmsd(i)=p(n+i,5)  
      IF(isi(i).EQ.1) THEN  
        IF(p(n+i,5).LE.pmth(3,kfld(i))) p(n+i,5)=pmth(1,kfld(i))    
      ENDIF 
  190 v(n+i,5)=p(n+i,5)**2  
    
C...Choose one of the daughters for evolution.  
  200 inum=0    
      IF(nep.EQ.1) inum=1   
      DO 210 i=1,nep    
  210 IF(inum.EQ.0.AND.isl(i).EQ.1) inum=i  
      DO 220 i=1,nep    
      IF(inum.EQ.0.AND.itry(i).EQ.0.AND.isi(i).EQ.1) THEN   
        IF(p(n+i,5).GE.pmth(2,kfld(i))) inum=i  
      ENDIF 
  220 CONTINUE  
      IF(inum.EQ.0) THEN    
        rmax=0. 
        DO 230 i=1,nep  
        IF(isi(i).EQ.1.AND.pmsd(i).GE.pmqth2) THEN  
          rpm=p(n+i,5)/pmsd(i)  
          IF(rpm.GT.rmax.AND.p(n+i,5).GE.pmth(2,kfld(i))) THEN  
            rmax=rpm    
            inum=i  
          ENDIF 
        ENDIF   
  230   CONTINUE    
      ENDIF 
    
C...Store information on choice of evolving daughter.   
      inum=max(1,inum)  
      iep(1)=n+inum 
      DO 240 i=2,nep    
      iep(i)=iep(i-1)+1 
  240 IF(iep(i).GT.n+nep) iep(i)=n+1    
      DO 250 i=1,nep    
  250 kfl(i)=iabs(k(iep(i),2))  
      itry(inum)=itry(inum)+1   
      IF(itry(inum).GT.200) THEN    
        CALL luerrm(14,'(LUSHOW:) caught in infinite loop') 
        IF(mstu(21).GE.1) RETURN    
      ENDIF 
      z=0.5 
      IF(kfl(1).EQ.0.OR.(kfl(1).GT.8.AND.kfl(1).NE.21)) goto 300    
      IF(p(iep(1),5).LT.pmth(2,kfl(1))) goto 300    
    
C...Calculate allowed z range.  
      IF(nep.EQ.1) THEN 
        pmed=ps(4)  
      ELSEIF(igm.EQ.0.OR.mstj(43).LE.2) THEN    
        pmed=p(im,5)    
      ELSE  
        IF(inum.EQ.1) pmed=v(im,1)*pem  
        IF(inum.EQ.2) pmed=(1.-v(im,1))*pem 
      ENDIF 
      IF(mod(mstj(43),2).EQ.1) THEN 
        zc=pmth(2,21)/pmed  
        zce=pmth(2,22)/pmed 
      ELSE  
        zc=0.5*(1.-sqrt(max(0.,1.-(2.*pmth(2,21)/pmed)**2)))    
        IF(zc.LT.1e-4) zc=(pmth(2,21)/pmed)**2  
        zce=0.5*(1.-sqrt(max(0.,1.-(2.*pmth(2,22)/pmed)**2)))   
        IF(zce.LT.1e-4) zce=(pmth(2,22)/pmed)**2    
      ENDIF 
      zc=min(zc,0.491)  
      zce=min(zce,0.491)    
      IF((mstj(41).EQ.1.AND.zc.GT.0.49).OR.(mstj(41).EQ.2.AND.  
     &min(zc,zce).GT.0.49)) THEN    
        p(iep(1),5)=pmth(1,kfl(1))  
        v(iep(1),5)=p(iep(1),5)**2  
        goto 300    
      ENDIF 
    
C...Integral of Altarelli-Parisi z kernel for QCD.  
      IF(mstj(49).EQ.0.AND.kfl(1).EQ.21) THEN   
        fbr=6.*log((1.-zc)/zc)+mstj(45)*(0.5-zc)    
      ELSEIF(mstj(49).EQ.0) THEN    
        fbr=(8./3.)*log((1.-zc)/zc) 
    
C...Integral of Altarelli-Parisi z kernel for scalar gluon. 
      ELSEIF(mstj(49).EQ.1.AND.kfl(1).EQ.21) THEN   
        fbr=(parj(87)+mstj(45)*parj(88))*(1.-2.*zc) 
      ELSEIF(mstj(49).EQ.1) THEN    
        fbr=(1.-2.*zc)/3.   
        IF(igm.EQ.0.AND.m3jc.EQ.1) fbr=4.*fbr   
    
C...Integral of Altarelli-Parisi z kernel for Abelian vector gluon. 
      ELSEIF(kfl(1).EQ.21) THEN 
        fbr=6.*mstj(45)*(0.5-zc)    
      ELSE  
        fbr=2.*log((1.-zc)/zc)  
      ENDIF 
    
C...Integral of Altarelli-Parisi kernel for photon emission.    
      IF(mstj(41).EQ.2.AND.kfl(1).GE.1.AND.kfl(1).LE.8) 
     &fbre=(kchg(kfl(1),1)/3.)**2*2.*log((1.-zce)/zce)  
    
C...Inner veto algorithm starts. Find maximum mass for evolution.   
  260 pms=v(iep(1),5)   
      IF(igm.GE.0) THEN 
        pm2=0.  
        DO 270 i=2,nep  
        pm=p(iep(i),5)  
        IF(kfl(i).GT.0.AND.(kfl(i).LE.8.OR.kfl(i).EQ.21)) pm=   
     &  pmth(2,kfl(i))  
  270   pm2=pm2+pm  
        pms=min(pms,(p(im,5)-pm2)**2)   
      ENDIF 
    
C...Select mass for daughter in QCD evolution.  
      b0=27./6. 
      DO 280 if=4,mstj(45)  
  280 IF(pms.GT.4.*pmth(2,if)**2) b0=(33.-2.*if)/6. 
      IF(mstj(44).LE.0) THEN    
        pmsqcd=pms*exp(max(-100.,log(rlu(0))*paru(2)/(paru(111)*fbr)))  
      ELSEIF(mstj(44).EQ.1) THEN    
        pmsqcd=4.*alams*(0.25*pms/alams)**(rlu(0)**(b0/fbr))    
      ELSE  
        pmsqcd=pms*rlu(0)**(alfm*b0/fbr)    
      ENDIF 
      IF(zc.GT.0.49.OR.pmsqcd.LE.pmth(4,kfl(1))**2) pmsqcd= 
     &pmth(2,kfl(1))**2 
      v(iep(1),5)=pmsqcd    
      mce=1 
    
C...Select mass for daughter in QED evolution.  
      IF(mstj(41).EQ.2.AND.kfl(1).GE.1.AND.kfl(1).LE.8) THEN    
        pmsqed=pms*exp(max(-100.,log(rlu(0))*paru(2)/(paru(101)*fbre))) 
        IF(zce.GT.0.49.OR.pmsqed.LE.pmth(5,kfl(1))**2) pmsqed=  
     &  pmth(2,kfl(1))**2   
        IF(pmsqed.GT.pmsqcd) THEN   
          v(iep(1),5)=pmsqed    
          mce=2 
        ENDIF   
      ENDIF 
    
C...Check whether daughter mass below cutoff.   
      p(iep(1),5)=sqrt(v(iep(1),5)) 
      IF(p(iep(1),5).LE.pmth(3,kfl(1))) THEN    
        p(iep(1),5)=pmth(1,kfl(1))  
        v(iep(1),5)=p(iep(1),5)**2  
        goto 300    
      ENDIF 
    
C...Select z value of branching: q -> qgamma.   
      IF(mce.EQ.2) THEN 
        z=1.-(1.-zce)*(zce/(1.-zce))**rlu(0)    
        IF(1.+z**2.LT.2.*rlu(0)) goto 260   
        k(iep(1),5)=22  
    
C...Select z value of branching: q -> qg, g -> gg, g -> qqbar.  
      ELSEIF(mstj(49).NE.1.AND.kfl(1).NE.21) THEN   
        z=1.-(1.-zc)*(zc/(1.-zc))**rlu(0)   
        IF(1.+z**2.LT.2.*rlu(0)) goto 260   
        k(iep(1),5)=21  
      ELSEIF(mstj(49).EQ.0.AND.mstj(45)*(0.5-zc).LT.rlu(0)*fbr) THEN    
        z=(1.-zc)*(zc/(1.-zc))**rlu(0)  
        IF(rlu(0).GT.0.5) z=1.-z    
        IF((1.-z*(1.-z))**2.LT.rlu(0)) goto 260 
        k(iep(1),5)=21  
      ELSEIF(mstj(49).NE.1) THEN    
        z=zc+(1.-2.*zc)*rlu(0)  
        IF(z**2+(1.-z)**2.LT.rlu(0)) goto 260   
        kflb=1+int(mstj(45)*rlu(0)) 
        pmq=4.*pmth(2,kflb)**2/v(iep(1),5)  
        IF(pmq.GE.1.) goto 260  
        pmq0=4.*pmth(2,21)**2/v(iep(1),5)   
        IF(mod(mstj(43),2).EQ.0.AND.(1.+0.5*pmq)*sqrt(1.-pmq).LT.   
     &  rlu(0)*(1.+0.5*pmq0)*sqrt(1.-pmq0)) goto 260    
        k(iep(1),5)=kflb    
    
C...Ditto for scalar gluon model.   
      ELSEIF(kfl(1).NE.21) THEN 
        z=1.-sqrt(zc**2+rlu(0)*(1.-2.*zc))  
        k(iep(1),5)=21  
      ELSEIF(rlu(0)*(parj(87)+mstj(45)*parj(88)).LE.parj(87)) THEN  
        z=zc+(1.-2.*zc)*rlu(0)  
        k(iep(1),5)=21  
      ELSE  
        z=zc+(1.-2.*zc)*rlu(0)  
        kflb=1+int(mstj(45)*rlu(0)) 
        pmq=4.*pmth(2,kflb)**2/v(iep(1),5)  
        IF(pmq.GE.1.) goto 260  
        k(iep(1),5)=kflb    
      ENDIF 
      IF(mce.EQ.1.AND.mstj(44).GE.2) THEN   
        IF(z*(1.-z)*v(iep(1),5).LT.pt2min) goto 260 
        IF(alfm/log(v(iep(1),5)*z*(1.-z)/alams).LT.rlu(0)) goto 260 
      ENDIF 
    
C...Check if z consistent with chosen m.    
      IF(kfl(1).EQ.21) THEN 
        kflgd1=iabs(k(iep(1),5))    
        kflgd2=kflgd1   
      ELSE  
        kflgd1=kfl(1)   
        kflgd2=iabs(k(iep(1),5))    
      ENDIF 
      IF(nep.EQ.1) THEN 
        ped=ps(4)   
      ELSEIF(nep.GE.3) THEN 
        ped=p(iep(1),4) 
      ELSEIF(igm.EQ.0.OR.mstj(43).LE.2) THEN    
        ped=0.5*(v(im,5)+v(iep(1),5)-pm2**2)/p(im,5)    
      ELSE  
        IF(iep(1).EQ.n+1) ped=v(im,1)*pem   
        IF(iep(1).EQ.n+2) ped=(1.-v(im,1))*pem  
      ENDIF 
      IF(mod(mstj(43),2).EQ.1) THEN 
        pmqth3=0.5*parj(82) 
        IF(kflgd2.EQ.22) pmqth3=0.5*parj(83)    
        pmq1=(pmth(1,kflgd1)**2+pmqth3**2)/v(iep(1),5)  
        pmq2=(pmth(1,kflgd2)**2+pmqth3**2)/v(iep(1),5)  
        zd=sqrt(max(0.,(1.-v(iep(1),5)/ped**2)*((1.-pmq1-pmq2)**2-  
     &  4.*pmq1*pmq2))) 
        zh=1.+pmq1-pmq2 
      ELSE  
        zd=sqrt(max(0.,1.-v(iep(1),5)/ped**2))  
        zh=1.   
      ENDIF 
      zl=0.5*(zh-zd)    
      zu=0.5*(zh+zd)    
      IF(z.LT.zl.OR.z.GT.zu) goto 260   
      IF(kfl(1).EQ.21) v(iep(1),3)=log(zu*(1.-zl)/max(1e-20,zl* 
     &(1.-zu))) 
      IF(kfl(1).NE.21) v(iep(1),3)=log((1.-zl)/max(1e-10,1.-zu))    
    
C...Three-jet matrix element correction.    
      IF(igm.EQ.0.AND.m3jc.EQ.1) THEN   
        x1=z*(1.+v(iep(1),5)/v(ns+1,5)) 
        x2=1.-v(iep(1),5)/v(ns+1,5) 
        x3=(1.-x1)+(1.-x2)  
        IF(mce.EQ.2) THEN   
          ki1=k(ipa(inum),2)    
          ki2=k(ipa(3-inum),2)  
          qf1=kchg(iabs(ki1),1)*isign(1,ki1)/3. 
          qf2=kchg(iabs(ki2),1)*isign(1,ki2)/3. 
          wshow=qf1**2*(1.-x1)/x3*(1.+(x1/(2.-x2))**2)+ 
     &    qf2**2*(1.-x2)/x3*(1.+(x2/(2.-x1))**2)    
          wme=(qf1*(1.-x1)/x3-qf2*(1.-x2)/x3)**2*(x1**2+x2**2)  
        ELSEIF(mstj(49).NE.1) THEN  
          wshow=1.+(1.-x1)/x3*(x1/(2.-x2))**2+  
     &    (1.-x2)/x3*(x2/(2.-x1))**2    
          wme=x1**2+x2**2   
        ELSE    
          wshow=4.*x3*((1.-x1)/(2.-x2)**2+(1.-x2)/(2.-x1)**2)   
          wme=x3**2 
        ENDIF   
        IF(wme.LT.rlu(0)*wshow) goto 260    
    
C...Impose angular ordering by rejection of nonordered emission.    
      ELSEIF(mce.EQ.1.AND.igm.GT.0.AND.mstj(42).GE.2) THEN  
        maom=1  
        zm=v(im,1)  
        IF(iep(1).EQ.n+2) zm=1.-v(im,1) 
        the2id=z*(1.-z)*(zm*p(im,4))**2/v(iep(1),5) 
        iaom=im 
  290   IF(k(iaom,5).EQ.22) THEN    
          iaom=k(iaom,3)    
          IF(k(iaom,3).LE.ns) maom=0    
          IF(maom.EQ.1) goto 290    
        ENDIF   
        IF(maom.EQ.1) THEN  
          the2im=v(iaom,1)*(1.-v(iaom,1))*p(iaom,4)**2/v(iaom,5)    
          IF(the2id.LT.the2im) goto 260 
        ENDIF   
      ENDIF 
    
C...Impose user-defined maximum angle at first branching.   
      IF(mstj(48).EQ.1) THEN    
        IF(nep.EQ.1.AND.im.EQ.ns) THEN  
          the2id=z*(1.-z)*ps(4)**2/v(iep(1),5)  
          IF(the2id.LT.1./parj(85)**2) goto 260 
        ELSEIF(nep.EQ.2.AND.iep(1).EQ.ns+2) THEN    
          the2id=z*(1.-z)*(0.5*p(im,4))**2/v(iep(1),5)  
          IF(the2id.LT.1./parj(85)**2) goto 260 
        ELSEIF(nep.EQ.2.AND.iep(1).EQ.ns+3) THEN    
          the2id=z*(1.-z)*(0.5*p(im,4))**2/v(iep(1),5)  
          IF(the2id.LT.1./parj(86)**2) goto 260 
        ENDIF   
      ENDIF 
    
C...End of inner veto algorithm. Check if only one leg evolved so far.  
  300 v(iep(1),1)=z 
      isl(1)=0  
      isl(2)=0  
      IF(nep.EQ.1) goto 330 
      IF(nep.EQ.2.AND.p(iep(1),5)+p(iep(2),5).GE.p(im,5)) goto 200  
      DO 310 i=1,nep    
      IF(itry(i).EQ.0.AND.kfld(i).GT.0.AND.(kfld(i).LE.8.OR.kfld(i).EQ. 
     &21)) THEN 
        IF(p(n+i,5).GE.pmth(2,kfld(i))) goto 200    
      ENDIF 
  310 CONTINUE  
    
C...Check if chosen multiplet m1,m2,z1,z2 is physical.  
      IF(nep.EQ.3) THEN 
        pa1s=(p(n+1,4)+p(n+1,5))*(p(n+1,4)-p(n+1,5))    
        pa2s=(p(n+2,4)+p(n+2,5))*(p(n+2,4)-p(n+2,5))    
        pa3s=(p(n+3,4)+p(n+3,5))*(p(n+3,4)-p(n+3,5))    
        pts=0.25*(2.*pa1s*pa2s+2.*pa1s*pa3s+2.*pa2s*pa3s-   
     &  pa1s**2-pa2s**2-pa3s**2)/pa1s   
        IF(pts.LE.0.) goto 200  
      ELSEIF(igm.EQ.0.OR.mstj(43).LE.2.OR.mod(mstj(43),2).EQ.0) THEN    
        DO 320 i1=n+1,n+2   
        kflda=iabs(k(i1,2)) 
        IF(kflda.EQ.0.OR.(kflda.GT.8.AND.kflda.NE.21)) goto 320 
        IF(p(i1,5).LT.pmth(2,kflda)) goto 320   
        IF(kflda.EQ.21) THEN    
          kflgd1=iabs(k(i1,5))  
          kflgd2=kflgd1 
        ELSE    
          kflgd1=kflda  
          kflgd2=iabs(k(i1,5))  
        ENDIF   
        i2=2*n+3-i1 
        IF(igm.EQ.0.OR.mstj(43).LE.2) THEN  
          ped=0.5*(v(im,5)+v(i1,5)-v(i2,5))/p(im,5) 
        ELSE    
          IF(i1.EQ.n+1) zm=v(im,1)  
          IF(i1.EQ.n+2) zm=1.-v(im,1)   
          pml=sqrt((v(im,5)-v(n+1,5)-v(n+2,5))**2-  
     &    4.*v(n+1,5)*v(n+2,5)) 
          ped=pem*(0.5*(v(im,5)-pml+v(i1,5)-v(i2,5))+pml*zm)/v(im,5)    
        ENDIF   
        IF(mod(mstj(43),2).EQ.1) THEN   
          pmqth3=0.5*parj(82)   
          IF(kflgd2.EQ.22) pmqth3=0.5*parj(83)  
          pmq1=(pmth(1,kflgd1)**2+pmqth3**2)/v(i1,5)    
          pmq2=(pmth(1,kflgd2)**2+pmqth3**2)/v(i1,5)    
          zd=sqrt(max(0.,(1.-v(i1,5)/ped**2)*((1.-pmq1-pmq2)**2-    
     &    4.*pmq1*pmq2)))   
          zh=1.+pmq1-pmq2   
        ELSE    
          zd=sqrt(max(0.,1.-v(i1,5)/ped**2))    
          zh=1. 
        ENDIF   
        zl=0.5*(zh-zd)  
        zu=0.5*(zh+zd)  
        IF(i1.EQ.n+1.AND.(v(i1,1).LT.zl.OR.v(i1,1).GT.zu)) isl(1)=1 
        IF(i1.EQ.n+2.AND.(v(i1,1).LT.zl.OR.v(i1,1).GT.zu)) isl(2)=1 
        IF(kflda.EQ.21) v(i1,4)=log(zu*(1.-zl)/max(1e-20,zl*(1.-zu)))   
        IF(kflda.NE.21) v(i1,4)=log((1.-zl)/max(1e-10,1.-zu))   
  320   CONTINUE    
        IF(isl(1).EQ.1.AND.isl(2).EQ.1.AND.islm.NE.0) THEN  
          isl(3-islm)=0 
          islm=3-islm   
        ELSEIF(isl(1).EQ.1.AND.isl(2).EQ.1) THEN    
          zdr1=max(0.,v(n+1,3)/v(n+1,4)-1.) 
          zdr2=max(0.,v(n+2,3)/v(n+2,4)-1.) 
          IF(zdr2.GT.rlu(0)*(zdr1+zdr2)) isl(1)=0   
          IF(isl(1).EQ.1) isl(2)=0  
          IF(isl(1).EQ.0) islm=1    
          IF(isl(2).EQ.0) islm=2    
        ENDIF   
        IF(isl(1).EQ.1.OR.isl(2).EQ.1) goto 200 
      ENDIF 
      IF(igm.GT.0.AND.mod(mstj(43),2).EQ.1.AND.(p(n+1,5).GE.    
     &pmth(2,kfld(1)).OR.p(n+2,5).GE.pmth(2,kfld(2)))) THEN 
        pmq1=v(n+1,5)/v(im,5)   
        pmq2=v(n+2,5)/v(im,5)   
        zd=sqrt(max(0.,(1.-v(im,5)/pem**2)*((1.-pmq1-pmq2)**2-  
     &  4.*pmq1*pmq2))) 
        zh=1.+pmq1-pmq2 
        zl=0.5*(zh-zd)  
        zu=0.5*(zh+zd)  
        IF(v(im,1).LT.zl.OR.v(im,1).GT.zu) goto 200 
      ENDIF 
    
C...Accepted branch. Construct four-momentum for initial partons.   
  330 mazip=0   
      mazic=0   
      IF(nep.EQ.1) THEN 
        p(n+1,1)=0. 
        p(n+1,2)=0. 
        p(n+1,3)=sqrt(max(0.,(p(ipa(1),4)+p(n+1,5))*(p(ipa(1),4)-   
     &  p(n+1,5)))) 
        p(n+1,4)=p(ipa(1),4)    
        v(n+1,2)=p(n+1,4)   
      ELSEIF(igm.EQ.0.AND.nep.EQ.2) THEN    
        ped1=0.5*(v(im,5)+v(n+1,5)-v(n+2,5))/p(im,5)    
        p(n+1,1)=0. 
        p(n+1,2)=0. 
        p(n+1,3)=sqrt(max(0.,(ped1+p(n+1,5))*(ped1-p(n+1,5))))  
        p(n+1,4)=ped1   
        p(n+2,1)=0. 
        p(n+2,2)=0. 
        p(n+2,3)=-p(n+1,3)  
        p(n+2,4)=p(im,5)-ped1   
        v(n+1,2)=p(n+1,4)   
        v(n+2,2)=p(n+2,4)   
      ELSEIF(nep.EQ.3) THEN 
        p(n+1,1)=0. 
        p(n+1,2)=0. 
        p(n+1,3)=sqrt(max(0.,pa1s)) 
        p(n+2,1)=sqrt(pts)  
        p(n+2,2)=0. 
        p(n+2,3)=0.5*(pa3s-pa2s-pa1s)/p(n+1,3)  
        p(n+3,1)=-p(n+2,1)  
        p(n+3,2)=0. 
        p(n+3,3)=-(p(n+1,3)+p(n+2,3))   
        v(n+1,2)=p(n+1,4)   
        v(n+2,2)=p(n+2,4)   
        v(n+3,2)=p(n+3,4)   
    
C...Construct transverse momentum for ordinary branching in shower. 
      ELSE  
        zm=v(im,1)  
        pzm=sqrt(max(0.,(pem+p(im,5))*(pem-p(im,5))))   
        pmls=(v(im,5)-v(n+1,5)-v(n+2,5))**2-4.*v(n+1,5)*v(n+2,5)    
        IF(pzm.LE.0.) THEN  
          pts=0.    
        ELSEIF(mod(mstj(43),2).EQ.1) THEN   
          pts=(pem**2*(zm*(1.-zm)*v(im,5)-(1.-zm)*v(n+1,5)- 
     &    zm*v(n+2,5))-0.25*pmls)/pzm**2    
        ELSE    
          pts=pmls*(zm*(1.-zm)*pem**2/v(im,5)-0.25)/pzm**2  
        ENDIF   
        pt=sqrt(max(0.,pts))    
    
C...Find coefficient of azimuthal asymmetry due to gluon polarization.  
        hazip=0.    
        IF(mstj(49).NE.1.AND.mod(mstj(46),2).EQ.1.AND.k(im,2).EQ.21.    
     &  and.iau.NE.0) THEN  
          IF(k(igm,3).NE.0) mazip=1 
          zau=v(igm,1)  
          IF(iau.EQ.im+1) zau=1.-v(igm,1)   
          IF(mazip.EQ.0) zau=0. 
          IF(k(igm,2).NE.21) THEN   
            hazip=2.*zau/(1.+zau**2)    
          ELSE  
            hazip=(zau/(1.-zau*(1.-zau)))**2    
          ENDIF 
          IF(k(n+1,2).NE.21) THEN   
            hazip=hazip*(-2.*zm*(1.-zm))/(1.-2.*zm*(1.-zm)) 
          ELSE  
            hazip=hazip*(zm*(1.-zm)/(1.-zm*(1.-zm)))**2 
          ENDIF 
        ENDIF   
    
C...Find coefficient of azimuthal asymmetry due to soft gluon   
C...interference.   
        hazic=0.    
        IF(mstj(46).GE.2.AND.(k(n+1,2).EQ.21.OR.k(n+2,2).EQ.21).    
     &  and.iau.NE.0) THEN  
          IF(k(igm,3).NE.0) mazic=n+1   
          IF(k(igm,3).NE.0.AND.k(n+1,2).NE.21) mazic=n+2    
          IF(k(igm,3).NE.0.AND.k(n+1,2).EQ.21.AND.k(n+2,2).EQ.21.AND.   
     &    zm.GT.0.5) mazic=n+2  
          IF(k(iau,2).EQ.22) mazic=0    
          zs=zm 
          IF(mazic.EQ.n+2) zs=1.-zm 
          zgm=v(igm,1)  
          IF(iau.EQ.im-1) zgm=1.-v(igm,1)   
          IF(mazic.EQ.0) zgm=1. 
          hazic=(p(im,5)/p(igm,5))*sqrt((1.-zs)*(1.-zgm)/(zs*zgm))  
          hazic=min(0.95,hazic) 
        ENDIF   
      ENDIF 
    
C...Construct kinematics for ordinary branching in shower.  
  340 IF(nep.EQ.2.AND.igm.GT.0) THEN    
        IF(mod(mstj(43),2).EQ.1) THEN   
          p(n+1,4)=pem*v(im,1)  
        ELSE    
          p(n+1,4)=pem*(0.5*(v(im,5)-sqrt(pmls)+v(n+1,5)-v(n+2,5))+ 
     &    sqrt(pmls)*zm)/v(im,5)    
        ENDIF   
        phi=paru(2)*rlu(0)  
        p(n+1,1)=pt*cos(phi)    
        p(n+1,2)=pt*sin(phi)    
        IF(pzm.GT.0.) THEN  
          p(n+1,3)=0.5*(v(n+2,5)-v(n+1,5)-v(im,5)+2.*pem*p(n+1,4))/pzm  
        ELSE    
          p(n+1,3)=0.   
        ENDIF   
        p(n+2,1)=-p(n+1,1)  
        p(n+2,2)=-p(n+1,2)  
        p(n+2,3)=pzm-p(n+1,3)   
        p(n+2,4)=pem-p(n+1,4)   
        IF(mstj(43).LE.2) THEN  
          v(n+1,2)=(pem*p(n+1,4)-pzm*p(n+1,3))/p(im,5)  
          v(n+2,2)=(pem*p(n+2,4)-pzm*p(n+2,3))/p(im,5)  
        ENDIF   
      ENDIF 
    
C...Rotate and boost daughters. 
      IF(igm.GT.0) THEN 
        IF(mstj(43).LE.2) THEN  
          bex=p(igm,1)/p(igm,4) 
          bey=p(igm,2)/p(igm,4) 
          bez=p(igm,3)/p(igm,4) 
          ga=p(igm,4)/p(igm,5)  
          gabep=ga*(ga*(bex*p(im,1)+bey*p(im,2)+bez*p(im,3))/(1.+ga)-   
     &    p(im,4))  
        ELSE    
          bex=0.    
          bey=0.    
          bez=0.    
          ga=1. 
          gabep=0.  
        ENDIF   
        the=ulangl(p(im,3)+gabep*bez,sqrt((p(im,1)+gabep*bex)**2+   
     &  (p(im,2)+gabep*bey)**2))    
        phi=ulangl(p(im,1)+gabep*bex,p(im,2)+gabep*bey) 
        DO 350 i=n+1,n+2    
        dp(1)=cos(the)*cos(phi)*p(i,1)-sin(phi)*p(i,2)+ 
     &  sin(the)*cos(phi)*p(i,3)    
        dp(2)=cos(the)*sin(phi)*p(i,1)+cos(phi)*p(i,2)+ 
     &  sin(the)*sin(phi)*p(i,3)    
        dp(3)=-sin(the)*p(i,1)+cos(the)*p(i,3)  
        dp(4)=p(i,4)    
        dbp=bex*dp(1)+bey*dp(2)+bez*dp(3)   
        dgabp=ga*(ga*dbp/(1d0+ga)+dp(4))    
        p(i,1)=dp(1)+dgabp*bex  
        p(i,2)=dp(2)+dgabp*bey  
        p(i,3)=dp(3)+dgabp*bez  
  350   p(i,4)=ga*(dp(4)+dbp)   
      ENDIF 
    
C...Weight with azimuthal distribution, if required.    
      IF(mazip.NE.0.OR.mazic.NE.0) THEN 
        DO 360 j=1,3    
        dpt(1,j)=p(im,j)    
        dpt(2,j)=p(iau,j)   
  360   dpt(3,j)=p(n+1,j)   
        dpma=dpt(1,1)*dpt(2,1)+dpt(1,2)*dpt(2,2)+dpt(1,3)*dpt(2,3)  
        dpmd=dpt(1,1)*dpt(3,1)+dpt(1,2)*dpt(3,2)+dpt(1,3)*dpt(3,3)  
        dpmm=dpt(1,1)**2+dpt(1,2)**2+dpt(1,3)**2    
        DO 370 j=1,3    
        dpt(4,j)=dpt(2,j)-dpma*dpt(1,j)/dpmm    
  370   dpt(5,j)=dpt(3,j)-dpmd*dpt(1,j)/dpmm    
        dpt(4,4)=sqrt(dpt(4,1)**2+dpt(4,2)**2+dpt(4,3)**2)  
        dpt(5,4)=sqrt(dpt(5,1)**2+dpt(5,2)**2+dpt(5,3)**2)  
        IF(min(dpt(4,4),dpt(5,4)).GT.0.1*parj(82)) THEN 
          cad=(dpt(4,1)*dpt(5,1)+dpt(4,2)*dpt(5,2)+ 
     &    dpt(4,3)*dpt(5,3))/(dpt(4,4)*dpt(5,4))    
          IF(mazip.NE.0) THEN   
            IF(1.+hazip*(2.*cad**2-1.).LT.rlu(0)*(1.+abs(hazip)))   
     &      goto 340    
          ENDIF 
          IF(mazic.NE.0) THEN   
            IF(mazic.EQ.n+2) cad=-cad   
            IF((1.-hazic)*(1.-hazic*cad)/(1.+hazic**2-2.*hazic*cad).    
     &      lt.rlu(0)) goto 340 
          ENDIF 
        ENDIF   
      ENDIF 
    
C...Continue loop over partons that may branch, until none left.    
      IF(igm.GE.0) k(im,1)=14   
      n=n+nep   
      nep=2 
      IF(n.GT.mstu(4)-mstu(32)-5) THEN  
        CALL luerrm(11,'(LUSHOW:) no more memory left in LUJETS')   
        IF(mstu(21).GE.1) n=ns  
        IF(mstu(21).GE.1) RETURN    
      ENDIF 
      goto 140  
    
C...Set information on imagined shower initiator.   
  380 IF(npa.GE.2) THEN 
        k(ns+1,1)=11    
        k(ns+1,2)=94    
        k(ns+1,3)=ip1   
        IF(ip2.GT.0.AND.ip2.LT.ip1) k(ns+1,3)=ip2   
        k(ns+1,4)=ns+2  
        k(ns+1,5)=ns+1+npa  
        iim=1   
      ELSE  
        iim=0   
      ENDIF 
    
C...Reconstruct string drawing information. 
      DO 390 i=ns+1+iim,n   
      IF(k(i,1).LE.10.AND.k(i,2).EQ.22) THEN    
        k(i,1)=1    
      ELSEIF(k(i,1).LE.10) THEN 
        k(i,4)=mstu(5)*(k(i,4)/mstu(5)) 
        k(i,5)=mstu(5)*(k(i,5)/mstu(5)) 
      ELSEIF(k(mod(k(i,4),mstu(5))+1,2).NE.22) THEN 
        id1=mod(k(i,4),mstu(5)) 
        IF(k(i,2).GE.1.AND.k(i,2).LE.8) id1=mod(k(i,4),mstu(5))+1   
        id2=2*mod(k(i,4),mstu(5))+1-id1 
        k(i,4)=mstu(5)*(k(i,4)/mstu(5))+id1 
        k(i,5)=mstu(5)*(k(i,5)/mstu(5))+id2 
        k(id1,4)=k(id1,4)+mstu(5)*i 
        k(id1,5)=k(id1,5)+mstu(5)*id2   
        k(id2,4)=k(id2,4)+mstu(5)*id1   
        k(id2,5)=k(id2,5)+mstu(5)*i 
      ELSE  
        id1=mod(k(i,4),mstu(5)) 
        id2=id1+1   
        k(i,4)=mstu(5)*(k(i,4)/mstu(5))+id1 
        k(i,5)=mstu(5)*(k(i,5)/mstu(5))+id1 
        k(id1,4)=k(id1,4)+mstu(5)*i 
        k(id1,5)=k(id1,5)+mstu(5)*i 
        k(id2,4)=0  
        k(id2,5)=0  
      ENDIF 
  390 CONTINUE  
    
C...Transformation from CM frame.   
      IF(npa.GE.2) THEN 
        bex=ps(1)/ps(4) 
        bey=ps(2)/ps(4) 
        bez=ps(3)/ps(4) 
        ga=ps(4)/ps(5)  
        gabep=ga*(ga*(bex*p(ipa(1),1)+bey*p(ipa(1),2)+bez*p(ipa(1),3))  
     &  /(1.+ga)-p(ipa(1),4))   
      ELSE  
        bex=0.  
        bey=0.  
        bez=0.  
        gabep=0.    
      ENDIF 
      the=ulangl(p(ipa(1),3)+gabep*bez,sqrt((p(ipa(1),1)    
     &+gabep*bex)**2+(p(ipa(1),2)+gabep*bey)**2))   
      phi=ulangl(p(ipa(1),1)+gabep*bex,p(ipa(1),2)+gabep*bey)   
      IF(npa.EQ.3) THEN 
        chi=ulangl(cos(the)*cos(phi)*(p(ipa(2),1)+gabep*bex)+cos(the)*  
     &  sin(phi)*(p(ipa(2),2)+gabep*bey)-sin(the)*(p(ipa(2),3)+gabep*   
     &  bez),-sin(phi)*(p(ipa(2),1)+gabep*bex)+cos(phi)*(p(ipa(2),2)+   
     &  gabep*bey)) 
        CALL ludbrb(ns+1,n,0.,chi,0d0,0d0,0d0)  
      ENDIF 
      dbex=dble(bex)    
      dbey=dble(bey)    
      dbez=dble(bez)    
      CALL ludbrb(ns+1,n,the,phi,dbex,dbey,dbez)    
    
C...Decay vertex of shower. 
      DO 400 i=ns+1,n   
      DO 400 j=1,5  
  400 v(i,j)=v(ip1,j)   
    
C...Delete trivial shower, else connect initiators. 
      IF(n.EQ.ns+npa+iim) THEN  
        n=ns    
      ELSE  
        DO 410 ip=1,npa 
        k(ipa(ip),1)=14 
        k(ipa(ip),4)=k(ipa(ip),4)+ns+iim+ip 
        k(ipa(ip),5)=k(ipa(ip),5)+ns+iim+ip 
        k(ns+iim+ip,3)=ipa(ip)  
        IF(iim.EQ.1.AND.mstu(16).NE.2) k(ns+iim+ip,3)=ns+1  
        k(ns+iim+ip,4)=mstu(5)*ipa(ip)+k(ns+iim+ip,4)   
  410   k(ns+iim+ip,5)=mstu(5)*ipa(ip)+k(ns+iim+ip,5)   
      ENDIF 
    
      RETURN    
      END