pyptfs.f

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

 
C*********************************************************************
 
C...PYPTFS
C...Generates pT-ordered timelike final-state parton showers.
 
C...MODE defines how to find radiators and recoilers.
C... = 0 : based on colour flow between undecayed partons.
C... = 1 : for IPART <= NPARTD only consider primary partons,
C...       whether decayed or not; else as above.
C... = 2 : based on common history, whether decayed or not.
 
      SUBROUTINE pyptfs(MODE,PTMAX,PTMIN,PTGEN)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
C...Parameter statement to help give large particle numbers.
      parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
     &kexcit=4000000,kdimen=5000000)
C...Parameter statement for maximum size of showers.
      parameter(maxnur=1000)
C...Commonblocks.
      common/pypart/npart,npartd,ipart(maxnur),ptpart(maxnur)
      common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
      common/pyctag/nct,mct(4000,2)
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      common/pypars/mstp(200),parp(200),msti(200),pari(200)
      common/pyint1/mint(400),vint(400)
      SAVE /pypart/,/pyjets/,/pyctag/,/pydat1/,/pydat2/,/pypars/,
     &/pyint1/
C...Local arrays.
      dimension ipos(2*maxnur),irec(2*maxnur),iflg(2*maxnur),
     &iscol(2*maxnur),ischg(2*maxnur),ptsca(2*maxnur),imesav(2*maxnur),
     &pt2sav(2*maxnur),zsav(2*maxnur),shtsav(2*maxnur),
     &mesys(maxnur,0:2),psum(5),dpt(5,4)
C...Statement functions.
      shat(i,j)=(p(i,4)+p(j,4))**2-(p(i,1)+p(j,1))**2-
     &(p(i,2)+p(j,2))**2-(p(i,3)+p(j,3))**2
 
C...Initial values. Check that valid system.
      ptgen=0d0
      IF(mstj(41).NE.1.AND.mstj(41).NE.2.AND.mstj(41).NE.11.AND.
     &mstj(41).NE.12) RETURN
      IF(npart.LE.0) THEN
        CALL pyerrm(2,'(PYPTFS:) showering system too small')
        RETURN
      ENDIF
      pt2cmx=ptmax**2
 
C...Mass thresholds and Lambda for QCD evolution.
      pmb=pmas(5,1)
      pmc=pmas(4,1)
      alam5=parj(81)
      alam4=alam5*(pmb/alam5)**(2d0/25d0)
      alam3=alam4*(pmc/alam4)**(2d0/27d0)
      pmbs=pmb**2
      pmcs=pmc**2
      alam5s=alam5**2
      alam4s=alam4**2
      alam3s=alam3**2
 
C...Cutoff scale for QCD evolution. Starting pT2.
      nflav=max(0,min(5,mstj(45)))
      pt0c=0.5d0*parj(82)
      pt2cmn=max(ptmin,pt0c,1.1d0*alam3)**2
 
C...Parameters for QED evolution.
      aem2pi=paru(101)/paru(2)
      pt0eq=0.5d0*parj(83)
      pt0el=0.5d0*parj(90)

C...Reset. Remove irrelevant colour tags.
      nevol=0
      DO 100 j=1,4
        psum(j)=0d0
  100 CONTINUE
      DO 110 i=mint(84)+1,n
        IF(k(i,2).GT.0.AND.k(i,2).LT.6) THEN
          k(i,5)=0
          mct(i,2)=0
        ENDIF
        IF(k(i,2).LT.0.AND.k(i,2).GT.-6) THEN
          k(i,4)=0
          mct(i,1)=0
        ENDIF
  110 CONTINUE
      nparts=npart
 
C...Begin loop to set up showering partons. Sum four-momenta.
      DO 210 ip=1,npart
        i=ipart(ip)
        IF(mode.NE.1.OR.i.GT.npartd) THEN
          IF(k(i,1).GT.10) goto 210
        ELSEIF(k(i,3).GT.mint(84)) THEN
          IF(k(i,3).GT.mint(84)+2) goto 210
        ELSE
          IF(k(k(i,3),3).GT.mint(83)+6) goto 210
        ENDIF
        DO 120 j=1,4
          psum(j)=psum(j)+p(i,j)
  120   CONTINUE
 
C...Find colour and charge, but skip diquarks.
        IF(iabs(k(i,2)).GT.1000.AND.iabs(k(i,2)).LT.10000) goto 210
        kcol=isign(kchg(pycomp(k(i,2)),2),k(i,2))
        kcha=isign(kchg(pycomp(k(i,2)),1),k(i,2))
 
C...Either colour or anticolour charge radiates; for gluon both.
        DO 160 jsgcol=1,-1,-2
          IF(kcol.EQ.jsgcol.OR.kcol.EQ.2) THEN
            jcol=4+(1-jsgcol)/2
            jcolr=9-jcol
 
C...Basic info about radiating parton.
            nevol=nevol+1
            ipos(nevol)=i
            iflg(nevol)=0
            iscol(nevol)=jsgcol
            ischg(nevol)=0
            ptsca(nevol)=ptpart(ip)
 
C...Begin search for colour recoiler when MODE = 0 or 1.
            IF(mode.LE.1) THEN
C...Find sister with matching anticolour to the radiating parton.
              irold=i
              irnew=k(irold,jcol)/mstu(5)
              move=1
 
C...The following will add MCT colour tracing for unprepped events
C...If not done, trace Les Houches colour tags for this dipole
C              IF (MCT(I,JCOL-3).EQ.0) THEN 
C                CALL PYCTTR(I,JCOL,INEW)
C...Clean up mother/daughter 'read' tags set by PYCTTR
C                DO 125 IR=1,N
C                  K(IR,4)=MOD(K(IR,4),MSTU(5)**2)
C                  K(IR,5)=MOD(K(IR,5),MSTU(5)**2)
C 125            CONTINUE
C              ENDIF

C...Skip radiation off loose colour ends.
  130         IF(irnew.EQ.0) THEN
                nevol=nevol-1
                goto 160
 
C...Optionally skip radiation on dipole to beam remnant.
              ELSEIF(mstp(72).LE.1.AND.irnew.GT.mint(53)) THEN
                nevol=nevol-1
                goto 160
 
C...For now always skip radiation on dipole to junction.
              ELSEIF(k(irnew,2).EQ.88) THEN
                nevol=nevol-1
                goto 160
 
