pysspa.f

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

 
C*********************************************************************
 
C...PYSSPA
C...Generates spacelike parton showers.
 
      SUBROUTINE pysspa(IPU1,IPU2)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Commonblocks.
      common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
      common/pyint1/mint(400),vint(400)
      common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
      common/pyint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
      SAVE /pyjets/,/pydat1/,/pydat2/,/pysubs/,/pypars/,/pyint1/,
     &/pyint2/,/pyint3/
C...Local arrays and data.
      dimension kfls(4),is(2),xs(2),zs(2),q2s(2),tevcsv(2),tevesv(2),
     &xfs(2,-25:25),xfa(-25:25),xfb(-25:25),xfn(-25:25),wtapc(-25:25),
     &wtape(-25:25),wtsf(-25:25),the2(2),alam(2),dq2(3),dpc(3),dpd(4),
     &dpb(4),robo(5),more(2),kfbeam(2),q2mncs(2),kcfi(2),nfis(2),
     &thefis(2,2),isfi(2),dphi(2),mcesv(2)
      DATA is/2*0/
 
C...Read out basic information; set global Q^2 scale.
      ipus1=ipu1
      ipus2=ipu2
      isub=mint(1)
      q2mx=vint(56)
      vint2r=vint(2)*vint(143)*vint(144)
      IF(iset(isub).EQ.2.OR.iset(isub).EQ.9.OR.iset(isub).EQ.11) q2mx=
     &min(vint2r,parp(67)*vint(56))
      fcq2mx=1d0
 
C...Define which processes ME corrections have been implemented for.
      mecor=0
      IF(mstp(68).EQ.1.OR.mstp(68).EQ.3) THEN
        IF(isub.EQ.1.OR.isub.EQ.2.OR.isub.EQ.141.OR.isub.EQ.142.OR.
     &  isub.EQ.144) mecor=1
        IF(isub.EQ.102.OR.isub.EQ.152.OR.isub.EQ.157) mecor=2
        IF(isub.EQ.3.OR.isub.EQ.151.OR.isub.EQ.156) mecor=3
      ENDIF
 
C...Initialize QCD evolution and check phase space.
      q2mnc=parp(62)**2
      q2mncs(1)=q2mnc
      q2mncs(2)=q2mnc
      IF(mint(107).EQ.2.AND.mstp(66).EQ.2) THEN
        q0s=parp(15)**2
        ps=vint(3)**2
        q2eff=vint(54)*((q0s+ps)/(vint(54)+ps))*
     &  exp(ps*(vint(54)-q0s)/((vint(54)+ps)*(q0s+ps)))
        q2int=sqrt(q0s*q2eff)
        q2mncs(1)=max(q2mnc,q2int)
      ELSEIF(mint(107).EQ.3.AND.mstp(66).GE.1) THEN
        q2mncs(1)=max(q2mnc,vint(283))
      ENDIF
      IF(mint(108).EQ.2.AND.mstp(66).EQ.2) THEN
        q0s=parp(15)**2
        ps=vint(4)**2
        q2eff=vint(54)*((q0s+ps)/(vint(54)+ps))*
     &  exp(ps*(vint(54)-q0s)/((vint(54)+ps)*(q0s+ps)))
        q2int=sqrt(q0s*q2eff)
        q2mncs(2)=max(q2mnc,q2int)
      ELSEIF(mint(108).EQ.3.AND.mstp(66).GE.1) THEN
        q2mncs(2)=max(q2mnc,vint(284))
      ENDIF
      mcev=0
      alams=paru(112)
      paru(112)=parp(61)
      fq2c=1d0
      tcmx=0d0
      IF(mint(47).GE.2.AND.(mint(47).LT.5.OR.mstp(12).GE.1)) THEN
        mcev=1
        IF(mstp(64).EQ.1) fq2c=parp(63)
        IF(mstp(64).EQ.2) fq2c=parp(64)
        tcmx=log(fq2c*q2mx/parp(61)**2)
        IF(q2mx.LT.max(q2mnc,2d0*parp(61)**2).OR.tcmx.LT.0.2d0)
     &  mcev=0
      ENDIF
 
C...Initialize QED evolution and check phase space.
      meev=0
      xee=1d-10
      spme=pmas(11,1)**2
      IF(iabs(mint(11)).EQ.13.OR.iabs(mint(12)).EQ.13)
     &spme=pmas(13,1)**2
      IF(iabs(mint(11)).EQ.15.OR.iabs(mint(12)).EQ.15)
     &spme=pmas(15,1)**2
      q2mne=max(parp(68)**2,2d0*spme)
      temx=0d0
      fwte=10d0
      IF(mint(45).EQ.3.OR.mint(46).EQ.3) THEN
        meev=1
        temx=log(q2mx/spme)
        IF(q2mx.LE.q2mne.OR.temx.LT.0.2d0) meev=0
      ENDIF
      IF(mstp(61).GE.2.AND.mcev.EQ.1.AND.meev.EQ.0) THEN
        meev=2
        temx=tcmx
        fwte=1d0
      ENDIF
      IF(mcev.EQ.0.AND.meev.EQ.0) RETURN
 
C...Loopback point in case of failure to reconstruct kinematics.
      ns=n
      loop=0
      mnt352=mint(352)
      mnt353=mint(353)
      vnt352=vint(352)
      vnt353=vint(353)
  100 loop=loop+1
      IF(loop.GT.100) THEN
        mint(51)=1
        RETURN
      ENDIF
      n=ns
      mint(352)=mnt352
      mint(353)=mnt353
      vint(352)=vnt352
      vint(353)=vnt353
 
