pythiaeRHIC/pythia/pyptfs.f

0001
0002 C*********************************************************************
0003
0004 C...PYPTFS
0005 C...Generates pT-ordered timelike final-state parton showers.
0006
0007 C...MODE defines how to find radiators and recoilers.
0008 C... = 0 : based on colour flow between undecayed partons.
0009 C... = 1 : for IPART <= NPARTD only consider primary partons,
0010 C...       whether decayed or not; else as above.
0011 C... = 2 : based on common history, whether decayed or not.
0012
0013       SUBROUTINE PYPTFS(MODE,PTMAX,PTMIN,PTGEN)
0014
0015 C...Double precision and integer declarations.
0016       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0017       IMPLICIT INTEGER(I-N)
0018       INTEGER PYK,PYCHGE,PYCOMP
0019 C...Parameter statement to help give large particle numbers.
0020       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0021      &KEXCIT=4000000,KDIMEN=5000000)
0022 C...Parameter statement for maximum size of showers.
0023       PARAMETER (MAXNUR=1000)
0024 C...Commonblocks.
0025       COMMON/PYPART/NPART,NPARTD,IPART(MAXNUR),PTPART(MAXNUR)
0026       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0027       COMMON/PYCTAG/NCT,MCT(4000,2)
0028       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0029       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0030       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0031       COMMON/PYINT1/MINT(400),VINT(400)
0032       SAVE /PYPART/,/PYJETS/,/PYCTAG/,/PYDAT1/,/PYDAT2/,/PYPARS/,
0033      &/PYINT1/
0034 C...Local arrays.
0035       DIMENSION IPOS(2*MAXNUR),IREC(2*MAXNUR),IFLG(2*MAXNUR),
0036      &ISCOL(2*MAXNUR),ISCHG(2*MAXNUR),PTSCA(2*MAXNUR),IMESAV(2*MAXNUR),
0037      &PT2SAV(2*MAXNUR),ZSAV(2*MAXNUR),SHTSAV(2*MAXNUR),
0038      &MESYS(MAXNUR,0:2),PSUM(5),DPT(5,4)
0039 C...Statement functions.
0040       SHAT(I,J)=(P(I,4)+P(J,4))**2-(P(I,1)+P(J,1))**2-
0041      &(P(I,2)+P(J,2))**2-(P(I,3)+P(J,3))**2
0042
0043 C...Initial values. Check that valid system.
0044       PTGEN=0D0
0045       IF(MSTJ(41).NE.1.AND.MSTJ(41).NE.2.AND.MSTJ(41).NE.11.AND.
0046      &MSTJ(41).NE.12) RETURN
0047       IF(NPART.LE.0) THEN
0048         CALL PYERRM(2,'(PYPTFS:) showering system too small')
0049         RETURN
0050       ENDIF
0051       PT2CMX=PTMAX**2
0052
0053 C...Mass thresholds and Lambda for QCD evolution.
0054       PMB=PMAS(5,1)
0055       PMC=PMAS(4,1)
0056       ALAM5=PARJ(81)
0057       ALAM4=ALAM5*(PMB/ALAM5)**(2D0/25D0)
0058       ALAM3=ALAM4*(PMC/ALAM4)**(2D0/27D0)
0059       PMBS=PMB**2
0060       PMCS=PMC**2
0061       ALAM5S=ALAM5**2
0062       ALAM4S=ALAM4**2
0063       ALAM3S=ALAM3**2
0064
0065 C...Cutoff scale for QCD evolution. Starting pT2.
0066       NFLAV=MAX(0,MIN(5,MSTJ(45)))
0067       PT0C=0.5D0*PARJ(82)
0068       PT2CMN=MAX(PTMIN,PT0C,1.1D0*ALAM3)**2
0069
0070 C...Parameters for QED evolution.
0071       AEM2PI=PARU(101)/PARU(2)
0072       PT0EQ=0.5D0*PARJ(83)
0073       PT0EL=0.5D0*PARJ(90)
0074
0075 C...Reset. Remove irrelevant colour tags.
0076       NEVOL=0
0077       DO 100 J=1,4
0078         PSUM(J)=0D0
0079   100 CONTINUE
0080       DO 110 I=MINT(84)+1,N
0081         IF(K(I,2).GT.0.AND.K(I,2).LT.6) THEN
0082           K(I,5)=0
0083           MCT(I,2)=0
0084         ENDIF
0085         IF(K(I,2).LT.0.AND.K(I,2).GT.-6) THEN
0086           K(I,4)=0
0087           MCT(I,1)=0
0088         ENDIF
0089   110 CONTINUE
0090       NPARTS=NPART
0091
0092 C...Begin loop to set up showering partons. Sum four-momenta.
0093       DO 210 IP=1,NPART
0094         I=IPART(IP)
0095         IF(MODE.NE.1.OR.I.GT.NPARTD) THEN
0096           IF(K(I,1).GT.10) GOTO 210
0097         ELSEIF(K(I,3).GT.MINT(84)) THEN
0098           IF(K(I,3).GT.MINT(84)+2) GOTO 210
0099         ELSE
0100           IF(K(K(I,3),3).GT.MINT(83)+6) GOTO 210
0101         ENDIF
0102         DO 120 J=1,4
0103           PSUM(J)=PSUM(J)+P(I,J)
0104   120   CONTINUE
0105
0106 C...Find colour and charge, but skip diquarks.
0107         IF(IABS(K(I,2)).GT.1000.AND.IABS(K(I,2)).LT.10000) GOTO 210
0108         KCOL=ISIGN(KCHG(PYCOMP(K(I,2)),2),K(I,2))
0109         KCHA=ISIGN(KCHG(PYCOMP(K(I,2)),1),K(I,2))
0110
0111 C...Either colour or anticolour charge radiates; for gluon both.
0112         DO 160 JSGCOL=1,-1,-2
0113           IF(KCOL.EQ.JSGCOL.OR.KCOL.EQ.2) THEN
0114             JCOL=4+(1-JSGCOL)/2
0115             JCOLR=9-JCOL
0116
0117 C...Basic info about radiating parton.
0118             NEVOL=NEVOL+1
0119             IPOS(NEVOL)=I
0120             IFLG(NEVOL)=0
0121             ISCOL(NEVOL)=JSGCOL
0122             ISCHG(NEVOL)=0
0123             PTSCA(NEVOL)=PTPART(IP)
0124
0125 C...Begin search for colour recoiler when MODE = 0 or 1.
0126             IF(MODE.LE.1) THEN
0127 C...Find sister with matching anticolour to the radiating parton.
0128               IROLD=I
0129               IRNEW=K(IROLD,JCOL)/MSTU(5)
0130               MOVE=1
0131
0132 C...The following will add MCT colour tracing for unprepped events
0133 C...If not done, trace Les Houches colour tags for this dipole
0134 C              IF (MCT(I,JCOL-3).EQ.0) THEN
0135 C                CALL PYCTTR(I,JCOL,INEW)
0136 C...Clean up mother/daughter 'read' tags set by PYCTTR
0137 C                DO 125 IR=1,N
0138 C                  K(IR,4)=MOD(K(IR,4),MSTU(5)**2)
0139 C                  K(IR,5)=MOD(K(IR,5),MSTU(5)**2)
0140 C 125            CONTINUE
0141 C              ENDIF
0142
0143 C...Skip radiation off loose colour ends.
0144   130         IF(IRNEW.EQ.0) THEN
0145                 NEVOL=NEVOL-1
0146                 GOTO 160
0147
0148 C...Optionally skip radiation on dipole to beam remnant.
0149               ELSEIF(MSTP(72).LE.1.AND.IRNEW.GT.MINT(53)) THEN
0150                 NEVOL=NEVOL-1
0151                 GOTO 160
0152
0153 C...For now always skip radiation on dipole to junction.
0154               ELSEIF(K(IRNEW,2).EQ.88) THEN
0155                 NEVOL=NEVOL-1
0156                 GOTO 160
0157
0158 C...For MODE=1: if reached primary then done.
0159               ELSEIF(MODE.EQ.1.AND.IRNEW.GT.MINT(84)+2.AND.
0160      &        IRNEW.LE.NPARTD) THEN
0161
0162 C...If sister stable and points back then done.
0163               ELSEIF(MOVE.EQ.1.AND.K(IRNEW,JCOLR)/MSTU(5).EQ.IROLD)
0164      &        THEN
0165                 IF(K(IRNEW,1).LT.10) THEN
0166
0167 C...If sister unstable then go to her daughter.
0168                 ELSE
0169                   IROLD=IRNEW
0170                   IRNEW=MOD(K(IRNEW,JCOLR),MSTU(5))
0171                   MOVE=2
0172                   GOTO 130
0173                ENDIF
0174
0175 C...If found mother then look for aunt.
0176               ELSEIF(MOVE.EQ.1.AND.MOD(K(IRNEW,JCOL),MSTU(5)).EQ.
0177      &        IROLD) THEN
0178                 IROLD=IRNEW
0179                 IRNEW=K(IROLD,JCOL)/MSTU(5)
0180                 GOTO 130
0181
0182 C...If daughter stable then done.
0183               ELSEIF(MOVE.EQ.2.AND.K(IRNEW,JCOLR)/MSTU(5).EQ.IROLD)
0184      &        THEN
0185                 IF(K(IRNEW,1).LT.10) THEN
0186
0187 C...If daughter unstable then go to granddaughter.
0188                 ELSE
0189                   IROLD=IRNEW
0190                   IRNEW=MOD(K(IRNEW,JCOLR),MSTU(5))
0191                   MOVE=2
0192                   GOTO 130
0193                 ENDIF
0194
0195 C...If daughter points to another daughter then done or move up.
0196               ELSEIF(MOVE.EQ.2.AND.MOD(K(IRNEW,JCOL),MSTU(5)).EQ.
0197      &        IROLD) THEN
0198                 IF(K(IRNEW,1).LT.10) THEN
0199                 ELSE
0200                   IROLD=IRNEW
0201                   IRNEW=K(IRNEW,JCOL)/MSTU(5)
0202                   MOVE=1
0203                   GOTO 130
0204                 ENDIF
0205               ENDIF
0206
0207 C...Begin search for colour recoiler when MODE = 2.
0208             ELSE
0209               IROLD=I
0210               IRNEW=K(IROLD,JCOL)/MSTU(5)
0211   140         IF(K(IRNEW,JCOLR)/MSTU(5).NE.IROLD) THEN
0212 C...Step up to mother if radiating parton already branched.
0213                 IF(K(IRNEW,2).EQ.K(IROLD,2)) THEN
0214                   IROLD=IRNEW
0215                   IRNEW=K(IROLD,JCOL)/MSTU(5)
0216                   GOTO 140
0217 C...Pick sister by history if no anticolour available.
0218                 ELSE
0219                   IF(IROLD.GT.1.AND.K(IROLD-1,3).EQ.K(IROLD,3)) THEN
0220                     IRNEW=IROLD-1
0221                   ELSEIF(IROLD.LT.N.AND.K(IROLD+1,3).EQ.K(IROLD,3))
0222      &            THEN
0223                     IRNEW=IROLD+1
0224 C...Last resort: pick at random among other primaries.
0225                   ELSE
0226                     ISTEP=MAX(1,MIN(NPART-1,INT(1D0+(NPART-1)*PYR(0))))
0227                     IRNEW=IPART(1+MOD(IP+ISTEP-1,NPART))
0228                   ENDIF
0229                 ENDIF
0230               ENDIF
0231 C...Trace down if sister branched.
0232   150         IF(K(IRNEW,1).GT.10) THEN
0233                 IRNEW=MOD(K(IRNEW,JCOLR),MSTU(5))
0234                 GOTO 150
0235               ENDIF
0236             ENDIF
0237
0238 C...Now found other end of colour dipole.
0239             IREC(NEVOL)=IRNEW
0240           ENDIF
0241   160   CONTINUE
0242
0243 C...Also electrical charge may radiate; so far only quarks and leptons.
0244         IF((MSTJ(41).EQ.2.OR.MSTJ(41).EQ.12).AND.KCHA.NE.0.AND.
0245      &  IABS(K(I,2)).LE.18) THEN
0246
0247 C...Basic info about radiating parton.
0248           NEVOL=NEVOL+1
0249           IPOS(NEVOL)=I
0250           IFLG(NEVOL)=0
0251           ISCOL(NEVOL)=0
0252           ISCHG(NEVOL)=KCHA
0253           PTSCA(NEVOL)=PTPART(IP)
0254
0255 C...Pick nearest (= smallest invariant mass) charged particle
0256 C...as recoiler when MODE = 0 or 1 (but for latter among primaries).
0257           IF(MODE.LE.1) THEN
0258             IRNEW=0
0259             PM2MIN=VINT(2)
0260             DO 170 IP2=1,NPART+N-MINT(53)
0261               IF(IP2.EQ.IP) GOTO 170
0262               IF(IP2.LE.NPART) THEN
0263                 I2=IPART(IP2)
0264                 IF(MODE.NE.1.OR.I2.GT.NPARTD) THEN
0265                   IF(K(I2,1).GT.10) GOTO 170
0266                 ELSEIF(K(I2,3).GT.MINT(84)) THEN
0267                   IF(K(I2,3).GT.MINT(84)+2) GOTO 170
0268                 ELSE
0269                   IF(K(K(I2,3),3).GT.MINT(83)+6) GOTO 170
0270                 ENDIF
0271               ELSE
0272                 I2=MINT(53)+IP2-NPART
0273               ENDIF
0274               IF(KCHG(PYCOMP(K(I2,2)),1).EQ.0) GOTO 170
0275               PM2INV=(P(I,4)+P(I2,4))**2-(P(I,1)+P(I2,1))**2-
0276      &        (P(I,2)+P(I2,2))**2-(P(I,3)+P(I2,3))**2
0277               IF(PM2INV.LT.PM2MIN) THEN
0278                 IRNEW=I2
0279                 PM2MIN=PM2INV
0280               ENDIF
0281   170       CONTINUE
0282             IF(IRNEW.EQ.0) THEN
0283               NEVOL=NEVOL-1
0284               GOTO 210
0285             ENDIF
0286
0287 C...Begin search for charge recoiler when MODE = 2.
0288           ELSE
0289             IROLD=I
0290 C...Pick sister by history; step up if parton already branched.
0291   180       IF(K(IROLD,3).GT.0.AND.K(K(IROLD,3),2).EQ.K(IROLD,2)) THEN
0292               IROLD=K(IROLD,3)
0293               GOTO 180
0294             ENDIF
0295             IF(IROLD.GT.1.AND.K(IROLD-1,3).EQ.K(IROLD,3)) THEN
0296               IRNEW=IROLD-1
0297             ELSEIF(IROLD.LT.N.AND.K(IROLD+1,3).EQ.K(IROLD,3)) THEN
0298               IRNEW=IROLD+1
0299 C...Last resort: pick at random among other primaries.
0300             ELSE
0301               ISTEP=MAX(1,MIN(NPART-1,INT(1D0+(NPART-1)*PYR(0))))
0302               IRNEW=IPART(1+MOD(IP+ISTEP-1,NPART))
0303             ENDIF
0304 C...Trace down if sister branched.
0305   190       IF(K(IRNEW,1).GT.10) THEN
0306               DO 200 IR=IRNEW+1,N
0307                 IF(K(IR,3).EQ.IRNEW.AND.K(IR,2).EQ.K(IRNEW,2)) THEN
0308                   IRNEW=IR
0309                   GOTO 190
0310                 ENDIF
0311   200         CONTINUE
0312             ENDIF
0313           ENDIF
0314           IREC(NEVOL)=IRNEW
0315         ENDIF
0316
0317 C...End loop to set up showering partons. System invariant mass.
0318   210 CONTINUE
0319       IF(NEVOL.LE.0) RETURN
0320       PSUM(5)=SQRT(MAX(0D0,PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2))
0321
0322 C...Check if 3-jet matrix elements to be used.
0323       M3JC=0
0324       ALPHA=0.5D0
0325       NMESYS=0
0326       IF(MSTJ(47).GE.1) THEN
0327
0328 C...Identify source: q(1), ~q(2), V(3), S(4), chi(5), ~g(6), unknown(0).
0329         KFSRCE=0
0330         IPART1=K(IPART(1),3)
0331         IPART2=K(IPART(2),3)
0332   220   IF(IPART1.EQ.IPART2.AND.IPART1.GT.0) THEN
0333           KFSRCE=IABS(K(IPART1,2))
0334         ELSEIF(IPART1.GT.IPART2.AND.IPART2.GT.0) THEN
0335           IPART1=K(IPART1,3)
0336           GOTO 220
0337         ELSEIF(IPART2.GT.IPART1.AND.IPART1.GT.0) THEN
0338           IPART2=K(IPART2,3)
0339           GOTO 220
0340         ENDIF
0341         ITYPES=0
0342         IF(KFSRCE.GE.1.AND.KFSRCE.LE.8) ITYPES=1
0343         IF(KFSRCE.GE.KSUSY1+1.AND.KFSRCE.LE.KSUSY1+8) ITYPES=2
0344         IF(KFSRCE.GE.KSUSY2+1.AND.KFSRCE.LE.KSUSY2+8) ITYPES=2
0345         IF(KFSRCE.GE.21.AND.KFSRCE.LE.24) ITYPES=3
0346         IF(KFSRCE.GE.32.AND.KFSRCE.LE.34) ITYPES=3
0347         IF(KFSRCE.EQ.25.OR.(KFSRCE.GE.35.AND.KFSRCE.LE.37)) ITYPES=4
0348         IF(KFSRCE.GE.KSUSY1+22.AND.KFSRCE.LE.KSUSY1+37) ITYPES=5
0349         IF(KFSRCE.EQ.KSUSY1+21) ITYPES=6
0350
0351 C...Identify two primary showerers.
0352         KFLA1=IABS(K(IPART(1),2))
0353         ITYPE1=0
0354         IF(KFLA1.GE.1.AND.KFLA1.LE.8) ITYPE1=1
0355         IF(KFLA1.GE.KSUSY1+1.AND.KFLA1.LE.KSUSY1+8) ITYPE1=2
0356         IF(KFLA1.GE.KSUSY2+1.AND.KFLA1.LE.KSUSY2+8) ITYPE1=2
0357         IF(KFLA1.GE.21.AND.KFLA1.LE.24) ITYPE1=3
0358         IF(KFLA1.GE.32.AND.KFLA1.LE.34) ITYPE1=3
0359         IF(KFLA1.EQ.25.OR.(KFLA1.GE.35.AND.KFLA1.LE.37)) ITYPE1=4
0360         IF(KFLA1.GE.KSUSY1+22.AND.KFLA1.LE.KSUSY1+37) ITYPE1=5
0361         IF(KFLA1.EQ.KSUSY1+21) ITYPE1=6
0362         KFLA2=IABS(K(IPART(2),2))
0363         ITYPE2=0
0364         IF(KFLA2.GE.1.AND.KFLA2.LE.8) ITYPE2=1
0365         IF(KFLA2.GE.KSUSY1+1.AND.KFLA2.LE.KSUSY1+8) ITYPE2=2
0366         IF(KFLA2.GE.KSUSY2+1.AND.KFLA2.LE.KSUSY2+8) ITYPE2=2
0367         IF(KFLA2.GE.21.AND.KFLA2.LE.24) ITYPE2=3
0368         IF(KFLA2.GE.32.AND.KFLA2.LE.34) ITYPE2=3
0369         IF(KFLA2.EQ.25.OR.(KFLA2.GE.35.AND.KFLA2.LE.37)) ITYPE2=4
0370         IF(KFLA2.GE.KSUSY1+22.AND.KFLA2.LE.KSUSY1+37) ITYPE2=5
0371         IF(KFLA2.EQ.KSUSY1+21) ITYPE2=6
0372
0373 C...Order of showerers. Presence of gluino.
0374         ITYPMN=MIN(ITYPE1,ITYPE2)
0375         ITYPMX=MAX(ITYPE1,ITYPE2)
0376         IORD=1
0377         IF(ITYPE1.GT.ITYPE2) IORD=2
0378         IGLUI=0
0379         IF(ITYPE1.EQ.6.OR.ITYPE2.EQ.6) IGLUI=1
0380
0381 C...Require exactly two primary showerers for ME corrections.
0382         NPRIM=0
0383         DO 230 I=1,N
0384           IF(K(I,3).EQ.IPART1.AND.K(I,2).NE.K(IPART1,2)) NPRIM=NPRIM+1
0385   230   CONTINUE
0386         IF(NPRIM.NE.2) THEN
0387
0388 C...Predetermined and default matrix element kinds.
0389         ELSEIF(MSTJ(38).NE.0) THEN
0390           M3JC=MSTJ(38)
0391           ALPHA=PARJ(80)
0392           MSTJ(38)=0
0393         ELSEIF(MSTJ(47).GE.6) THEN
0394           M3JC=MSTJ(47)
0395         ELSE
0396           ICLASS=1
0397           ICOMBI=4
0398
0399 C...Vector/axial vector -> q + qbar; q -> q + V.
0400           IF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.(ITYPES.EQ.0.OR.
0401      &    ITYPES.EQ.3)) THEN
0402             ICLASS=2
0403             IF(KFSRCE.EQ.21.OR.KFSRCE.EQ.22) THEN
0404               ICOMBI=1
0405             ELSEIF(KFSRCE.EQ.23.OR.(KFSRCE.EQ.0.AND.
0406      &      K(IPART(1),2)+K(IPART(2),2).EQ.0)) THEN
0407 C...gamma*/Z0: assume e+e- initial state if unknown.
0408               EI=-1D0
0409               IF(KFSRCE.EQ.23) THEN
0410                 IANNFL=IPART1
0411                 IF(K(IANNFL,2).EQ.23) IANNFL=K(IANNFL,3)
0412                 IF(IANNFL.GT.0) THEN
0413                   IF(K(IANNFL,2).EQ.23) IANNFL=K(IANNFL,3)
0414                 ENDIF
0415                 IF(IANNFL.NE.0) THEN
0416                   KANNFL=IABS(K(IANNFL,2))
0417                   IF(KANNFL.GE.1.AND.KANNFL.LE.18) EI=KCHG(KANNFL,1)/3D0
0418                 ENDIF
0419               ENDIF
0420               AI=SIGN(1D0,EI+0.1D0)
0421               VI=AI-4D0*EI*PARU(102)
0422               EF=KCHG(KFLA1,1)/3D0
0423               AF=SIGN(1D0,EF+0.1D0)
0424               VF=AF-4D0*EF*PARU(102)
0425               XWC=1D0/(16D0*PARU(102)*(1D0-PARU(102)))
0426               SH=PSUM(5)**2
0427               SQMZ=PMAS(23,1)**2
0428               SQWZ=PSUM(5)*PMAS(23,2)
0429               SBWZ=1D0/((SH-SQMZ)**2+SQWZ**2)
0430               VECT=EI**2*EF**2+2D0*EI*VI*EF*VF*XWC*SH*(SH-SQMZ)*SBWZ+
0431      &        (VI**2+AI**2)*VF**2*XWC**2*SH**2*SBWZ
0432               AXIV=(VI**2+AI**2)*AF**2*XWC**2*SH**2*SBWZ
0433               ICOMBI=3
0434               ALPHA=VECT/(VECT+AXIV)
0435             ELSEIF(KFSRCE.EQ.24.OR.KFSRCE.EQ.0) THEN
0436               ICOMBI=4
0437             ENDIF
0438 C...For chi -> chi q qbar, use V/A -> q qbar as first approximation.
0439           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.ITYPES.EQ.5) THEN
0440             ICLASS=2
0441           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.3.AND.(ITYPES.EQ.0.OR.
0442      &    ITYPES.EQ.1)) THEN
0443             ICLASS=3
0444
0445 C...Scalar/pseudoscalar -> q + qbar; q -> q + S.
0446           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.ITYPES.EQ.4) THEN
0447             ICLASS=4
0448             IF(KFSRCE.EQ.25.OR.KFSRCE.EQ.35.OR.KFSRCE.EQ.37) THEN
0449               ICOMBI=1
0450             ELSEIF(KFSRCE.EQ.36) THEN
0451               ICOMBI=2
0452             ENDIF
0453           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.4.AND.(ITYPES.EQ.0.OR.
0454      &    ITYPES.EQ.1)) THEN
0455             ICLASS=5
0456
0457 C...V -> ~q + ~qbar; ~q -> ~q + V; S -> ~q + ~qbar; ~q -> ~q + S.
0458           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.2.AND.(ITYPES.EQ.0.OR.
0459      &    ITYPES.EQ.3)) THEN
0460             ICLASS=6
0461           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.3.AND.(ITYPES.EQ.0.OR.
0462      &    ITYPES.EQ.2)) THEN
0463             ICLASS=7
0464           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.2.AND.ITYPES.EQ.4) THEN
0465             ICLASS=8
0466           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.4.AND.(ITYPES.EQ.0.OR.
0467      &    ITYPES.EQ.2)) THEN
0468             ICLASS=9
0469
0470 C...chi -> q + ~qbar; ~q -> q + chi; q -> ~q + chi.
0471           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.2.AND.(ITYPES.EQ.0.OR.
0472      &    ITYPES.EQ.5)) THEN
0473             ICLASS=10
0474           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.5.AND.(ITYPES.EQ.0.OR.
0475      &    ITYPES.EQ.2)) THEN
0476             ICLASS=11
0477           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.5.AND.(ITYPES.EQ.0.OR.
0478      &    ITYPES.EQ.1)) THEN
0479             ICLASS=12
0480
0481 C...~g -> q + ~qbar; ~q -> q + ~g; q -> ~q + ~g.
0482           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.2.AND.ITYPES.EQ.6) THEN
0483             ICLASS=13
0484           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.6.AND.(ITYPES.EQ.0.OR.
0485      &    ITYPES.EQ.2)) THEN
0486             ICLASS=14
0487           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.6.AND.(ITYPES.EQ.0.OR.
0488      &    ITYPES.EQ.1)) THEN
0489             ICLASS=15
0490
0491 C...g -> ~g + ~g (eikonal approximation).
0492           ELSEIF(ITYPMN.EQ.6.AND.ITYPMX.EQ.6.AND.ITYPES.EQ.0) THEN
0493             ICLASS=16
0494           ENDIF
0495           M3JC=5*ICLASS+ICOMBI
0496         ENDIF
0497
0498 C...Store pair that together define matrix element treatment.
0499         IF(M3JC.NE.0) THEN
0500           NMESYS=1
0501           MESYS(NMESYS,0)=M3JC
0502           MESYS(NMESYS,1)=IPART(1)
0503           MESYS(NMESYS,2)=IPART(2)
0504         ENDIF
0505
0506 C...Store qqbar or l+l- pairs for QED radiation.
0507         IF(KFLA1.LE.18.AND.KFLA2.LE.18) THEN
0508           NMESYS=NMESYS+1
0509           MESYS(NMESYS,0)=101
0510           IF(K(IPART(1),2)+K(IPART(2),2).EQ.0) MESYS(NMESYS,0)=102
0511           MESYS(NMESYS,1)=IPART(1)
0512           MESYS(NMESYS,2)=IPART(2)
0513         ENDIF
0514
0515 C...Store other qqbar/l+l- pairs from g/gamma branchings.
0516         DO 270 I1=1,N
0517           IF(K(I1,1).GT.10.OR.IABS(K(I1,2)).GT.18) GOTO 270
0518           I1M=K(I1,3)
0519   240     IF(I1M.GT.0.AND.K(I1M,2).EQ.K(I1,2)) THEN
0520             I1M=K(I1M,3)
0521             GOTO 240
0522           ENDIF
0523 C...Move up this check to avoid out-of-bounds.
0524           IF(I1M.EQ.0) GOTO 270
0525           IF(K(I1M,2).NE.21.AND.K(I1M,2).NE.22) GOTO 270
0526           DO 260 I2=I1+1,N
0527             IF(K(I2,1).GT.10.OR.K(I2,2)+K(I1,2).NE.0) GOTO 260
0528             I2M=K(I2,3)
0529   250       IF(I2M.GT.0.AND.K(I2M,2).EQ.K(I2,2)) THEN
0530               I2M=K(I2M,3)
0531               GOTO 250
0532             ENDIF
0533             IF(I1M.EQ.I2M.AND.I1M.GT.0) THEN
0534               NMESYS=NMESYS+1
0535               MESYS(NMESYS,0)=66
0536               MESYS(NMESYS,1)=I1
0537               MESYS(NMESYS,2)=I2
0538               NMESYS=NMESYS+1
0539               MESYS(NMESYS,0)=102
0540               MESYS(NMESYS,1)=I1
0541               MESYS(NMESYS,2)=I2
0542             ENDIF
0543   260     CONTINUE
0544   270   CONTINUE
0545       ENDIF
0546
0547 C..Loopback point for counting number of emissions.
0548       NGEN=0
0549   280 NGEN=NGEN+1
0550
0551 C...Begin loop to evolve all existing partons, if required.
0552   290 IMX=0
0553       PT2MX=0D0
0554       DO 360 IEVOL=1,NEVOL
0555         IF(IFLG(IEVOL).EQ.0) THEN
0556
0557 C...Basic info on radiator and recoil.
0558           I=IPOS(IEVOL)
0559           IR=IREC(IEVOL)
0560           SHT=SHAT(I,IR)
0561           PM2I=P(I,5)**2
0562           PM2R=P(IR,5)**2
0563
0564 C...Invariant mass of "dipole".Starting value for pT evolution.
0565           SHTCOR=(SQRT(SHT)-P(IR,5))**2-PM2I
0566           PT2=MIN(PT2CMX,0.25D0*SHTCOR,PTSCA(IEVOL)**2)
0567
0568 C...Case of evolution by QCD branching.
0569           IF(ISCOL(IEVOL).NE.0) THEN
0570
0571 C...Parton-by-parton maximum scale from initial conditions.
0572           IF(MSTP(72).EQ.0) THEN
0573             DO 300 IPRT=1,NPARTS
0574               IF(IR.EQ.IPART(IPRT)) PT2=MIN(PT2,PTPART(IPRT)**2)
0575   300       CONTINUE
0576           ENDIF
0577
0578 C...If kinematically impossible then do not evolve.
0579             IF(PT2.LT.PT2CMN) THEN
0580               IFLG(IEVOL)=-1
0581               GOTO 360
0582             ENDIF
0583
0584 C...Check if part of system for which ME corrections should be applied.
0585             IMESYS=0
0586             DO 310 IME=1,NMESYS
0587               IF((I.EQ.MESYS(IME,1).OR.I.EQ.MESYS(IME,2)).AND.
0588      &        MESYS(IME,0).LT.100) IMESYS=IME
0589   310       CONTINUE
0590
0591 C...Special flag for colour octet states.
0592             MOCT=0
0593             IF(K(I,2).EQ.21) MOCT=1
0594             IF(K(I,2).EQ.KSUSY1+21) MOCT=2
0595
0596 C...Upper estimate for matrix element weighting and colour factor.
0597 C...Note that g->gg and g->qqbar is split on two sides = "dipoles".
0598             WTPSGL=2D0
0599             COLFAC=4D0/3D0
0600             IF(MOCT.GE.1) COLFAC=3D0/2D0
0601             IF(IGLUI.EQ.1.AND.IMESYS.EQ.1.AND.MOCT.EQ.0) COLFAC=3D0
0602             WTPSQQ=0.5D0*0.5D0*NFLAV
0603
0604 C...Determine overestimated z range: switch at c and b masses.
0605   320       IZRG=1
0606             PT2MNE=PT2CMN
0607             B0=27D0/6D0
0608             ALAMS=ALAM3S
0609             IF(PT2.GT.1.01D0*PMCS) THEN
0610               IZRG=2
0611               PT2MNE=PMCS
0612               B0=25D0/6D0
0613               ALAMS=ALAM4S
0614             ENDIF
0615             IF(PT2.GT.1.01D0*PMBS) THEN
0616               IZRG=3
0617               PT2MNE=PMBS
0618               B0=23D0/6D0
0619               ALAMS=ALAM5S
0620             ENDIF
0621             ZMNCUT=0.5D0-SQRT(MAX(0D0,0.25D0-PT2MNE/SHTCOR))
0622             IF(ZMNCUT.LT.1D-8) ZMNCUT=PT2MNE/SHTCOR
0623
0624 C...Find evolution coefficients for q->qg/g->gg and g->qqbar.
0625             EVEMGL=WTPSGL*COLFAC*LOG(1D0/ZMNCUT-1D0)/B0
0626             EVCOEF=EVEMGL
0627             IF(MOCT.EQ.1) THEN
0628               EVEMQQ=WTPSQQ*(1D0-2D0*ZMNCUT)/B0
0629               EVCOEF=EVCOEF+EVEMQQ
0630             ENDIF
0631
0632 C...Pick pT2 (in overestimated z range).
0633   330       PT2=ALAMS*(PT2/ALAMS)**(PYR(0)**(1D0/EVCOEF))
0634
0635 C...Loopback if crossed c/b mass thresholds.
0636             IF(IZRG.EQ.3.AND.PT2.LT.PMBS) THEN
0637               PT2=PMBS
0638               GOTO 320
0639             ENDIF
0640             IF(IZRG.EQ.2.AND.PT2.LT.PMCS) THEN
0641               PT2=PMCS
0642               GOTO 320
0643             ENDIF
0644
0645 C...Finish if below lower cutoff.
0646             IF(PT2.LT.PT2CMN) THEN
0647               IFLG(IEVOL)=-1
0648               GOTO 360
0649             ENDIF
0650
0651 C...Pick kind of branching: q->qg/g->gg/X->Xg or g->qqbar.
0652             IFLAG=1
0653             IF(MOCT.EQ.1.AND.EVEMGL.LT.PYR(0)*EVCOEF) IFLAG=2
0654
0655 C...Pick z: dz/(1-z) or dz.
0656             IF(IFLAG.EQ.1) THEN
0657               Z=1D0-ZMNCUT*(1D0/ZMNCUT-1D0)**PYR(0)
0658             ELSE
0659               Z=ZMNCUT+PYR(0)*(1D0-2D0*ZMNCUT)
0660             ENDIF
0661
0662 C...Loopback if outside allowed range for given pT2.
0663             ZMNNOW=0.5D0-SQRT(MAX(0D0,0.25D0-PT2/SHTCOR))
0664             IF(ZMNNOW.LT.1D-8) ZMNNOW=PT2/SHTCOR
0665             IF(Z.LE.ZMNNOW.OR.Z.GE.1D0-ZMNNOW) GOTO 330
0666             PM2=PM2I+PT2/(Z*(1D0-Z))
0667             IF(Z*(1D0-Z).LE.PM2*SHT/(SHT+PM2-PM2R)**2) GOTO 330
0668
0669 C...No weighting for primary partons; to be done later on.
0670             IF(IMESYS.GT.0) THEN
0671
0672 C...Weighting of q->qg/X->Xg branching.
0673             ELSEIF(IFLAG.EQ.1.AND.MOCT.NE.1) THEN
0674               IF(1D0+Z**2.LT.WTPSGL*PYR(0)) GOTO 330
0675
0676 C...Weighting of g->gg branching.
0677             ELSEIF(IFLAG.EQ.1) THEN
0678               IF(1D0+Z**3.LT.WTPSGL*PYR(0)) GOTO 330
0679
0680 C...Flavour choice and weighting of g->qqbar branching.
0681             ELSE
0682               KFQ=MIN(5,1+INT(NFLAV*PYR(0)))
0683               PMQ=PMAS(KFQ,1)
0684               ROOTQQ=SQRT(MAX(0D0,1D0-4D0*PMQ**2/PM2))
0685               WTME=ROOTQQ*(Z**2+(1D0-Z)**2)
0686               IF(WTME.LT.PYR(0)) GOTO 330
0687               IFLAG=10+KFQ
0688             ENDIF
0689
0690 C...Case of evolution by QED branching.
0691           ELSEIF(ISCHG(IEVOL).NE.0) THEN
0692
0693 C...If kinematically impossible then do not evolve.
0694             PT2EMN=PT0EQ**2
0695             IF(IABS(K(I,2)).GT.10) PT2EMN=PT0EL**2
0696             IF(PT2.LT.PT2EMN) THEN
0697               IFLG(IEVOL)=-1
0698               GOTO 360
0699             ENDIF
0700
0701 C...Check if part of system for which ME corrections should be applied.
0702            IMESYS=0
0703             DO 340 IME=1,NMESYS
0704               IF((I.EQ.MESYS(IME,1).OR.I.EQ.MESYS(IME,2)).AND.
0705      &        MESYS(IME,0).GT.100) IMESYS=IME
0706   340      CONTINUE
0707
0708 C...Charge. Matrix element weighting factor.
0709             CHG=ISCHG(IEVOL)/3D0
0710             WTPSGA=2D0
0711
0712 C...Determine overestimated z range. Find evolution coefficient.
0713             ZMNCUT=0.5D0-SQRT(MAX(0D0,0.25D0-PT2EMN/SHTCOR))
0714             IF(ZMNCUT.LT.1D-8) ZMNCUT=PT2EMN/SHTCOR
0715             EVCOEF=AEM2PI*CHG**2*WTPSGA*LOG(1D0/ZMNCUT-1D0)
0716
0717 C...Pick pT2 (in overestimated z range).
0718   350       PT2=PT2*PYR(0)**(1D0/EVCOEF)
0719
0720 C...Finish if below lower cutoff.
0721             IF(PT2.LT.PT2EMN) THEN
0722               IFLG(IEVOL)=-1
0723               GOTO 360
0724             ENDIF
0725
0726 C...Pick z: dz/(1-z).
0727             Z=1D0-ZMNCUT*(1D0/ZMNCUT-1D0)**PYR(0)
0728
0729 C...Loopback if outside allowed range for given pT2.
0730             ZMNNOW=0.5D0-SQRT(MAX(0D0,0.25D0-PT2/SHTCOR))
0731             IF(ZMNNOW.LT.1D-8) ZMNNOW=PT2/SHTCOR
0732             IF(Z.LE.ZMNNOW.OR.Z.GE.1D0-ZMNNOW) GOTO 350
0733             PM2=PM2I+PT2/(Z*(1D0-Z))
0734             IF(Z*(1D0-Z).LE.PM2*SHT/(SHT+PM2-PM2R)**2) GOTO 350
0735
0736 C...Weighting by branching kernel, except if ME weighting later.
0737             IF(IMESYS.EQ.0) THEN
0738               IF(1D0+Z**2.LT.WTPSGA*PYR(0)) GOTO 350
0739             ENDIF
0740             IFLAG=3
0741           ENDIF
0742
0743 C...Save acceptable branching.
0744           IFLG(IEVOL)=IFLAG
0745           IMESAV(IEVOL)=IMESYS
0746           PT2SAV(IEVOL)=PT2
0747           ZSAV(IEVOL)=Z
0748           SHTSAV(IEVOL)=SHT
0749         ENDIF
0750
0751 C...Check if branching has highest pT.
0752         IF(IFLG(IEVOL).GE.1.AND.PT2SAV(IEVOL).GT.PT2MX) THEN
0753           IMX=IEVOL
0754           PT2MX=PT2SAV(IEVOL)
0755         ENDIF
0756   360 CONTINUE
0757
0758 C...Finished if no more branchings to be done.
0759       IF(IMX.EQ.0) GOTO 480
0760
0761 C...Restore info on hardest branching to be processed.
0762       I=IPOS(IMX)
0763       IR=IREC(IMX)
0764       KCOL=ISCOL(IMX)
0765       KCHA=ISCHG(IMX)
0766       IMESYS=IMESAV(IMX)
0767       PT2=PT2SAV(IMX)
0768       Z=ZSAV(IMX)
0769       SHT=SHTSAV(IMX)
0770       PM2I=P(I,5)**2
0771       PM2R=P(IR,5)**2
0772       PM2=PM2I+PT2/(Z*(1D0-Z))
0773
0774 C...Special flag for colour octet states.
0775       MOCT=0
0776       IF(K(I,2).EQ.21) MOCT=1
0777       IF(K(I,2).EQ.KSUSY1+21) MOCT=2
0778
0779 C...Restore further info for g->qqbar branching.
0780       KFQ=0
0781       IF(IFLG(IMX).GT.10) THEN
0782         KFQ=IFLG(IMX)-10
0783         PMQ=PMAS(KFQ,1)
0784         ROOTQQ=SQRT(MAX(0D0,1D0-4D0*PMQ**2/PM2))
0785       ENDIF
0786
0787 C...For branching g include azimuthal asymmetries from polarization.
0788       ASYPOL=0D0
0789       IF(MOCT.EQ.1.AND.MOD(MSTJ(46),2).EQ.1) THEN
0790 C...Trace grandmother via intermediate recoil copies.
0791         KFGM=0
0792         IM=I
0793   370   IF(K(IM,3).NE.K(IM-1,3).AND.K(IM,3).NE.K(IM+1,3).AND.
0794      &  K(IM,3).GT.0) THEN
0795           IM=K(IM,3)
0796           IF(IM.GT.MINT(84)) GOTO 370
0797         ENDIF
0798         IGM=K(IM,3)
0799         IF(IGM.GT.MINT(84).AND.IGM.LT.IM.AND.IM.LE.I)
0800      &  KFGM=IABS(K(IGM,2))
0801 C...Define approximate energy sharing by identifying aunt.
0802         IAU=IM+1
0803         IF(IAU.GT.N-3.OR.K(IAU,3).NE.IGM) IAU=IM-1
0804         IF(KFGM.NE.0.AND.(KFGM.LE.6.OR.KFGM.EQ.21)) THEN
0805           ZOLD=P(IM,4)/(P(IM,4)+P(IAU,4))
0806 C...Coefficient from gluon production.
0807           IF(KFGM.LE.6) THEN
0808             ASYPOL=2D0*(1D0-ZOLD)/(1D0+(1D0-ZOLD)**2)
0809           ELSE
0810             ASYPOL=((1D0-ZOLD)/(1D0-ZOLD*(1D0-ZOLD)))**2
0811           ENDIF
0812 C...Coefficient from gluon decay.
0813           IF(KFQ.EQ.0) THEN
0814             ASYPOL=ASYPOL*(Z*(1D0-Z)/(1D0-Z*(1D0-Z)))**2
0815           ELSE
0816             ASYPOL=-ASYPOL*2D0*Z*(1D0-Z)/(1D0-2D0*Z*(1D0-Z))
0817           ENDIF
0818         ENDIF
0819       ENDIF
0820
0821 C...Create new slots for branching products and recoil.
0822       INEW=N+1
0823       IGNEW=N+2
0824       IRNEW=N+3
0825       N=N+3
0826
0827 C...Set status, flavour and mother of new ones.
0828       K(INEW,1)=K(I,1)
0829       K(IGNEW,1)=3
0830       IF(KCHA.NE.0)  K(IGNEW,1)=1
0831       K(IRNEW,1)=K(IR,1)
0832       IF(KFQ.EQ.0) THEN
0833         K(INEW,2)=K(I,2)
0834         K(IGNEW,2)=21
0835         IF(KCHA.NE.0)  K(IGNEW,2)=22
0836       ELSE
0837         K(INEW,2)=-ISIGN(KFQ,KCOL)
0838         K(IGNEW,2)=-K(INEW,2)
0839       ENDIF
0840       K(IRNEW,2)=K(IR,2)
0841       K(INEW,3)=I
0842       K(IGNEW,3)=I
0843       K(IRNEW,3)=IR
0844
0845 C...Find rest frame and angles of branching+recoil.
0846       DO 380 J=1,5
0847         P(INEW,J)=P(I,J)
0848         P(IGNEW,J)=0D0
0849         P(IRNEW,J)=P(IR,J)
0850   380 CONTINUE
0851       BETAX=(P(INEW,1)+P(IRNEW,1))/(P(INEW,4)+P(IRNEW,4))
0852       BETAY=(P(INEW,2)+P(IRNEW,2))/(P(INEW,4)+P(IRNEW,4))
0853       BETAZ=(P(INEW,3)+P(IRNEW,3))/(P(INEW,4)+P(IRNEW,4))
0854       CALL PYROBO(INEW,IRNEW,0D0,0D0,-BETAX,-BETAY,-BETAZ)
0855       PHI=PYANGL(P(INEW,1),P(INEW,2))
0856       THETA=PYANGL(P(INEW,3),SQRT(P(INEW,1)**2+P(INEW,2)**2))
0857
0858 C...Derive kinematics of branching: generics (like g->gg).
0859       DO 390 J=1,4
0860         P(INEW,J)=0D0
0861         P(IRNEW,J)=0D0
0862   390 CONTINUE
0863       PEM=0.5D0*(SHT+PM2-PM2R)/SQRT(SHT)
0864       PZM=0.5D0*SQRT(MAX(0D0,(SHT-PM2-PM2R)**2-4D0*PM2*PM2R)/SHT)
0865       PT2COR=PM2*(PEM**2*Z*(1D0-Z)-0.25D0*PM2)/PZM**2
0866       PTCOR=SQRT(MAX(0D0,PT2COR))
0867       PZN=(PEM**2*Z-0.5D0*PM2)/PZM
0868       PZG=(PEM**2*(1D0-Z)-0.5D0*PM2)/PZM
0869 C...Specific kinematics reduction for q->qg with m_q > 0.
0870       IF(MOCT.NE.1) THEN
0871         PTCOR=(1D0-PM2I/PM2)*PTCOR
0872         PZN=PZN+PM2I*PZG/PM2
0873         PZG=(1D0-PM2I/PM2)*PZG
0874 C...Specific kinematics reduction for g->qqbar with m_q > 0.
0875       ELSEIF(KFQ.NE.0) THEN
0876         P(INEW,5)=PMQ
0877         P(IGNEW,5)=PMQ
0878         PTCOR=ROOTQQ*PTCOR
0879         PZN=0.5D0*((1D0+ROOTQQ)*PZN+(1D0-ROOTQQ)*PZG)
0880         PZG=PZM-PZN
0881       ENDIF
0882
0883 C...Pick phi and construct kinematics of branching.
0884   400 PHIROT=PARU(2)*PYR(0)
0885       P(INEW,1)=PTCOR*COS(PHIROT)
0886       P(INEW,2)=PTCOR*SIN(PHIROT)
0887       P(INEW,3)=PZN
0888       P(INEW,4)=SQRT(PTCOR**2+P(INEW,3)**2+P(INEW,5)**2)
0889       P(IGNEW,1)=-P(INEW,1)
0890       P(IGNEW,2)=-P(INEW,2)
0891       P(IGNEW,3)=PZG
0892       P(IGNEW,4)=SQRT(PTCOR**2+P(IGNEW,3)**2+P(IGNEW,5)**2)
0893       P(IRNEW,1)=0D0
0894       P(IRNEW,2)=0D0
0895       P(IRNEW,3)=-PZM
0896       P(IRNEW,4)=0.5D0*(SHT+PM2R-PM2)/SQRT(SHT)
0897
0898 C...Boost branching system to lab frame.
0899       CALL PYROBO(INEW,IRNEW,THETA,PHI,BETAX,BETAY,BETAZ)
0900
0901 C...Renew choice of phi angle according to polarization asymmetry.
0902       IF(ABS(ASYPOL).GT.1D-3) THEN
0903         DO 410 J=1,3
0904           DPT(1,J)=P(I,J)
0905           DPT(2,J)=P(IAU,J)
0906           DPT(3,J)=P(INEW,J)
0907   410   CONTINUE
0908         DPMA=DPT(1,1)*DPT(2,1)+DPT(1,2)*DPT(2,2)+DPT(1,3)*DPT(2,3)
0909         DPMD=DPT(1,1)*DPT(3,1)+DPT(1,2)*DPT(3,2)+DPT(1,3)*DPT(3,3)
0910         DPMM=DPT(1,1)**2+DPT(1,2)**2+DPT(1,3)**2
0911         DO 420 J=1,3
0912           DPT(4,J)=DPT(2,J)-DPMA*DPT(1,J)/MAX(1D-10,DPMM)
0913           DPT(5,J)=DPT(3,J)-DPMD*DPT(1,J)/MAX(1D-10,DPMM)
0914   420   CONTINUE
0915         DPT(4,4)=SQRT(DPT(4,1)**2+DPT(4,2)**2+DPT(4,3)**2)
0916         DPT(5,4)=SQRT(DPT(5,1)**2+DPT(5,2)**2+DPT(5,3)**2)
0917         IF(MIN(DPT(4,4),DPT(5,4)).GT.0.1D0*PARJ(82)) THEN
0918           CAD=(DPT(4,1)*DPT(5,1)+DPT(4,2)*DPT(5,2)+
0919      &    DPT(4,3)*DPT(5,3))/(DPT(4,4)*DPT(5,4))
0920           IF(1D0+ASYPOL*(2D0*CAD**2-1D0).LT.PYR(0)*(1D0+ABS(ASYPOL)))
0921      &    GOTO 400
0922         ENDIF
0923       ENDIF
0924
0925 C...Matrix element corrections for primary partons when requested.
0926       IF(IMESYS.GT.0) THEN
0927         M3JC=MESYS(IMESYS,0)
0928
0929 C...Identify recoiling partner and set up three-body kinematics.
0930         IRP=MESYS(IMESYS,1)
0931         IF(IRP.EQ.I) IRP=MESYS(IMESYS,2)
0932         IF(IRP.EQ.IR) IRP=IRNEW
0933         DO 430 J=1,4
0934           PSUM(J)=P(INEW,J)+P(IRP,J)+P(IGNEW,J)
0935   430   CONTINUE
0936         PSUM(5)=SQRT(MAX(0D0,PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-
0937      &  PSUM(3)**2))
0938         X1=2D0*(PSUM(4)*P(INEW,4)-PSUM(1)*P(INEW,1)-PSUM(2)*P(INEW,2)-
0939      &  PSUM(3)*P(INEW,3))/PSUM(5)**2
0940         X2=2D0*(PSUM(4)*P(IRP,4)-PSUM(1)*P(IRP,1)-PSUM(2)*P(IRP,2)-
0941      &  PSUM(3)*P(IRP,3))/PSUM(5)**2
0942         X3=2D0-X1-X2
0943         R1ME=P(INEW,5)/PSUM(5)
0944         R2ME=P(IRP,5)/PSUM(5)
0945
0946 C...Matrix elements for gluon emission.
0947         IF(M3JC.LT.100) THEN
0948
0949 C...Call ME, with right order important for two inequivalent showerers.
0950           IF(MESYS(IMESYS,IORD).EQ.I) THEN
0951             WME=PYMAEL(M3JC,X1,X2,R1ME,R2ME,ALPHA)
0952           ELSE
0953             WME=PYMAEL(M3JC,X2,X1,R2ME,R1ME,ALPHA)
0954           ENDIF
0955
0956 C...Split up total ME when two radiating partons.
0957           ISPRAD=1
0958           IF((M3JC.GE.16.AND.M3JC.LE.19).OR.(M3JC.GE.26.AND.M3JC.LE.29)
0959      &    .OR.(M3JC.GE.36.AND.M3JC.LE.39).OR.(M3JC.GE.46.AND.M3JC.LE.49)
0960      &    .OR.(M3JC.GE.56.AND.M3JC.LE.64)) ISPRAD=0
0961           IF(ISPRAD.EQ.1) WME=WME*MAX(1D-10,1D0+R1ME**2-R2ME**2-X1)/
0962      &    MAX(1D-10,2D0-X1-X2)
0963
0964 C...Evaluate shower rate.
0965           WPS=2D0/(MAX(1D-10,2D0-X1-X2)*
0966      &    MAX(1D-10,1D0+R2ME**2-R1ME**2-X2))
0967           IF(IGLUI.EQ.1) WPS=(9D0/4D0)*WPS
0968
0969 C...Matrix elements for photon emission: still rather primitive.
0970         ELSE
0971
0972 C...For generic charge combination currently only massless expression.
0973           IF(M3JC.EQ.101) THEN
0974             CHG1=KCHG(PYCOMP(K(I,2)),1)*ISIGN(1,K(I,2))/3D0
0975             CHG2=KCHG(PYCOMP(K(IRP,2)),1)*ISIGN(1,K(IRP,2))/3D0
0976             WME=(CHG1*(1D0-X1)/X3-CHG2*(1D0-X2)/X3)**2*(X1**2+X2**2)
0977             WPS=2D0*(CHG1**2*(1D0-X1)/X3+CHG2**2*(1D0-X2)/X3)
0978
0979 C...For flavour neutral system assume vector source and include masses.
0980           ELSE
0981             WME=PYMAEL(11,X1,X2,R1ME,R2ME,0D0)*MAX(1D-10,
0982      &      1D0+R1ME**2-R2ME**2-X1)/MAX(1D-10,2D0-X1-X2)
0983             WPS=2D0/(MAX(1D-10,2D0-X1-X2)*
0984      &      MAX(1D-10,1D0+R2ME**2-R1ME**2-X2))
0985           ENDIF
0986         ENDIF
0987
0988 C...Perform weighting with W_ME/W_PS.
0989         IF(WME.LT.PYR(0)*WPS) THEN
0990           N=N-3
0991           IFLG(IMX)=0
0992           GOTO 290
0993         ENDIF
0994       ENDIF
0995
0996 C...Now for sure accepted branching. Save highest pT.
0997       IF(NGEN.EQ.1) PTGEN=SQRT(PT2)
0998
0999 C...Update status for obsolete ones. Bookkkep the moved original parton
1000 C...and new daughter (arbitrary choice for g->gg or g->qqbar).
1001 C...Do not bookkeep radiated photon, since it cannot radiate further.
1002       K(I,1)=K(I,1)+10
1003       K(IR,1)=K(IR,1)+10
1004       DO 440 IP=1,NPART
1005         IF(IPART(IP).EQ.I) IPART(IP)=INEW
1006         IF(IPART(IP).EQ.IR) IPART(IP)=IRNEW
1007   440 CONTINUE
1008       IF(KCHA.EQ.0) THEN
1009         NPART=NPART+1
1010         IPART(NPART)=IGNEW
1011       ENDIF
1012
1013 C...Initialize colour flow of branching.
1014 C...Use both old and new style colour tags for flexibility.
1015       K(INEW,4)=0
1016       K(IGNEW,4)=0
1017       K(INEW,5)=0
1018       K(IGNEW,5)=0
1019       JCOLP=4+(1-KCOL)/2
1020       JCOLN=9-JCOLP
1021       MCT(INEW,1)=0
1022       MCT(INEW,2)=0
1023       MCT(IGNEW,1)=0
1024       MCT(IGNEW,2)=0
1025       MCT(IRNEW,1)=0
1026       MCT(IRNEW,2)=0
1027
1028 C...Trivial colour flow for l->lgamma and q->qgamma.
1029       IF(IABS(KCHA).EQ.3) THEN
1030         K(I,4)=INEW
1031         K(I,5)=IGNEW
1032       ELSEIF(KCHA.NE.0) THEN
1033         IF(K(I,4).NE.0) THEN
1034           K(I,4)=K(I,4)+INEW
1035           K(INEW,4)=MSTU(5)*I
1036           MCT(INEW,1)=MCT(I,1)
1037         ENDIF
1038         IF(K(I,5).NE.0) THEN
1039           K(I,5)=K(I,5)+INEW
1040           K(INEW,5)=MSTU(5)*I
1041           MCT(INEW,2)=MCT(I,2)
1042         ENDIF
1043
1044 C...Set colour flow for q->qg and g->gg.
1045       ELSEIF(KFQ.EQ.0) THEN
1046         K(I,JCOLP)=K(I,JCOLP)+IGNEW
1047         K(IGNEW,JCOLP)=MSTU(5)*I
1048         K(INEW,JCOLP)=MSTU(5)*IGNEW
1049         K(IGNEW,JCOLN)=MSTU(5)*INEW
1050         MCT(IGNEW,JCOLP-3)=MCT(I,JCOLP-3)
1051         NCT=NCT+1
1052         MCT(INEW,JCOLP-3)=NCT
1053         MCT(IGNEW,JCOLN-3)=NCT
1054         IF(MOCT.GE.1) THEN
1055           K(I,JCOLN)=K(I,JCOLN)+INEW
1056           K(INEW,JCOLN)=MSTU(5)*I
1057           MCT(INEW,JCOLN-3)=MCT(I,JCOLN-3)
1058         ENDIF
1059
1060 C...Set colour flow for g->qqbar.
1061       ELSE
1062         K(I,JCOLN)=K(I,JCOLN)+INEW
1063         K(INEW,JCOLN)=MSTU(5)*I
1064         K(I,JCOLP)=K(I,JCOLP)+IGNEW
1065         K(IGNEW,JCOLP)=MSTU(5)*I
1066         MCT(INEW,JCOLN-3)=MCT(I,JCOLN-3)
1067         MCT(IGNEW,JCOLP-3)=MCT(I,JCOLP-3)
1068       ENDIF
1069
1070 C...Daughter info for colourless recoiling parton.
1071       IF(K(IR,4).EQ.0.AND.K(IR,5).EQ.0) THEN
1072         K(IR,4)=IRNEW
1073         K(IR,5)=IRNEW
1074         K(IRNEW,4)=0
1075         K(IRNEW,5)=0
1076
1077 C...Colour of recoiling parton sails through unchanged.
1078       ELSE
1079         IF(K(IR,4).NE.0) THEN
1080           K(IR,4)=K(IR,4)+IRNEW
1081           K(IRNEW,4)=MSTU(5)*IR
1082           MCT(IRNEW,1)=MCT(IR,1)
1083         ENDIF
1084         IF(K(IR,5).NE.0) THEN
1085           K(IR,5)=K(IR,5)+IRNEW
1086           K(IRNEW,5)=MSTU(5)*IR
1087           MCT(IRNEW,2)=MCT(IR,2)
1088         ENDIF
1089       ENDIF
1090
1091 C...Vertex information trivial.
1092       DO 450 J=1,5
1093         V(INEW,J)=V(I,J)
1094         V(IGNEW,J)=V(I,J)
1095         V(IRNEW,J)=V(IR,J)
1096   450 CONTINUE
1097
1098 C...Update list of old radiators.
1099         DO 460 IEVOL=1,NEVOL
1100           IF(IPOS(IEVOL).EQ.I.AND.IREC(IEVOL).EQ.IR) THEN
1101             IPOS(IEVOL)=INEW
1102             IF(KCOL.NE.0.AND.ISCOL(IEVOL).EQ.KCOL) IPOS(IEVOL)=IGNEW
1103             IREC(IEVOL)=IRNEW
1104             IFLG(IEVOL)=0
1105           ELSEIF(IPOS(IEVOL).EQ.I) THEN
1106             IPOS(IEVOL)=INEW
1107             IFLG(IEVOL)=0
1108           ELSEIF(IPOS(IEVOL).EQ.IR.AND.IREC(IEVOL).EQ.I) THEN
1109             IPOS(IEVOL)=IRNEW
1110             IREC(IEVOL)=INEW
1111             IF(KCOL.NE.0.AND.ISCOL(IEVOL).NE.KCOL) IREC(IEVOL)=IGNEW
1112             IFLG(IEVOL)=0
1113           ELSEIF(IPOS(IEVOL).EQ.IR) THEN
1114             IPOS(IEVOL)=IRNEW
1115             IFLG(IEVOL)=0
1116           ENDIF
1117 C...Update links of old connected partons.
1118           IF(IREC(IEVOL).EQ.I) THEN
1119             IREC(IEVOL)=INEW
1120             IFLG(IEVOL)=0
1121           ELSEIF(IREC(IEVOL).EQ.IR) THEN
1122             IREC(IEVOL)=IRNEW
1123             IFLG(IEVOL)=0
1124           ENDIF
1125   460   CONTINUE
1126
1127 C...q->qg or g->gg: create new gluon radiators.
1128       IF(KCOL.NE.0.AND.KFQ.EQ.0) THEN
1129         NEVOL=NEVOL+1
1130         IPOS(NEVOL)=INEW
1131         IREC(NEVOL)=IGNEW
1132         IFLG(NEVOL)=0
1133         ISCOL(NEVOL)=KCOL
1134         ISCHG(NEVOL)=0
1135         PTSCA(NEVOL)=SQRT(PT2)
1136         NEVOL=NEVOL+1
1137         IPOS(NEVOL)=IGNEW
1138         IREC(NEVOL)=INEW
1139         IFLG(NEVOL)=0
1140         ISCOL(NEVOL)=-KCOL
1141         ISCHG(NEVOL)=0
1142         PTSCA(NEVOL)=PTSCA(NEVOL-1)
1143       ENDIF
1144
1145 C...Update matrix elements parton list and add new for g/gamma->qqbar.
1146       DO 470 IME=1,NMESYS
1147         IF(MESYS(IME,1).EQ.I) MESYS(IME,1)=INEW
1148         IF(MESYS(IME,2).EQ.I) MESYS(IME,2)=INEW
1149         IF(MESYS(IME,1).EQ.IR) MESYS(IME,1)=IRNEW
1150         IF(MESYS(IME,2).EQ.IR) MESYS(IME,2)=IRNEW
1151   470 CONTINUE
1152       IF(KFQ.NE.0) THEN
1153         NMESYS=NMESYS+1
1154         MESYS(NMESYS,0)=66
1155         MESYS(NMESYS,1)=INEW
1156         MESYS(NMESYS,2)=IGNEW
1157         NMESYS=NMESYS+1
1158         MESYS(NMESYS,0)=102
1159         MESYS(NMESYS,1)=INEW
1160         MESYS(NMESYS,2)=IGNEW
1161       ENDIF
1162
1163 C...Global statistics.
1164       MINT(353)=MINT(353)+1
1165       VINT(353)=VINT(353)+PTCOR
1166       IF (MINT(353).EQ.1) VINT(358)=PTCOR
1167
1168 C...Loopback for more emissions if enough space.
1169       PT2CMX=PT2
1170       IF(NPART.LT.MAXNUR-1.AND.NEVOL.LT.2*MAXNUR-2.AND.
1171      &NMESYS.LT.MAXNUR-2.AND.N.LT.MSTU(4)-MSTU(32)-5) THEN
1172         GOTO 280
1173       ELSE
1174         CALL PYERRM(11,'(PYPTFS:) no more memory left for shower')
1175       ENDIF
1176
1177 C...Done.
1178   480 CONTINUE
1179
1180       RETURN
1181       END