C...For MODE=1: if reached primary then done.
              ELSEIF(mode.EQ.1.AND.irnew.GT.mint(84)+2.AND.
     &        irnew.LE.npartd) THEN
 
C...If sister stable and points back then done.
              ELSEIF(move.EQ.1.AND.k(irnew,jcolr)/mstu(5).EQ.irold)
     &        THEN
                IF(k(irnew,1).LT.10) THEN
 
C...If sister unstable then go to her daughter.
                ELSE
                  irold=irnew
                  irnew=mod(k(irnew,jcolr),mstu(5))
                  move=2
                  goto 130
               ENDIF
 
C...If found mother then look for aunt.
              ELSEIF(move.EQ.1.AND.mod(k(irnew,jcol),mstu(5)).EQ.
     &        irold) THEN
                irold=irnew
                irnew=k(irold,jcol)/mstu(5)
                goto 130
 
C...If daughter stable then done.
              ELSEIF(move.EQ.2.AND.k(irnew,jcolr)/mstu(5).EQ.irold)
     &        THEN
                IF(k(irnew,1).LT.10) THEN
 
C...If daughter unstable then go to granddaughter.
                ELSE
                  irold=irnew
                  irnew=mod(k(irnew,jcolr),mstu(5))
                  move=2
                  goto 130
                ENDIF
 
C...If daughter points to another daughter then done or move up.
              ELSEIF(move.EQ.2.AND.mod(k(irnew,jcol),mstu(5)).EQ.
     &        irold) THEN
                IF(k(irnew,1).LT.10) THEN
                ELSE
                  irold=irnew
                  irnew=k(irnew,jcol)/mstu(5)
                  move=1
                  goto 130
                ENDIF
              ENDIF
 
C...Begin search for colour recoiler when MODE = 2.
            ELSE
              irold=i
              irnew=k(irold,jcol)/mstu(5)
  140         IF(k(irnew,jcolr)/mstu(5).NE.irold) THEN
C...Step up to mother if radiating parton already branched.
                IF(k(irnew,2).EQ.k(irold,2)) THEN
                  irold=irnew
                  irnew=k(irold,jcol)/mstu(5)
                  goto 140
C...Pick sister by history if no anticolour available.
                ELSE
                  IF(irold.GT.1.AND.k(irold-1,3).EQ.k(irold,3)) THEN
                    irnew=irold-1
                  ELSEIF(irold.LT.n.AND.k(irold+1,3).EQ.k(irold,3))
     &            THEN
                    irnew=irold+1
C...Last resort: pick at random among other primaries.
                  ELSE
                    istep=max(1,min(npart-1,int(1d0+(npart-1)*pyr(0))))
                    irnew=ipart(1+mod(ip+istep-1,npart))
                  ENDIF
                ENDIF
              ENDIF
C...Trace down if sister branched.
  150         IF(k(irnew,1).GT.10) THEN
                irnew=mod(k(irnew,jcolr),mstu(5))
                goto 150
              ENDIF
            ENDIF
 
C...Now found other end of colour dipole.
            irec(nevol)=irnew
          ENDIF
  160   CONTINUE
 
C...Also electrical charge may radiate; so far only quarks and leptons.
        IF((mstj(41).EQ.2.OR.mstj(41).EQ.12).AND.kcha.NE.0.AND.
     &  iabs(k(i,2)).LE.18) THEN
 
C...Basic info about radiating parton.
          nevol=nevol+1
          ipos(nevol)=i
          iflg(nevol)=0
          iscol(nevol)=0
          ischg(nevol)=kcha
          ptsca(nevol)=ptpart(ip)
 
C...Pick nearest (= smallest invariant mass) charged particle
C...as recoiler when MODE = 0 or 1 (but for latter among primaries).
          IF(mode.LE.1) THEN
            irnew=0
            pm2min=vint(2)
            DO 170 ip2=1,npart+n-mint(53)
              IF(ip2.EQ.ip) goto 170
              IF(ip2.LE.npart) THEN
                i2=ipart(ip2)
                IF(mode.NE.1.OR.i2.GT.npartd) THEN
                  IF(k(i2,1).GT.10) goto 170
                ELSEIF(k(i2,3).GT.mint(84)) THEN
                  IF(k(i2,3).GT.mint(84)+2) goto 170
                ELSE
                  IF(k(k(i2,3),3).GT.mint(83)+6) goto 170
                ENDIF
              ELSE
                i2=mint(53)+ip2-npart
              ENDIF
              IF(kchg(pycomp(k(i2,2)),1).EQ.0) goto 170
              pm2inv=(p(i,4)+p(i2,4))**2-(p(i,1)+p(i2,1))**2-
     &        (p(i,2)+p(i2,2))**2-(p(i,3)+p(i2,3))**2
              IF(pm2inv.LT.pm2min) THEN
                irnew=i2
                pm2min=pm2inv
              ENDIF
  170       CONTINUE
            IF(irnew.EQ.0) THEN
              nevol=nevol-1
              goto 210
            ENDIF
 
C...Begin search for charge recoiler when MODE = 2.
          ELSE
            irold=i
C...Pick sister by history; step up if parton already branched.
  180       IF(k(irold,3).GT.0.AND.k(k(irold,3),2).EQ.k(irold,2)) THEN
              irold=k(irold,3)
              goto 180
            ENDIF
            IF(irold.GT.1.AND.k(irold-1,3).EQ.k(irold,3)) THEN
              irnew=irold-1
            ELSEIF(irold.LT.n.AND.k(irold+1,3).EQ.k(irold,3)) THEN
              irnew=irold+1
C...Last resort: pick at random among other primaries.
            ELSE
              istep=max(1,min(npart-1,int(1d0+(npart-1)*pyr(0))))
              irnew=ipart(1+mod(ip+istep-1,npart))
            ENDIF
C...Trace down if sister branched.
  190       IF(k(irnew,1).GT.10) THEN
              DO 200 ir=irnew+1,n
                IF(k(ir,3).EQ.irnew.AND.k(ir,2).EQ.k(irnew,2)) THEN
                  irnew=ir
                  goto 190
                ENDIF
  200         CONTINUE
            ENDIF
          ENDIF
          irec(nevol)=irnew
        ENDIF
 
C...End loop to set up showering partons. System invariant mass.
  210 CONTINUE
      IF(nevol.LE.0) RETURN
      psum(5)=sqrt(max(0d0,psum(4)**2-psum(1)**2-psum(2)**2-psum(3)**2))
 