C...Initial values: flavours, momenta, virtualities.
      DO 120 jt=1,2
        more(jt)=1
        kfbeam(jt)=mint(10+jt)
        IF(mint(18+jt).EQ.1)kfbeam(jt)=22
        kfls(jt)=mint(14+jt)
        kfls(jt+2)=kfls(jt)
        xs(jt)=vint(40+jt)
        IF(mint(18+jt).EQ.1) xs(jt)=vint(40+jt)/vint(154+jt)
        IF(mint(31).GE.2) xs(jt)=xs(jt)/vint(142+jt)
        zs(jt)=1d0
        q2s(jt)=fcq2mx*q2mx
        dq2(jt)=0d0
        tevcsv(jt)=tcmx
        alam(jt)=parp(61)
        the2(jt)=1d0
        tevesv(jt)=temx
        mcesv(jt)=0
C...Calculate initial parton distribution weights.
        mint(105)=mint(102+jt)
        mint(109)=mint(106+jt)
        vint(120)=vint(2+jt)
        IF(xs(jt).LT.1d0-xee) THEN
          IF(mint(31).GE.2) mint(30)=jt
          IF(mstp(57).LE.1) THEN
            CALL pypdfu(kfbeam(jt),xs(jt),q2s(jt),xfb)
          ELSE
            CALL pypdfl(kfbeam(jt),xs(jt),q2s(jt),xfb)
          ENDIF
        ENDIF
        DO 110 kfl=-25,25
          xfs(jt,kfl)=xfb(kfl)
  110   CONTINUE
C...Special kinematics check for c/b quarks (that g -> c cbar or
C...b bbar kinematically possible).
      kflcb=iabs(kfls(jt))
      IF(kfbeam(jt).NE.22.AND.(kflcb.EQ.4.OR.kflcb.EQ.5)) THEN
        IF(xs(jt).GT.0.9d0*q2s(jt)/(pmas(kflcb,1)**2+q2s(jt))) THEN
          mint(51)=1
          RETURN
        ENDIF
      ENDIF
  120 CONTINUE
      dsh=vint(44)
      IF(iset(isub).GE.3.AND.iset(isub).LE.5) dsh=vint(26)*vint(2)
 
C...Find if interference with final state partons.
      mfis=0
      IF(mstp(67).GE.1.AND.mstp(67).LE.3) mfis=mstp(67)
      IF(mfis.NE.0) THEN
        DO 140 i=1,2
          kcfi(i)=0
          kca=pycomp(iabs(kfls(i)))
          IF(kca.NE.0) kcfi(i)=kchg(kca,2)*isign(1,kfls(i))
          nfis(i)=0
          IF(kcfi(i).NE.0) THEN
            IF(i.EQ.1) ipfs=ipus1
            IF(i.EQ.2) ipfs=ipus2
            DO 130 j=1,2
              icsi=mod(k(ipfs,3+j),mstu(5))
              IF(icsi.GT.0.AND.icsi.NE.ipus1.AND.icsi.NE.ipus2.AND.
     &        (kcfi(i).EQ.(-1)**(j+1).OR.kcfi(i).EQ.2)) THEN
                nfis(i)=nfis(i)+1
                thefis(i,nfis(i))=pyangl(p(icsi,3),sqrt(p(icsi,1)**2+
     &          p(icsi,2)**2))
                IF(i.EQ.2) thefis(i,nfis(i))=paru(1)-thefis(i,nfis(i))
              ENDIF
  130       CONTINUE
          ENDIF
  140   CONTINUE
        IF(nfis(1)+nfis(2).EQ.0) mfis=0
      ENDIF
 
C...Pick up leg with highest virtuality.
      jtold=1
  150 n=n+1
      jt=1
      IF(n.GT.ns+1.AND.q2s(2).GT.q2s(1)) jt=2
      IF(n.EQ.ns+2.AND.jt.EQ.jtold) jt=3-jt
      IF(more(jt).EQ.0) jt=3-jt
      jtold=jt
      kflb=kfls(jt)
      xb=xs(jt)
      DO 160 kfl=-25,25
        xfb(kfl)=xfs(jt,kfl)
  160 CONTINUE
      dshr=2d0*sqrt(dsh)
      dshz=dsh/zs(jt)
 
C...Check if allowed to branch.
      mcev=0
      IF(iabs(kflb).LE.10.OR.kflb.EQ.21) THEN
        mcev=1
        xec=max(parp(65)*dshr/vint2r,xb*(1d0/(1d0-parp(66))-1d0))
        IF(xb.GE.1d0-2d0*xec) mcev=0
      ENDIF
      meev=0
      IF(mint(44+jt).EQ.3) THEN
        meev=1
        IF(xb.GE.1d0-2d0*xee) meev=0
        IF((iabs(kflb).LE.10.OR.kflb.EQ.21).AND.xb.GE.1d0-2d0*xec)
     &  meev=0
C***Currently kill QED shower for resolved photoproduction.
        IF(mint(18+jt).EQ.1) meev=0
C***Currently kill shower for W inside electron.
        IF(iabs(kflb).EQ.24) THEN
          mcev=0
          meev=0
        ENDIF
      ENDIF
      IF(mstp(61).GE.2.AND.mcev.EQ.1.AND.meev.EQ.0.AND.iabs(kflb).LE.10)
     &meev=2
      IF(mcev.EQ.0.AND.meev.EQ.0) THEN
        q2b=0d0
        goto 260
      ENDIF
 
C...Maximum Q2 with or without Q2 ordering. Effective Lambda and n_f.
      q2b=q2s(jt)
      tevcb=tevcsv(jt)
      teveb=tevesv(jt)
      IF(mstp(62).LE.1) THEN
        IF(zs(jt).GT.0.99999d0) THEN
          q2b=q2s(jt)
        ELSE
          q2b=0.5d0*(1d0/zs(jt)+1d0)*q2s(jt)+0.5d0*(1d0/zs(jt)-1d0)*
     &    (q2s(3-jt)-dsh+sqrt((dsh+q2s(1)+q2s(2))**2+
     &    8d0*q2s(1)*q2s(2)*zs(jt)/(1d0-zs(jt))))
        ENDIF
        IF(mcev.EQ.1) tevcb=log(fq2c*q2b/alam(jt)**2)
        IF(meev.EQ.1) teveb=log(q2b/spme)
      ENDIF
      IF(mcev.EQ.1) THEN
        alsdum=pyalps(fq2c*q2b)
        tevcb=tevcb+2d0*log(alam(jt)/paru(117))
        alam(jt)=paru(117)
        b0=(33d0-2d0*mstu(118))/6d0
      ENDIF
      IF(meev.EQ.2) teveb=tevcb
      tevcbs=tevcb
      tevebs=teveb
 
C...Select side for interference with final state partons.
      IF(mfis.GE.1.AND.n.LE.ns+2) THEN
        ifi=n-ns
        isfi(ifi)=0
        IF(iabs(kcfi(ifi)).EQ.1.AND.nfis(ifi).EQ.1) THEN
          isfi(ifi)=1
        ELSEIF(kcfi(ifi).EQ.2.AND.nfis(ifi).EQ.1) THEN
          IF(pyr(0).GT.0.5d0) isfi(ifi)=1
        ELSEIF(kcfi(ifi).EQ.2.AND.nfis(ifi).EQ.2) THEN
          isfi(ifi)=1
          IF(pyr(0).GT.0.5d0) isfi(ifi)=2
        ENDIF
      ENDIF
 
C...Calculate preweighting factor for ME-corrected processes.
      IF(mecor.GE.1) CALL pymemx(mecor,wtff,wtgf,wtfg,wtgg)
 
C...Calculate Altarelli-Parisi weights.
      DO 170 kfl=-25,25
        wtapc(kfl)=0d0
        wtape(kfl)=0d0
        wtsf(kfl)=0d0
  170 CONTINUE
C...q -> q (g or gamma emission), g -> q.
      IF(iabs(kflb).LE.10) THEN
        wtapc(kflb)=(8d0/3d0)*log((1d0-xec-xb)*(xb+xec)/(xec*(1d0-xec)))
        wtapc(21)=0.5d0*(xb/(xb+xec)-xb/(1d0-xec))
        eq2=1d0/9d0
        IF(mod(iabs(kflb),2).EQ.0) eq2=4d0*eq2
        IF(meev.EQ.2) wtape(kflb)=2.*eq2*log((1d0-xec-xb)*(xb+xec)/
     &  (xec*(1d0-xec)))
        IF(mecor.GE.1.AND.(n.EQ.ns+1.OR.n.EQ.ns+2)) THEN
          wtapc(kflb)=wtff*wtapc(kflb)
          wtapc(21)=wtgf*wtapc(21)
          wtape(kflb)=wtff*wtape(kflb)
        ENDIF
C...f -> f, gamma -> f.
      ELSEIF(iabs(kflb).LE.20) THEN
        wtapf1=log((1d0-xee-xb)*(xb+xee)/(xee*(1d0-xee)))
        wtapf2=log((1d0-xee-xb)*(1d0-xee)/(xee*(xb+xee)))
        wtape(kflb)=2d0*(wtapf1+wtapf2)
        IF(mstp(12).GE.1) wtape(22)=xb/(xb+xee)-xb/(1d0-xee)
        IF(mecor.GE.1.AND.(n.EQ.ns+1.OR.n.EQ.ns+2)) THEN
          wtape(kflb)=wtff*wtape(kflb)
          wtape(22)=wtgf*wtape(22)
        ENDIF
C...f -> g, g -> g.
      ELSEIF(kflb.EQ.21) THEN
        wtapq=(16d0/3d0)*(sqrt((1d0-xec)/xb)-sqrt((xb+xec)/xb))
        DO 180 kfl=1,mstp(58)
          wtapc(kfl)=wtapq
          wtapc(-kfl)=wtapq
  180   CONTINUE
        wtapc(21)=6d0*log((1d0-xec-xb)/xec)
        IF(mecor.GE.1.AND.(n.EQ.ns+1.OR.n.EQ.ns+2)) THEN
          DO 190 kfl=1,mstp(58)
            wtapc(kfl)=wtfg*wtapc(kfl)
            wtapc(-kfl)=wtfg*wtapc(-kfl)
  190     CONTINUE
          wtapc(21)=wtgg*wtapc(21)
        ENDIF
C...f -> gamma, W+, W-.
      ELSEIF(kflb.EQ.22) THEN
        wtapf=log((1d0-xee-xb)*(1d0-xee)/(xee*(xb+xee)))/xb
        wtape(11)=wtapf
        wtape(-11)=wtapf
        IF(mecor.GE.1.AND.(n.EQ.ns+1.OR.n.EQ.ns+2)) THEN
          wtape(11)=wtfg*wtape(11)
          wtape(-11)=wtfg*wtape(-11)
        ENDIF
      ELSEIF(kflb.EQ.24) THEN
        wtape(-11)=1d0/(4d0*paru(102))*log((1d0-xee-xb)*(1d0-xee)/
     &  (xee*(xb+xee)))/xb
      ELSEIF(kflb.EQ.-24) THEN
        wtape(11)=1d0/(4d0*paru(102))*log((1d0-xee-xb)*(1d0-xee)/
     &  (xee*(xb+xee)))/xb
      ENDIF
 