C...Check if 3-jet matrix elements to be used.
      m3jc=0
      alpha=0.5d0
      nmesys=0
      IF(mstj(47).GE.1) THEN
 
C...Identify source: q(1), ~q(2), V(3), S(4), chi(5), ~g(6), unknown(0).
        kfsrce=0
        ipart1=k(ipart(1),3)
        ipart2=k(ipart(2),3)
  220   IF(ipart1.EQ.ipart2.AND.ipart1.GT.0) THEN
          kfsrce=iabs(k(ipart1,2))
        ELSEIF(ipart1.GT.ipart2.AND.ipart2.GT.0) THEN
          ipart1=k(ipart1,3)
          goto 220
        ELSEIF(ipart2.GT.ipart1.AND.ipart1.GT.0) THEN
          ipart2=k(ipart2,3)
          goto 220
        ENDIF
        itypes=0
        IF(kfsrce.GE.1.AND.kfsrce.LE.8) itypes=1
        IF(kfsrce.GE.ksusy1+1.AND.kfsrce.LE.ksusy1+8) itypes=2
        IF(kfsrce.GE.ksusy2+1.AND.kfsrce.LE.ksusy2+8) itypes=2
        IF(kfsrce.GE.21.AND.kfsrce.LE.24) itypes=3
        IF(kfsrce.GE.32.AND.kfsrce.LE.34) itypes=3
        IF(kfsrce.EQ.25.OR.(kfsrce.GE.35.AND.kfsrce.LE.37)) itypes=4
        IF(kfsrce.GE.ksusy1+22.AND.kfsrce.LE.ksusy1+37) itypes=5
        IF(kfsrce.EQ.ksusy1+21) itypes=6
 
C...Identify two primary showerers.
        kfla1=iabs(k(ipart(1),2))
        itype1=0
        IF(kfla1.GE.1.AND.kfla1.LE.8) itype1=1
        IF(kfla1.GE.ksusy1+1.AND.kfla1.LE.ksusy1+8) itype1=2
        IF(kfla1.GE.ksusy2+1.AND.kfla1.LE.ksusy2+8) itype1=2
        IF(kfla1.GE.21.AND.kfla1.LE.24) itype1=3
        IF(kfla1.GE.32.AND.kfla1.LE.34) itype1=3
        IF(kfla1.EQ.25.OR.(kfla1.GE.35.AND.kfla1.LE.37)) itype1=4
        IF(kfla1.GE.ksusy1+22.AND.kfla1.LE.ksusy1+37) itype1=5
        IF(kfla1.EQ.ksusy1+21) itype1=6
        kfla2=iabs(k(ipart(2),2))
        itype2=0
        IF(kfla2.GE.1.AND.kfla2.LE.8) itype2=1
        IF(kfla2.GE.ksusy1+1.AND.kfla2.LE.ksusy1+8) itype2=2
        IF(kfla2.GE.ksusy2+1.AND.kfla2.LE.ksusy2+8) itype2=2
        IF(kfla2.GE.21.AND.kfla2.LE.24) itype2=3
        IF(kfla2.GE.32.AND.kfla2.LE.34) itype2=3
        IF(kfla2.EQ.25.OR.(kfla2.GE.35.AND.kfla2.LE.37)) itype2=4
        IF(kfla2.GE.ksusy1+22.AND.kfla2.LE.ksusy1+37) itype2=5
        IF(kfla2.EQ.ksusy1+21) itype2=6
 
C...Order of showerers. Presence of gluino.
        itypmn=min(itype1,itype2)
        itypmx=max(itype1,itype2)
        iord=1
        IF(itype1.GT.itype2) iord=2
        iglui=0
        IF(itype1.EQ.6.OR.itype2.EQ.6) iglui=1
 
C...Require exactly two primary showerers for ME corrections.
        nprim=0
        DO 230 i=1,n
          IF(k(i,3).EQ.ipart1.AND.k(i,2).NE.k(ipart1,2)) nprim=nprim+1
  230   CONTINUE
        IF(nprim.NE.2) THEN
 
C...Predetermined and default matrix element kinds.
        ELSEIF(mstj(38).NE.0) THEN
          m3jc=mstj(38)
          alpha=parj(80)
          mstj(38)=0
        ELSEIF(mstj(47).GE.6) THEN
          m3jc=mstj(47)
        ELSE
          iclass=1
          icombi=4
 
C...Vector/axial vector -> q + qbar; q -> q + V.
          IF(itypmn.EQ.1.AND.itypmx.EQ.1.AND.(itypes.EQ.0.OR.
     &    itypes.EQ.3)) THEN
            iclass=2
            IF(kfsrce.EQ.21.OR.kfsrce.EQ.22) THEN
              icombi=1
            ELSEIF(kfsrce.EQ.23.OR.(kfsrce.EQ.0.AND.
     &      k(ipart(1),2)+k(ipart(2),2).EQ.0)) THEN
C...gamma*/Z0: assume e+e- initial state if unknown.
              ei=-1d0
              IF(kfsrce.EQ.23) THEN
                iannfl=ipart1
                IF(k(iannfl,2).EQ.23) iannfl=k(iannfl,3)
                IF(iannfl.GT.0) THEN
                  IF(k(iannfl,2).EQ.23) iannfl=k(iannfl,3)
                ENDIF
                IF(iannfl.NE.0) THEN
                  kannfl=iabs(k(iannfl,2))
                  IF(kannfl.GE.1.AND.kannfl.LE.18) ei=kchg(kannfl,1)/3d0
                ENDIF
              ENDIF
              ai=sign(1d0,ei+0.1d0)
              vi=ai-4d0*ei*paru(102)
              ef=kchg(kfla1,1)/3d0
              af=sign(1d0,ef+0.1d0)
              vf=af-4d0*ef*paru(102)
              xwc=1d0/(16d0*paru(102)*(1d0-paru(102)))
              sh=psum(5)**2
              sqmz=pmas(23,1)**2
              sqwz=psum(5)*pmas(23,2)
              sbwz=1d0/((sh-sqmz)**2+sqwz**2)
              vect=ei**2*ef**2+2d0*ei*vi*ef*vf*xwc*sh*(sh-sqmz)*sbwz+
     &        (vi**2+ai**2)*vf**2*xwc**2*sh**2*sbwz
              axiv=(vi**2+ai**2)*af**2*xwc**2*sh**2*sbwz
              icombi=3
              alpha=vect/(vect+axiv)
            ELSEIF(kfsrce.EQ.24.OR.kfsrce.EQ.0) THEN
              icombi=4
            ENDIF