C...Calculate parton distribution weights and sum.
      ntry=0
  200 ntry=ntry+1
      IF(ntry.GT.500) THEN
        mint(51)=1
        RETURN
      ENDIF
      wtsumc=0d0
      wtsume=0d0
      xfbo=max(1d-10,xfb(kflb))
      DO 210 kfl=-25,25
        wtsf(kfl)=xfb(kfl)/xfbo
        wtsumc=wtsumc+wtapc(kfl)*wtsf(kfl)
        wtsume=wtsume+wtape(kfl)*wtsf(kfl)
  210 CONTINUE
      wtsumc=max(0.0001d0,wtsumc)
      wtsume=max(0.0001d0/fwte,wtsume)
 
C...Choose new t: fix alpha_s, alpha_s(Q^2), alpha_s(k_T^2).
      ntry2=0
  220 ntry2=ntry2+1
      IF(ntry2.GT.500) THEN
        mint(51)=1
        RETURN
      ENDIF
      IF(mcev.EQ.1) THEN
        IF(mstp(64).LE.0) THEN
          tevcb=tevcb+log(pyr(0))*paru(2)/(paru(111)*wtsumc)
        ELSEIF(mstp(64).EQ.1) THEN
          tevcb=tevcb*exp(max(-50d0,log(pyr(0))*b0/wtsumc))
        ELSE
          tevcb=tevcb*exp(max(-50d0,log(pyr(0))*b0/(5d0*wtsumc)))
        ENDIF
      ENDIF
      IF(meev.EQ.1) THEN
        teveb=teveb*exp(max(-50d0,log(pyr(0))*paru(2)/
     &  (paru(101)*fwte*wtsume*temx)))
      ELSEIF(meev.EQ.2) THEN
        teveb=teveb+log(pyr(0))*paru(2)/(paru(101)*wtsume)
      ENDIF
 
C...Translate t into Q2 scale; choose between QCD and QED evolution.
  230 IF(mcev.EQ.1) q2cb=alam(jt)**2*exp(max(-50d0,tevcb))/fq2c
      IF(meev.EQ.1) q2eb=spme*exp(max(-50d0,teveb))
      IF(meev.EQ.2) q2eb=alam(jt)**2*exp(max(-50d0,teveb))/fq2c
C...Ensure that Q2 is above threshold for charm/bottom.
      kflcb=iabs(kflb)
      IF(kfbeam(jt).NE.22.AND.(kflcb.EQ.4.OR.kflcb.EQ.5).AND.
     &mcev.EQ.1) THEN
        IF(q2cb.LT.pmas(kflcb,1)**2) THEN
          q2cb=1.1d0*pmas(kflcb,1)**2
          tevcb=log(fq2c*q2b/alam(jt)**2)
          fcq2mx=min(2d0,1.05d0*fcq2mx)
        ENDIF
      ENDIF
      IF(kfbeam(jt).NE.22.AND.(kflcb.EQ.4.OR.kflcb.EQ.5).AND.
     &meev.EQ.2) THEN
        IF(q2eb.LT.pmas(kflcb,1)**2) meev=0
      ENDIF
      mce=0
      IF(mcev.EQ.0.AND.meev.EQ.0) THEN
      ELSEIF(mcev.EQ.1.AND.meev.EQ.0) THEN
        IF(q2cb.GT.q2mncs(jt)) mce=1
      ELSEIF(mcev.EQ.0.AND.meev.EQ.1) THEN
        IF(q2eb.GT.q2mne) mce=2
      ELSEIF(mcev.EQ.0.AND.meev.EQ.2) THEN
        IF(q2eb.GT.q2mncs(jt)) mce=2
      ELSEIF(mcev.EQ.1.AND.meev.EQ.2) THEN
        IF(q2cb.GT.q2eb.AND.q2cb.GT.q2mncs(jt)) mce=1
        IF(q2eb.GT.q2cb.AND.q2eb.GT.q2mncs(jt)) mce=2
      ELSEIF(q2mncs(jt).GT.q2mne) THEN
        mce=1
        IF(q2eb.GT.q2cb.OR.q2cb.LE.q2mncs(jt)) mce=2
        IF(mce.EQ.2.AND.q2eb.LE.q2mne) mce=0
      ELSE
        mce=2
        IF(q2cb.GT.q2eb.OR.q2eb.LE.q2mne) mce=1
        IF(mce.EQ.1.AND.q2cb.LE.q2mncs(jt)) mce=0
      ENDIF
 
C...Evolution possibly ended. Update t values.
      IF(mce.EQ.0) THEN
        q2b=0d0
        goto 260
      ELSEIF(mce.EQ.1) THEN
        q2b=q2cb
        q2ref=fq2c*q2b
        IF(meev.EQ.1) teveb=log(q2b/spme)
        IF(meev.EQ.2) teveb=log(fq2c*q2b/alam(jt)**2)
      ELSE
        q2b=q2eb
        q2ref=q2b
        IF(mcev.EQ.1) tevcb=log(fq2c*q2b/alam(jt)**2)
      ENDIF
 
C...Select flavour for branching parton.
      IF(mce.EQ.1) wtran=pyr(0)*wtsumc
      IF(mce.EQ.2) wtran=pyr(0)*wtsume
      kfla=-25
  240 kfla=kfla+1
      IF(mce.EQ.1) wtran=wtran-wtapc(kfla)*wtsf(kfla)
      IF(mce.EQ.2) wtran=wtran-wtape(kfla)*wtsf(kfla)
      IF(kfla.LE.24.AND.wtran.GT.0d0) goto 240
      IF(kfla.EQ.25) THEN
        q2b=0d0
        goto 260
      ENDIF
 
C...Choose z value and corrective weight.
      wtz=0d0
C...q -> q + g or q -> q + gamma.
      IF(iabs(kfla).LE.10.AND.iabs(kflb).LE.10) THEN
        z=1d0-((1d0-xb-xec)/(1d0-xec))*
     &  (xec*(1d0-xec)/((xb+xec)*(1d0-xb-xec)))**pyr(0)
        wtz=0.5d0*(1d0+z**2)
C...q -> g + q.
      ELSEIF(iabs(kfla).LE.10.AND.kflb.EQ.21) THEN
        z=xb/(sqrt(xb+xec)+pyr(0)*(sqrt(1d0-xec)-sqrt(xb+xec)))**2
        wtz=0.5d0*(1d0+(1d0-z)**2)*sqrt(z)
C...f -> f + gamma.
      ELSEIF(iabs(kfla).LE.20.AND.iabs(kflb).LE.20) THEN
        IF(wtapf1.GT.pyr(0)*(wtapf1+wtapf2)) THEN
          z=1d0-((1d0-xb-xee)/(1d0-xee))*
     &    (xee*(1d0-xee)/((xb+xee)*(1d0-xb-xee)))**pyr(0)
        ELSE
          z=xb+xb*(xee/(1d0-xee))*
     &    ((1d0-xb-xee)*(1d0-xee)/(xee*(xb+xee)))**pyr(0)
        ENDIF
        wtz=0.5d0*(1d0+z**2)*(z-xb)/(1d0-xb)
C...f -> gamma + f.
      ELSEIF(iabs(kfla).LE.20.AND.kflb.EQ.22) THEN
        z=xb+xb*(xee/(1d0-xee))*
     &  ((1d0-xb-xee)*(1d0-xee)/(xee*(xb+xee)))**pyr(0)
        wtz=0.5d0*(1d0+(1d0-z)**2)*xb*(z-xb)/z
C...f -> W+- + f.
      ELSEIF(iabs(kfla).LE.20.AND.iabs(kflb).EQ.24) THEN
        z=xb+xb*(xee/(1d0-xee))*
     &  ((1d0-xb-xee)*(1d0-xee)/(xee*(xb+xee)))**pyr(0)
        wtz=0.5d0*(1d0+(1d0-z)**2)*(xb*(z-xb)/z)*
     &  (q2b/(q2b+pmas(24,1)**2))
C...g -> q + qbar.
      ELSEIF(kfla.EQ.21.AND.iabs(kflb).LE.10) THEN
        z=xb/(1d0-xec)+pyr(0)*(xb/(xb+xec)-xb/(1d0-xec))
        wtz=1d0-2d0*z*(1d0-z)
C...g -> g + g.
      ELSEIF(kfla.EQ.21.AND.kflb.EQ.21) THEN
        z=1d0/(1d0+((1d0-xec-xb)/xb)*(xec/(1d0-xec-xb))**pyr(0))
        wtz=(1d0-z*(1d0-z))**2
C...gamma -> f + fbar.
      ELSEIF(kfla.EQ.22.AND.iabs(kflb).LE.20) THEN
        z=xb/(1d0-xee)+pyr(0)*(xb/(xb+xee)-xb/(1d0-xee))
        wtz=1d0-2d0*z*(1d0-z)
      ENDIF
      IF(mce.EQ.2.AND.meev.EQ.1) wtz=(wtz/fwte)*(teveb/temx)
 
C...Option with resummation of soft gluon emission as effective z shift.
      IF(mce.EQ.1) THEN
        IF(mstp(65).GE.1) THEN
          rsoft=6d0
          IF(kflb.NE.21) rsoft=8d0/3d0
          z=z*(tevcb/tevcsv(jt))**(rsoft*xec/((xb+xec)*b0))
          IF(z.LE.xb) goto 220
        ENDIF
 
C...Option with alpha_s(k_T^2): demand k_T^2 > cutoff, reweight.
        IF(mstp(64).GE.2) THEN
          IF((1d0-z)*q2b.LT.q2mncs(jt)) goto 220
          alprat=tevcb/(tevcb+log(1d0-z))
          IF(alprat.LT.5d0*pyr(0)) goto 220
          IF(alprat.GT.5d0) wtz=wtz*alprat/5d0
        ENDIF
      ENDIF
 
C...Remove kinematically impossible branchings.
      uhat=q2b-dsh*(1d0-z)/z
      IF(mstp(68).GE.0.AND.uhat.GT.0d0) goto 220
 
C...Select phi angle of branching at random.
      phibr=paru(2)*pyr(0)
 
C...Matrix-element corrections for some processes.
      IF(mecor.GE.1.AND.(n.EQ.ns+1.OR.n.EQ.ns+2)) THEN
        IF(iabs(kfla).LE.20.AND.iabs(kflb).LE.20) THEN
          CALL pymewt(mecor,1,q2b,z,phibr,wtme)
          wtz=wtz*wtme/wtff
        ELSEIF((kfla.EQ.21.OR.kfla.EQ.22).AND.iabs(kflb).LE.20) THEN
          CALL pymewt(mecor,2,q2b,z,phibr,wtme)
          wtz=wtz*wtme/wtgf
        ELSEIF(iabs(kfla).LE.20.AND.(kflb.EQ.21.OR.kflb.EQ.22)) THEN
          CALL pymewt(mecor,3,q2b,z,phibr,wtme)
          wtz=wtz*wtme/wtfg
        ELSEIF(kfla.EQ.21.AND.kflb.EQ.21) THEN
          CALL pymewt(mecor,4,q2b,z,phibr,wtme)
          wtz=wtz*wtme/wtgg
        ENDIF
      ENDIF
 