C...For chi -> chi q qbar, use V/A -> q qbar as first approximation.
          ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.1.AND.itypes.EQ.5) THEN
            iclass=2
          ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.3.AND.(itypes.EQ.0.OR.
     &    itypes.EQ.1)) THEN
            iclass=3
 
C...Scalar/pseudoscalar -> q + qbar; q -> q + S.
          ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.1.AND.itypes.EQ.4) THEN
            iclass=4
            IF(kfsrce.EQ.25.OR.kfsrce.EQ.35.OR.kfsrce.EQ.37) THEN
              icombi=1
            ELSEIF(kfsrce.EQ.36) THEN
              icombi=2
            ENDIF
          ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.4.AND.(itypes.EQ.0.OR.
     &    itypes.EQ.1)) THEN
            iclass=5
 
C...V -> ~q + ~qbar; ~q -> ~q + V; S -> ~q + ~qbar; ~q -> ~q + S.
          ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.2.AND.(itypes.EQ.0.OR.
     &    itypes.EQ.3)) THEN
            iclass=6
          ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.3.AND.(itypes.EQ.0.OR.
     &    itypes.EQ.2)) THEN
            iclass=7
          ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.2.AND.itypes.EQ.4) THEN
            iclass=8
          ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.4.AND.(itypes.EQ.0.OR.
     &    itypes.EQ.2)) THEN
            iclass=9
 
C...chi -> q + ~qbar; ~q -> q + chi; q -> ~q + chi.
          ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.2.AND.(itypes.EQ.0.OR.
     &    itypes.EQ.5)) THEN
            iclass=10
          ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.5.AND.(itypes.EQ.0.OR.
     &    itypes.EQ.2)) THEN
            iclass=11
          ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.5.AND.(itypes.EQ.0.OR.
     &    itypes.EQ.1)) THEN
            iclass=12
 
C...~g -> q + ~qbar; ~q -> q + ~g; q -> ~q + ~g.
          ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.2.AND.itypes.EQ.6) THEN
            iclass=13
          ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.6.AND.(itypes.EQ.0.OR.
     &    itypes.EQ.2)) THEN
            iclass=14
          ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.6.AND.(itypes.EQ.0.OR.
     &    itypes.EQ.1)) THEN
            iclass=15
 
C...g -> ~g + ~g (eikonal approximation).
          ELSEIF(itypmn.EQ.6.AND.itypmx.EQ.6.AND.itypes.EQ.0) THEN
            iclass=16
          ENDIF
          m3jc=5*iclass+icombi
        ENDIF
 
C...Store pair that together define matrix element treatment.
        IF(m3jc.NE.0) THEN
          nmesys=1
          mesys(nmesys,0)=m3jc
          mesys(nmesys,1)=ipart(1)
          mesys(nmesys,2)=ipart(2)
        ENDIF
 
C...Store qqbar or l+l- pairs for QED radiation.
        IF(kfla1.LE.18.AND.kfla2.LE.18) THEN
          nmesys=nmesys+1
          mesys(nmesys,0)=101
          IF(k(ipart(1),2)+k(ipart(2),2).EQ.0) mesys(nmesys,0)=102
          mesys(nmesys,1)=ipart(1)
          mesys(nmesys,2)=ipart(2)
        ENDIF
 
C...Store other qqbar/l+l- pairs from g/gamma branchings.
        DO 270 i1=1,n
          IF(k(i1,1).GT.10.OR.iabs(k(i1,2)).GT.18) goto 270
          i1m=k(i1,3)
  240     IF(i1m.GT.0.AND.k(i1m,2).EQ.k(i1,2)) THEN
            i1m=k(i1m,3)
            goto 240
          ENDIF
C...Move up this check to avoid out-of-bounds.
          IF(i1m.EQ.0) goto 270
          IF(k(i1m,2).NE.21.AND.k(i1m,2).NE.22) goto 270
          DO 260 i2=i1+1,n
            IF(k(i2,1).GT.10.OR.k(i2,2)+k(i1,2).NE.0) goto 260
            i2m=k(i2,3)
  250       IF(i2m.GT.0.AND.k(i2m,2).EQ.k(i2,2)) THEN
              i2m=k(i2m,3)
              goto 250
            ENDIF
            IF(i1m.EQ.i2m.AND.i1m.GT.0) THEN
              nmesys=nmesys+1
              mesys(nmesys,0)=66
              mesys(nmesys,1)=i1
              mesys(nmesys,2)=i2
              nmesys=nmesys+1
              mesys(nmesys,0)=102
              mesys(nmesys,1)=i1
              mesys(nmesys,2)=i2
            ENDIF
  260     CONTINUE
  270   CONTINUE
      ENDIF
 
C..Loopback point for counting number of emissions.
      ngen=0
  280 ngen=ngen+1
 
C...Begin loop to evolve all existing partons, if required.
  290 imx=0
      pt2mx=0d0
      DO 360 ievol=1,nevol
        IF(iflg(ievol).EQ.0) THEN
 
C...Basic info on radiator and recoil.
          i=ipos(ievol)
          ir=irec(ievol)
          sht=shat(i,ir)
          pm2i=p(i,5)**2
          pm2r=p(ir,5)**2
 
C...Invariant mass of "dipole".Starting value for pT evolution.
          shtcor=(sqrt(sht)-p(ir,5))**2-pm2i
          pt2=min(pt2cmx,0.25d0*shtcor,ptsca(ievol)**2)
 
C...Case of evolution by QCD branching.
          IF(iscol(ievol).NE.0) THEN
 
C...Parton-by-parton maximum scale from initial conditions.
          IF(mstp(72).EQ.0) THEN
            DO 300 iprt=1,nparts
              IF(ir.EQ.ipart(iprt)) pt2=min(pt2,ptpart(iprt)**2)
  300       CONTINUE
          ENDIF
 
C...If kinematically impossible then do not evolve.
            IF(pt2.LT.pt2cmn) THEN
              iflg(ievol)=-1
              goto 360
            ENDIF
 
C...Check if part of system for which ME corrections should be applied.
            imesys=0
            DO 310 ime=1,nmesys
              IF((i.EQ.mesys(ime,1).OR.i.EQ.mesys(ime,2)).AND.
     &        mesys(ime,0).LT.100) imesys=ime
  310       CONTINUE
 
C...Special flag for colour octet states.
            moct=0
            IF(k(i,2).EQ.21) moct=1
            IF(k(i,2).EQ.ksusy1+21) moct=2
 
C...Upper estimate for matrix element weighting and colour factor.
C...Note that g->gg and g->qqbar is split on two sides = "dipoles".
            wtpsgl=2d0
            colfac=4d0/3d0
            IF(moct.GE.1) colfac=3d0/2d0
            IF(iglui.EQ.1.AND.imesys.EQ.1.AND.moct.EQ.0) colfac=3d0
            wtpsqq=0.5d0*0.5d0*nflav
 
C...Determine overestimated z range: switch at c and b masses.
  320       izrg=1
            pt2mne=pt2cmn
            b0=27d0/6d0
            alams=alam3s
            IF(pt2.GT.1.01d0*pmcs) THEN
              izrg=2
              pt2mne=pmcs
              b0=25d0/6d0
              alams=alam4s
            ENDIF
            IF(pt2.GT.1.01d0*pmbs) THEN
              izrg=3
              pt2mne=pmbs
              b0=23d0/6d0
              alams=alam5s
            ENDIF
            zmncut=0.5d0-sqrt(max(0d0,0.25d0-pt2mne/shtcor))
            IF(zmncut.LT.1d-8) zmncut=pt2mne/shtcor
 
C...Find evolution coefficients for q->qg/g->gg and g->qqbar.
            evemgl=wtpsgl*colfac*log(1d0/zmncut-1d0)/b0
            evcoef=evemgl
            IF(moct.EQ.1) THEN
              evemqq=wtpsqq*(1d0-2d0*zmncut)/b0
              evcoef=evcoef+evemqq
            ENDIF
 
C...Pick pT2 (in overestimated z range).
  330       pt2=alams*(pt2/alams)**(pyr(0)**(1d0/evcoef))
 
C...Loopback if crossed c/b mass thresholds.
            IF(izrg.EQ.3.AND.pt2.LT.pmbs) THEN
              pt2=pmbs
              goto 320
            ENDIF
            IF(izrg.EQ.2.AND.pt2.LT.pmcs) THEN
              pt2=pmcs
              goto 320
            ENDIF
 
C...Finish if below lower cutoff.
            IF(pt2.LT.pt2cmn) THEN
              iflg(ievol)=-1
              goto 360
            ENDIF
 
C...Pick kind of branching: q->qg/g->gg/X->Xg or g->qqbar.
            iflag=1
            IF(moct.EQ.1.AND.evemgl.LT.pyr(0)*evcoef) iflag=2
 
C...Pick z: dz/(1-z) or dz.
            IF(iflag.EQ.1) THEN
              z=1d0-zmncut*(1d0/zmncut-1d0)**pyr(0)
            ELSE
              z=zmncut+pyr(0)*(1d0-2d0*zmncut)
            ENDIF
 
C...Loopback if outside allowed range for given pT2.
            zmnnow=0.5d0-sqrt(max(0d0,0.25d0-pt2/shtcor))
            IF(zmnnow.LT.1d-8) zmnnow=pt2/shtcor
            IF(z.LE.zmnnow.OR.z.GE.1d0-zmnnow) goto 330
            pm2=pm2i+pt2/(z*(1d0-z))
            IF(z*(1d0-z).LE.pm2*sht/(sht+pm2-pm2r)**2) goto 330
 
C...No weighting for primary partons; to be done later on.
            IF(imesys.GT.0) THEN
 
C...Weighting of q->qg/X->Xg branching.
            ELSEIF(iflag.EQ.1.AND.moct.NE.1) THEN
              IF(1d0+z**2.LT.wtpsgl*pyr(0)) goto 330
 
C...Weighting of g->gg branching.
            ELSEIF(iflag.EQ.1) THEN
              IF(1d0+z**3.LT.wtpsgl*pyr(0)) goto 330
 
C...Flavour choice and weighting of g->qqbar branching.
            ELSE
              kfq=min(5,1+int(nflav*pyr(0)))
              pmq=pmas(kfq,1)
              rootqq=sqrt(max(0d0,1d0-4d0*pmq**2/pm2))
              wtme=rootqq*(z**2+(1d0-z)**2)
              IF(wtme.LT.pyr(0)) goto 330
              iflag=10+kfq
            ENDIF
 
C...Case of evolution by QED branching.
          ELSEIF(ischg(ievol).NE.0) THEN
 
C...If kinematically impossible then do not evolve.
            pt2emn=pt0eq**2
            IF(iabs(k(i,2)).GT.10) pt2emn=pt0el**2
            IF(pt2.LT.pt2emn) THEN
              iflg(ievol)=-1
              goto 360
            ENDIF
 
C...Check if part of system for which ME corrections should be applied.
           imesys=0
            DO 340 ime=1,nmesys
              IF((i.EQ.mesys(ime,1).OR.i.EQ.mesys(ime,2)).AND.
     &        mesys(ime,0).GT.100) imesys=ime
  340      CONTINUE
 
C...Charge. Matrix element weighting factor.
            chg=ischg(ievol)/3d0
            wtpsga=2d0
 
C...Determine overestimated z range. Find evolution coefficient.
            zmncut=0.5d0-sqrt(max(0d0,0.25d0-pt2emn/shtcor))
            IF(zmncut.LT.1d-8) zmncut=pt2emn/shtcor
            evcoef=aem2pi*chg**2*wtpsga*log(1d0/zmncut-1d0)
 
C...Pick pT2 (in overestimated z range).
  350       pt2=pt2*pyr(0)**(1d0/evcoef)
 
C...Finish if below lower cutoff.
            IF(pt2.LT.pt2emn) THEN
              iflg(ievol)=-1
              goto 360
            ENDIF
 
C...Pick z: dz/(1-z).
            z=1d0-zmncut*(1d0/zmncut-1d0)**pyr(0)
 
C...Loopback if outside allowed range for given pT2.
            zmnnow=0.5d0-sqrt(max(0d0,0.25d0-pt2/shtcor))
            IF(zmnnow.LT.1d-8) zmnnow=pt2/shtcor
            IF(z.LE.zmnnow.OR.z.GE.1d0-zmnnow) goto 350
            pm2=pm2i+pt2/(z*(1d0-z))
            IF(z*(1d0-z).LE.pm2*sht/(sht+pm2-pm2r)**2) goto 350
 