C...Impose angular constraint in first branching from interference
C...with final state partons.
      IF(mce.EQ.1) THEN
        IF(mfis.GE.1.AND.n.LE.ns+2.AND.ntry2.LT.200) THEN
          the2d=(4d0*q2b)/(dsh*(1d0-z))
          IF(n.EQ.ns+1.AND.isfi(1).GE.1) THEN
            IF(the2d.GT.thefis(1,isfi(1))**2) goto 220
          ELSEIF(n.EQ.ns+2.AND.isfi(2).GE.1) THEN
            IF(the2d.GT.thefis(2,isfi(2))**2) goto 220
          ENDIF
        ENDIF
 
C...Option with angular ordering requirement.
        IF(mstp(62).GE.3.AND.ntry2.LT.200) THEN
          the2t=(4d0*z**2*q2b)/(4d0*z**2*q2b+(1d0-z)*xb**2*vint2r)
          IF(the2t.GT.the2(jt)) goto 220
        ENDIF
      ENDIF
 
C...Weighting with new parton distributions.
      mint(105)=mint(102+jt)
      mint(109)=mint(106+jt)
      vint(120)=vint(2+jt)
      IF(mint(31).GE.2) mint(30)=jt
      IF(mstp(57).LE.1) THEN
        CALL pypdfu(kfbeam(jt),xb,q2ref,xfn)
      ELSE
        CALL pypdfl(kfbeam(jt),xb,q2ref,xfn)
      ENDIF
      xfbn=xfn(kflb)
      IF(xfbn.LT.1d-20) THEN
        IF(kfla.EQ.kflb) THEN
          tevcb=tevcbs
          teveb=tevebs
          wtapc(kflb)=0d0
          wtape(kflb)=0d0
          goto 200
        ELSEIF(mce.EQ.1.AND.tevcbs-tevcb.GT.0.2d0) THEN
          tevcb=0.5d0*(tevcbs+tevcb)
          goto 230
        ELSEIF(mce.EQ.2.AND.tevebs-teveb.GT.0.2d0) THEN
          teveb=0.5d0*(tevebs+teveb)
          goto 230
        ELSE
          xfbn=1d-10
          xfn(kflb)=xfbn
        ENDIF
      ENDIF
      DO 250 kfl=-25,25
        xfb(kfl)=xfn(kfl)
  250 CONTINUE
      xa=xb/z
      IF(mint(31).GE.2) mint(30)=jt
      IF(mstp(57).LE.1) THEN
        CALL pypdfu(kfbeam(jt),xa,q2ref,xfa)
      ELSE
        CALL pypdfl(kfbeam(jt),xa,q2ref,xfa)
      ENDIF
      xfan=xfa(kfla)
      IF(xfan.LT.1d-20) goto 200
      wtsfa=wtsf(kfla)
      IF(wtz*xfan/xfbn.LT.pyr(0)*wtsfa) goto 200
 
C...Define two hard scatterers in their CM-frame.
  260 IF(n.EQ.ns+2) THEN
        dq2(jt)=q2b
        dplcm=sqrt((dsh+dq2(1)+dq2(2))**2-4d0*dq2(1)*dq2(2))/dshr
        DO 280 jr=1,2
          i=ns+jr
          IF(jr.EQ.1) ipo=ipus1
          IF(jr.EQ.2) ipo=ipus2
          DO 270 j=1,5
            k(i,j)=0
            p(i,j)=0d0
            v(i,j)=0d0
  270     CONTINUE
          k(i,1)=14
          k(i,2)=kfls(jr+2)
          k(i,4)=ipo
          k(i,5)=ipo
          p(i,3)=dplcm*(-1)**(jr+1)
          p(i,4)=(dsh+dq2(3-jr)-dq2(jr))/dshr
          p(i,5)=-sqrt(dq2(jr))
          k(ipo,1)=14
          k(ipo,3)=i
          k(ipo,4)=mod(k(ipo,4),mstu(5))+mstu(5)*i
          k(ipo,5)=mod(k(ipo,5),mstu(5))+mstu(5)*i
  280   CONTINUE
 
C...Find maximum allowed mass of timelike parton.
      ELSEIF(n.GT.ns+2) THEN
        jr=3-jt
        dq2(3)=q2b
        dpc(1)=p(is(1),4)
        dpc(2)=p(is(2),4)
        dpc(3)=0.5d0*(abs(p(is(1),3))+abs(p(is(2),3)))
        dpd(1)=dsh+dq2(jr)+dq2(jt)
        dpd(2)=dshz+dq2(jr)+dq2(3)
        dpd(3)=sqrt(dpd(1)**2-4d0*dq2(jr)*dq2(jt))
        dpd(4)=sqrt(dpd(2)**2-4d0*dq2(jr)*dq2(3))
        ikin=0
        IF(q2s(jr).GE.0.25d0*q2mnc.AND.dpd(1)-dpd(3).GE.
     &  1d-10*dpd(1)) ikin=1
        IF(ikin.EQ.0) dmsma=(dq2(jt)/zs(jt)-dq2(3))*
     &  (dsh/(dsh+dq2(jt))-dsh/(dshz+dq2(3)))
        IF(ikin.EQ.1) dmsma=(dpd(1)*dpd(2)-dpd(3)*dpd(4))/
     &  (2d0*dq2(jr))-dq2(jt)-dq2(3)
 