C...Weighting by branching kernel, except if ME weighting later.
            IF(imesys.EQ.0) THEN
              IF(1d0+z**2.LT.wtpsga*pyr(0)) goto 350
            ENDIF
            iflag=3
          ENDIF
 
C...Save acceptable branching.
          iflg(ievol)=iflag
          imesav(ievol)=imesys
          pt2sav(ievol)=pt2
          zsav(ievol)=z
          shtsav(ievol)=sht
        ENDIF
 
C...Check if branching has highest pT.
        IF(iflg(ievol).GE.1.AND.pt2sav(ievol).GT.pt2mx) THEN
          imx=ievol
          pt2mx=pt2sav(ievol)
        ENDIF
  360 CONTINUE
 
C...Finished if no more branchings to be done.
      IF(imx.EQ.0) goto 480
 
C...Restore info on hardest branching to be processed.
      i=ipos(imx)
      ir=irec(imx)
      kcol=iscol(imx)
      kcha=ischg(imx)
      imesys=imesav(imx)
      pt2=pt2sav(imx)
      z=zsav(imx)
      sht=shtsav(imx)
      pm2i=p(i,5)**2
      pm2r=p(ir,5)**2
      pm2=pm2i+pt2/(z*(1d0-z))
 
C...Special flag for colour octet states.
      moct=0
      IF(k(i,2).EQ.21) moct=1
      IF(k(i,2).EQ.ksusy1+21) moct=2
 
C...Restore further info for g->qqbar branching.
      kfq=0
      IF(iflg(imx).GT.10) THEN
        kfq=iflg(imx)-10
        pmq=pmas(kfq,1)
        rootqq=sqrt(max(0d0,1d0-4d0*pmq**2/pm2))
      ENDIF
 
C...For branching g include azimuthal asymmetries from polarization.
      asypol=0d0
      IF(moct.EQ.1.AND.mod(mstj(46),2).EQ.1) THEN
C...Trace grandmother via intermediate recoil copies.
        kfgm=0
        im=i
  370   IF(k(im,3).NE.k(im-1,3).AND.k(im,3).NE.k(im+1,3).AND.
     &  k(im,3).GT.0) THEN
          im=k(im,3)
          IF(im.GT.mint(84)) goto 370
        ENDIF
        igm=k(im,3)
        IF(igm.GT.mint(84).AND.igm.LT.im.AND.im.LE.i)
     &  kfgm=iabs(k(igm,2))
C...Define approximate energy sharing by identifying aunt.
        iau=im+1
        IF(iau.GT.n-3.OR.k(iau,3).NE.igm) iau=im-1
        IF(kfgm.NE.0.AND.(kfgm.LE.6.OR.kfgm.EQ.21)) THEN
          zold=p(im,4)/(p(im,4)+p(iau,4))
C...Coefficient from gluon production.
          IF(kfgm.LE.6) THEN
            asypol=2d0*(1d0-zold)/(1d0+(1d0-zold)**2)
          ELSE
            asypol=((1d0-zold)/(1d0-zold*(1d0-zold)))**2
          ENDIF
C...Coefficient from gluon decay.
          IF(kfq.EQ.0) THEN
            asypol=asypol*(z*(1d0-z)/(1d0-z*(1d0-z)))**2
          ELSE
            asypol=-asypol*2d0*z*(1d0-z)/(1d0-2d0*z*(1d0-z))
          ENDIF
        ENDIF
      ENDIF
 
C...Create new slots for branching products and recoil.
      inew=n+1
      ignew=n+2
      irnew=n+3
      n=n+3
 
C...Set status, flavour and mother of new ones.
      k(inew,1)=k(i,1)
      k(ignew,1)=3
      IF(kcha.NE.0)  k(ignew,1)=1
      k(irnew,1)=k(ir,1)
      IF(kfq.EQ.0) THEN
        k(inew,2)=k(i,2)
        k(ignew,2)=21
        IF(kcha.NE.0)  k(ignew,2)=22
      ELSE
        k(inew,2)=-isign(kfq,kcol)
        k(ignew,2)=-k(inew,2)
      ENDIF
      k(irnew,2)=k(ir,2)
      k(inew,3)=i
      k(ignew,3)=i
      k(irnew,3)=ir
 
C...Find rest frame and angles of branching+recoil.
      DO 380 j=1,5
        p(inew,j)=p(i,j)
        p(ignew,j)=0d0
        p(irnew,j)=p(ir,j)
  380 CONTINUE
      betax=(p(inew,1)+p(irnew,1))/(p(inew,4)+p(irnew,4))
      betay=(p(inew,2)+p(irnew,2))/(p(inew,4)+p(irnew,4))
      betaz=(p(inew,3)+p(irnew,3))/(p(inew,4)+p(irnew,4))
      CALL pyrobo(inew,irnew,0d0,0d0,-betax,-betay,-betaz)
      phi=pyangl(p(inew,1),p(inew,2))
      theta=pyangl(p(inew,3),sqrt(p(inew,1)**2+p(inew,2)**2))
 
C...Derive kinematics of branching: generics (like g->gg).
      DO 390 j=1,4
        p(inew,j)=0d0
        p(irnew,j)=0d0
  390 CONTINUE
      pem=0.5d0*(sht+pm2-pm2r)/sqrt(sht)
      pzm=0.5d0*sqrt(max(0d0,(sht-pm2-pm2r)**2-4d0*pm2*pm2r)/sht)
      pt2cor=pm2*(pem**2*z*(1d0-z)-0.25d0*pm2)/pzm**2
      ptcor=sqrt(max(0d0,pt2cor))
      pzn=(pem**2*z-0.5d0*pm2)/pzm
      pzg=(pem**2*(1d0-z)-0.5d0*pm2)/pzm
C...Specific kinematics reduction for q->qg with m_q > 0.
      IF(moct.NE.1) THEN
        ptcor=(1d0-pm2i/pm2)*ptcor
        pzn=pzn+pm2i*pzg/pm2
        pzg=(1d0-pm2i/pm2)*pzg
C...Specific kinematics reduction for g->qqbar with m_q > 0.
      ELSEIF(kfq.NE.0) THEN
        p(inew,5)=pmq
        p(ignew,5)=pmq
        ptcor=rootqq*ptcor
        pzn=0.5d0*((1d0+rootqq)*pzn+(1d0-rootqq)*pzg)
        pzg=pzm-pzn
      ENDIF
 
C...Pick phi and construct kinematics of branching.
  400 phirot=paru(2)*pyr(0)
      p(inew,1)=ptcor*cos(phirot)
      p(inew,2)=ptcor*sin(phirot)
      p(inew,3)=pzn
      p(inew,4)=sqrt(ptcor**2+p(inew,3)**2+p(inew,5)**2)
      p(ignew,1)=-p(inew,1)
      p(ignew,2)=-p(inew,2)
      p(ignew,3)=pzg
      p(ignew,4)=sqrt(ptcor**2+p(ignew,3)**2+p(ignew,5)**2)
      p(irnew,1)=0d0
      p(irnew,2)=0d0
      p(irnew,3)=-pzm
      p(irnew,4)=0.5d0*(sht+pm2r-pm2)/sqrt(sht)
 
C...Boost branching system to lab frame.
      CALL pyrobo(inew,irnew,theta,phi,betax,betay,betaz)
 
C...Renew choice of phi angle according to polarization asymmetry.
      IF(abs(asypol).GT.1d-3) THEN
        DO 410 j=1,3
          dpt(1,j)=p(i,j)
          dpt(2,j)=p(iau,j)
          dpt(3,j)=p(inew,j)
  410   CONTINUE
        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 420 j=1,3
          dpt(4,j)=dpt(2,j)-dpma*dpt(1,j)/max(1d-10,dpmm)
          dpt(5,j)=dpt(3,j)-dpmd*dpt(1,j)/max(1d-10,dpmm)
  420   CONTINUE
        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.1d0*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(1d0+asypol*(2d0*cad**2-1d0).LT.pyr(0)*(1d0+abs(asypol)))
     &    goto 400
        ENDIF
      ENDIF
 
C...Matrix element corrections for primary partons when requested.
      IF(imesys.GT.0) THEN
        m3jc=mesys(imesys,0)
 
C...Identify recoiling partner and set up three-body kinematics.
        irp=mesys(imesys,1)
        IF(irp.EQ.i) irp=mesys(imesys,2)
        IF(irp.EQ.ir) irp=irnew
        DO 430 j=1,4
          psum(j)=p(inew,j)+p(irp,j)+p(ignew,j)
  430   CONTINUE
        psum(5)=sqrt(max(0d0,psum(4)**2-psum(1)**2-psum(2)**2-
     &  psum(3)**2))
        x1=2d0*(psum(4)*p(inew,4)-psum(1)*p(inew,1)-psum(2)*p(inew,2)-
     &  psum(3)*p(inew,3))/psum(5)**2
        x2=2d0*(psum(4)*p(irp,4)-psum(1)*p(irp,1)-psum(2)*p(irp,2)-
     &  psum(3)*p(irp,3))/psum(5)**2
        x3=2d0-x1-x2
        r1me=p(inew,5)/psum(5)
        r2me=p(irp,5)/psum(5)
 
C...Matrix elements for gluon emission.
        IF(m3jc.LT.100) THEN
 
C...Call ME, with right order important for two inequivalent showerers.
          IF(mesys(imesys,iord).EQ.i) THEN
            wme=pymael(m3jc,x1,x2,r1me,r2me,alpha)
          ELSE
            wme=pymael(m3jc,x2,x1,r2me,r1me,alpha)
          ENDIF
 
C...Split up total ME when two radiating partons.
          isprad=1
          IF((m3jc.GE.16.AND.m3jc.LE.19).OR.(m3jc.GE.26.AND.m3jc.LE.29)
     &    .OR.(m3jc.GE.36.AND.m3jc.LE.39).OR.(m3jc.GE.46.AND.m3jc.LE.49)
     &    .OR.(m3jc.GE.56.AND.m3jc.LE.64)) isprad=0
          IF(isprad.EQ.1) wme=wme*max(1d-10,1d0+r1me**2-r2me**2-x1)/
     &    max(1d-10,2d0-x1-x2)
 
C...Evaluate shower rate.
          wps=2d0/(max(1d-10,2d0-x1-x2)*
     &    max(1d-10,1d0+r2me**2-r1me**2-x2))
          IF(iglui.EQ.1) wps=(9d0/4d0)*wps
 
C...Matrix elements for photon emission: still rather primitive.
        ELSE
 
C...For generic charge combination currently only massless expression.
          IF(m3jc.EQ.101) THEN
            chg1=kchg(pycomp(k(i,2)),1)*isign(1,k(i,2))/3d0
            chg2=kchg(pycomp(k(irp,2)),1)*isign(1,k(irp,2))/3d0
            wme=(chg1*(1d0-x1)/x3-chg2*(1d0-x2)/x3)**2*(x1**2+x2**2)
            wps=2d0*(chg1**2*(1d0-x1)/x3+chg2**2*(1d0-x2)/x3)
 
C...For flavour neutral system assume vector source and include masses.
          ELSE
            wme=pymael(11,x1,x2,r1me,r2me,0d0)*max(1d-10,
     &      1d0+r1me**2-r2me**2-x1)/max(1d-10,2d0-x1-x2)
            wps=2d0/(max(1d-10,2d0-x1-x2)*
     &      max(1d-10,1d0+r2me**2-r1me**2-x2))
          ENDIF
        ENDIF
 
C...Perform weighting with W_ME/W_PS.
        IF(wme.LT.pyr(0)*wps) THEN
          n=n-3
          iflg(imx)=0
          goto 290
        ENDIF
      ENDIF
 
C...Now for sure accepted branching. Save highest pT.
      IF(ngen.EQ.1) ptgen=sqrt(pt2)
 
C...Update status for obsolete ones. Bookkkep the moved original parton
C...and new daughter (arbitrary choice for g->gg or g->qqbar).
C...Do not bookkeep radiated photon, since it cannot radiate further.
      k(i,1)=k(i,1)+10
      k(ir,1)=k(ir,1)+10
      DO 440 ip=1,npart
        IF(ipart(ip).EQ.i) ipart(ip)=inew
        IF(ipart(ip).EQ.ir) ipart(ip)=irnew
  440 CONTINUE
      IF(kcha.EQ.0) THEN
        npart=npart+1
        ipart(npart)=ignew
      ENDIF
 
C...Initialize colour flow of branching.
C...Use both old and new style colour tags for flexibility.
      k(inew,4)=0
      k(ignew,4)=0
      k(inew,5)=0
      k(ignew,5)=0
      jcolp=4+(1-kcol)/2
      jcoln=9-jcolp
      mct(inew,1)=0
      mct(inew,2)=0
      mct(ignew,1)=0
      mct(ignew,2)=0
      mct(irnew,1)=0
      mct(irnew,2)=0
 