C...Generate timelike parton shower (if required).
        it=n
        DO 290 j=1,5
          k(it,j)=0
          p(it,j)=0d0
          v(it,j)=0d0
  290   CONTINUE
C...f -> f + g (gamma).
        IF(iabs(kflb).LE.20.AND.iabs(kfls(jt+2)).LE.20) THEN
          k(it,2)=21
          IF(mcesv(jt).EQ.2.OR.iabs(kflb).GE.11) k(it,2)=22
C...f -> g (gamma, W+-) + f.
        ELSEIF(iabs(kflb).LE.20.AND.iabs(kfls(jt+2)).GT.20) THEN
          k(it,2)=kflb
          IF(kfls(jt+2).EQ.24) THEN
            k(it,2)=-12
          ELSEIF(kfls(jt+2).EQ.-24) THEN
            k(it,2)=12
          ENDIF
C...g (gamma) -> f + fbar, g + g.
        ELSE
          k(it,2)=-kfls(jt+2)
          IF(kfls(jt+2).GT.20) k(it,2)=kfls(jt+2)
        ENDIF
        k(it,1)=3
        IF((iabs(k(it,2)).GE.11.AND.iabs(k(it,2)).LE.18).OR.
     &  iabs(k(it,2)).EQ.22) k(it,1)=1
        p(it,5)=pymass(k(it,2))
        IF(dmsma.LE.p(it,5)**2) goto 100
        IF(mstp(63).GE.1.AND.mcesv(jt).EQ.1) THEN
          mstj48=mstj(48)
          parj85=parj(85)
          p(it,4)=(dshz-dsh-p(it,5)**2)/dshr
          p(it,3)=sqrt(p(it,4)**2-p(it,5)**2)
          IF(mstp(63).EQ.1) THEN
            q2tim=dmsma
          ELSEIF(mstp(63).EQ.2) THEN
            q2tim=min(dmsma,parp(71)*q2s(jt))
          ELSE
            q2tim=dmsma
            mstj(48)=1
            IF(ikin.EQ.0) dpt2=dmsma*(dshz+dq2(3))/(dsh+dq2(jt))
            IF(ikin.EQ.1) dpt2=dmsma*(0.5d0*dpd(1)*dpd(2)+0.5d0*dpd(3)*
     &      dpd(4)-dq2(jr)*(dq2(jt)+dq2(3)))/(4d0*dsh*dpc(3)**2)
            parj(85)=sqrt(max(0d0,dpt2))*
     &      (1d0/p(it,4)+1d0/p(is(jt),4))
          ENDIF
          CALL pyshow(it,0,sqrt(q2tim))
          mstj(48)=mstj48
          parj(85)=parj85
          IF(n.GE.it+1) p(it,5)=p(it+1,5)
        ENDIF
 
C...Reconstruct kinematics of branching: timelike parton shower.
        dms=p(it,5)**2
        IF(ikin.EQ.0) dpt2=(dmsma-dms)*(dshz+dq2(3))/(dsh+dq2(jt))
        IF(ikin.EQ.1) dpt2=(dmsma-dms)*(0.5d0*dpd(1)*dpd(2)+
     &  0.5d0*dpd(3)*dpd(4)-dq2(jr)*(dq2(jt)+dq2(3)+dms))/
     &  (4d0*dsh*dpc(3)**2)
        IF(dpt2.LT.0d0) goto 100
        dpb(1)=(0.5d0*dpd(2)-dpc(jr)*(dshz+dq2(jr)-dq2(jt)-dms)/
     &  dshr)/dpc(3)-dpc(3)
        p(it,1)=sqrt(dpt2)
        p(it,3)=dpb(1)*(-1)**(jt+1)
        p(it,4)=sqrt(dpt2+dpb(1)**2+dms)
        IF(n.GE.it+1) THEN
          dpb(1)=sqrt(dpb(1)**2+dpt2)
          dpb(2)=sqrt(dpb(1)**2+dms)
          dpb(3)=p(it+1,3)
          dpb(4)=sqrt(dpb(3)**2+dms)
          dbez=(dpb(4)*dpb(1)-dpb(3)*dpb(2))/(dpb(4)*dpb(2)-dpb(3)*
     &    dpb(1))
          CALL pyrobo(it+1,n,0d0,0d0,0d0,0d0,dbez)
          the=pyangl(p(it,3),p(it,1))
          CALL pyrobo(it+1,n,the,0d0,0d0,0d0,0d0)
        ENDIF
 
C...Reconstruct kinematics of branching: spacelike parton.
        DO 300 j=1,5
          k(n+1,j)=0
          p(n+1,j)=0d0
          v(n+1,j)=0d0
  300   CONTINUE
        k(n+1,1)=14
        k(n+1,2)=kflb
        p(n+1,1)=p(it,1)
        p(n+1,3)=p(it,3)+p(is(jt),3)
        p(n+1,4)=p(it,4)+p(is(jt),4)
        p(n+1,5)=-sqrt(dq2(3))
 
C...Define colour flow of branching.
        k(is(jt),3)=n+1
        k(it,3)=n+1
        im1=n+1
        im2=n+1
C...f -> f + gamma (Z, W).
        IF(iabs(k(it,2)).GE.22) THEN
          k(it,1)=1
          id1=is(jt)
          id2=is(jt)
C...f -> gamma (Z, W) + f.
        ELSEIF(iabs(k(is(jt),2)).GE.22) THEN
          id1=it
          id2=it
C...gamma -> q + qbar, g + g.
        ELSEIF(k(n+1,2).EQ.22) THEN
          id1=is(jt)
          id2=it
          im1=id2
          im2=id1