C...Trivial colour flow for l->lgamma and q->qgamma.
      IF(iabs(kcha).EQ.3) THEN
        k(i,4)=inew
        k(i,5)=ignew
      ELSEIF(kcha.NE.0) THEN
        IF(k(i,4).NE.0) THEN
          k(i,4)=k(i,4)+inew
          k(inew,4)=mstu(5)*i
          mct(inew,1)=mct(i,1)
        ENDIF
        IF(k(i,5).NE.0) THEN
          k(i,5)=k(i,5)+inew
          k(inew,5)=mstu(5)*i
          mct(inew,2)=mct(i,2)
        ENDIF
 
C...Set colour flow for q->qg and g->gg.
      ELSEIF(kfq.EQ.0) THEN
        k(i,jcolp)=k(i,jcolp)+ignew
        k(ignew,jcolp)=mstu(5)*i
        k(inew,jcolp)=mstu(5)*ignew
        k(ignew,jcoln)=mstu(5)*inew
        mct(ignew,jcolp-3)=mct(i,jcolp-3)
        nct=nct+1
        mct(inew,jcolp-3)=nct
        mct(ignew,jcoln-3)=nct
        IF(moct.GE.1) THEN
          k(i,jcoln)=k(i,jcoln)+inew
          k(inew,jcoln)=mstu(5)*i
          mct(inew,jcoln-3)=mct(i,jcoln-3)
        ENDIF
 
C...Set colour flow for g->qqbar.
      ELSE
        k(i,jcoln)=k(i,jcoln)+inew
        k(inew,jcoln)=mstu(5)*i
        k(i,jcolp)=k(i,jcolp)+ignew
        k(ignew,jcolp)=mstu(5)*i
        mct(inew,jcoln-3)=mct(i,jcoln-3)
        mct(ignew,jcolp-3)=mct(i,jcolp-3)
      ENDIF
 
C...Daughter info for colourless recoiling parton.
      IF(k(ir,4).EQ.0.AND.k(ir,5).EQ.0) THEN
        k(ir,4)=irnew
        k(ir,5)=irnew
        k(irnew,4)=0
        k(irnew,5)=0
 
C...Colour of recoiling parton sails through unchanged.
      ELSE
        IF(k(ir,4).NE.0) THEN
          k(ir,4)=k(ir,4)+irnew
          k(irnew,4)=mstu(5)*ir
          mct(irnew,1)=mct(ir,1)
        ENDIF
        IF(k(ir,5).NE.0) THEN
          k(ir,5)=k(ir,5)+irnew
          k(irnew,5)=mstu(5)*ir
          mct(irnew,2)=mct(ir,2)
        ENDIF
      ENDIF
 
C...Vertex information trivial.
      DO 450 j=1,5
        v(inew,j)=v(i,j)
        v(ignew,j)=v(i,j)
        v(irnew,j)=v(ir,j)
  450 CONTINUE
 
C...Update list of old radiators.
        DO 460 ievol=1,nevol
          IF(ipos(ievol).EQ.i.AND.irec(ievol).EQ.ir) THEN
            ipos(ievol)=inew
            IF(kcol.NE.0.AND.iscol(ievol).EQ.kcol) ipos(ievol)=ignew
            irec(ievol)=irnew
            iflg(ievol)=0
          ELSEIF(ipos(ievol).EQ.i) THEN
            ipos(ievol)=inew
            iflg(ievol)=0
          ELSEIF(ipos(ievol).EQ.ir.AND.irec(ievol).EQ.i) THEN
            ipos(ievol)=irnew
            irec(ievol)=inew
            IF(kcol.NE.0.AND.iscol(ievol).NE.kcol) irec(ievol)=ignew
            iflg(ievol)=0
          ELSEIF(ipos(ievol).EQ.ir) THEN
            ipos(ievol)=irnew
            iflg(ievol)=0
          ENDIF
C...Update links of old connected partons.
          IF(irec(ievol).EQ.i) THEN
            irec(ievol)=inew
            iflg(ievol)=0
          ELSEIF(irec(ievol).EQ.ir) THEN
            irec(ievol)=irnew
            iflg(ievol)=0
          ENDIF
  460   CONTINUE
 
C...q->qg or g->gg: create new gluon radiators.
      IF(kcol.NE.0.AND.kfq.EQ.0) THEN
        nevol=nevol+1
        ipos(nevol)=inew
        irec(nevol)=ignew
        iflg(nevol)=0
        iscol(nevol)=kcol
        ischg(nevol)=0
        ptsca(nevol)=sqrt(pt2)
        nevol=nevol+1
        ipos(nevol)=ignew
        irec(nevol)=inew
        iflg(nevol)=0
        iscol(nevol)=-kcol
        ischg(nevol)=0
        ptsca(nevol)=ptsca(nevol-1)
      ENDIF
 
C...Update matrix elements parton list and add new for g/gamma->qqbar.
      DO 470 ime=1,nmesys
        IF(mesys(ime,1).EQ.i) mesys(ime,1)=inew
        IF(mesys(ime,2).EQ.i) mesys(ime,2)=inew
        IF(mesys(ime,1).EQ.ir) mesys(ime,1)=irnew
        IF(mesys(ime,2).EQ.ir) mesys(ime,2)=irnew
  470 CONTINUE
      IF(kfq.NE.0) THEN
        nmesys=nmesys+1
        mesys(nmesys,0)=66
        mesys(nmesys,1)=inew
        mesys(nmesys,2)=ignew
        nmesys=nmesys+1
        mesys(nmesys,0)=102
        mesys(nmesys,1)=inew
        mesys(nmesys,2)=ignew
      ENDIF
 
C...Global statistics.
      mint(353)=mint(353)+1
      vint(353)=vint(353)+ptcor
      IF (mint(353).EQ.1) vint(358)=ptcor
 
C...Loopback for more emissions if enough space.
      pt2cmx=pt2
      IF(npart.LT.maxnur-1.AND.nevol.LT.2*maxnur-2.AND.
     &nmesys.LT.maxnur-2.AND.n.LT.mstu(4)-mstu(32)-5) THEN
        goto 280
      ELSE
        CALL pyerrm(11,'(PYPTFS:) no more memory left for shower')
      ENDIF
 
C...Done.
  480 CONTINUE
 
      RETURN
      END