C...q -> q + g.
        ELSEIF(k(n+1,2).GT.0.AND.k(n+1,2).NE.21.AND.k(it,2).EQ.21) THEN
          id1=it
          id2=is(jt)
C...q -> g + q.
        ELSEIF(k(n+1,2).GT.0.AND.k(n+1,2).NE.21) THEN
          id1=is(jt)
          id2=it
C...qbar -> qbar + g.
        ELSEIF(k(n+1,2).LT.0.AND.k(it,2).EQ.21) THEN
          id1=is(jt)
          id2=it
C...qbar -> g + qbar.
        ELSEIF(k(n+1,2).LT.0) THEN
          id1=it
          id2=is(jt)
C...g -> g + g; g -> q + qbar.
        ELSEIF((k(it,2).EQ.21.AND.pyr(0).GT.0.5d0).OR.k(it,2).LT.0) THEN
          id1=is(jt)
          id2=it
        ELSE
          id1=it
          id2=is(jt)
        ENDIF
        IF(im1.EQ.n+1) k(im1,4)=k(im1,4)+id1
        IF(im2.EQ.n+1) k(im2,5)=k(im2,5)+id2
        k(id1,4)=k(id1,4)+mstu(5)*im1
        k(id2,5)=k(id2,5)+mstu(5)*im2
        IF(id1.NE.id2) THEN
          k(id1,5)=k(id1,5)+mstu(5)*id2
          k(id2,4)=k(id2,4)+mstu(5)*id1
        ENDIF
        n=n+1
        IF(k(it,1).EQ.1) THEN
          k(it,4)=0
          k(it,5)=0
        ENDIF
 
C...Boost to new CM-frame.
        dbsvx=(p(n,1)+p(is(jr),1))/(p(n,4)+p(is(jr),4))
        dbsvz=(p(n,3)+p(is(jr),3))/(p(n,4)+p(is(jr),4))
        IF(dbsvx**2+dbsvz**2.GE.1d0) goto 100
        CALL pyrobo(ns+1,n,0d0,0d0,-dbsvx,0d0,-dbsvz)
        ir=n+(jt-1)*(is(1)-n)
        CALL pyrobo(ns+1,n,-pyangl(p(ir,3),p(ir,1)),dphi(jt),
     &  0d0,0d0,0d0)
 
C...Global statistics.
        mint(352)=mint(352)+1
        vint(352)=vint(352)+sqrt(p(it,1)**2+p(it,2)**2)
        IF (mint(352).EQ.1) vint(357)=sqrt(p(it,1)**2+p(it,2)**2)
      ENDIF
 
C...Update kinematics variables.
      is(jt)=n
      dq2(jt)=q2b
      IF(mstp(62).GE.3.AND.ntry2.LT.200) the2(jt)=the2t
      dsh=dshz
 
C...Save quantities; loop back.
      q2s(jt)=q2b
      dphi(jt)=phibr
      mcesv(jt)=mce
      IF((mcev.EQ.1.AND.q2b.GE.0.25d0*q2mnc).OR.
     &(meev.EQ.1.AND.q2b.GE.q2mne)) THEN
        kfls(jt+2)=kfls(jt)
        kfls(jt)=kfla
        xs(jt)=xa
        zs(jt)=z
        DO 310 kfl=-25,25
          xfs(jt,kfl)=xfa(kfl)
  310   CONTINUE
        tevcsv(jt)=tevcb
        tevesv(jt)=teveb
      ELSE
        more(jt)=0
        IF(jt.EQ.1) ipu1=n
        IF(jt.EQ.2) ipu2=n
      ENDIF
      IF(n.GT.mstu(4)-mstu(32)-10) THEN
        CALL pyerrm(11,'(PYSSPA:) no more memory left in PYJETS')
        IF(mstu(21).GE.1) n=ns
        IF(mstu(21).GE.1) RETURN
      ENDIF
      IF(more(1).EQ.1.OR.more(2).EQ.1) goto 150
 
C...Boost hard scattering partons to frame of shower initiators.
      DO 320 j=1,3
        robo(j+2)=(p(ns+1,j)+p(ns+2,j))/(p(ns+1,4)+p(ns+2,4))
  320 CONTINUE
      k(n+2,1)=1
      DO 330 j=1,5
        p(n+2,j)=p(ns+1,j)
  330 CONTINUE
      CALL pyrobo(n+2,n+2,0d0,0d0,-robo(3),-robo(4),-robo(5))
      robo(2)=pyangl(p(n+2,1),p(n+2,2))
      robo(1)=pyangl(p(n+2,3),sqrt(p(n+2,1)**2+p(n+2,2)**2))
      imin=mint(83)+5
      IF(mint(31).GE.2) imin=min(ipus1,ipus2)
      CALL pyrobo(imin,ns,0d0,-robo(2),0d0,0d0,0d0)
      CALL pyrobo(imin,ns,robo(1),robo(2),robo(3),robo(4),robo(5))
 
C...Store user information. Reset Lambda value.
      IF(mint(31).LE.1) THEN
        k(ipu1,3)=mint(83)+3
        k(ipu2,3)=mint(83)+4
      ELSE
        k(ipu1,3)=mint(83)+1
        k(ipu2,3)=mint(83)+2
      ENDIF
      DO 340 jt=1,2
        mint(12+jt)=kfls(jt)
        vint(140+jt)=xs(jt)
        IF(mint(18+jt).EQ.1) vint(140+jt)=vint(154+jt)*xs(jt)
        IF(mint(31).GE.2) vint(140+jt)=vint(140+jt)*vint(142+jt)
  340 CONTINUE
      paru(112)=alams
 
      RETURN
      END