Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:21:18 

0001  
0002 C*********************************************************************
0003  
0004 C...PYSHOW
0005 C...Generates timelike parton showers from given partons.
0006  
0007       SUBROUTINE PYSHOW(IP1,IP2,QMAX)
0008  
0009 C...Double precision and integer declarations.
0010       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0011       IMPLICIT INTEGER(I-N)
0012       INTEGER PYK,PYCHGE,PYCOMP
0013 C...Parameter statement to help give large particle numbers.
0014       PARAMETER (KSUSY1=1000000,KSUSY2=2000000,KTECHN=3000000,
0015      &KEXCIT=4000000,KDIMEN=5000000)
0016       PARAMETER (MAXNUR=1000)
0017 C...Commonblocks.
0018       COMMON/PYPART/NPART,NPARTD,IPART(MAXNUR),PTPART(MAXNUR)
0019       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0020       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0021       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0022       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0023       COMMON/PYINT1/MINT(400),VINT(400)
0024       SAVE /PYPART/,/PYJETS/,/PYDAT1/,/PYDAT2/,/PYPARS/,/PYINT1/
0025 C...Local arrays.
0026       DIMENSION PMTH(5,140),PS(5),PMA(100),PMSD(100),IEP(100),IPA(100),
0027      &KFLA(100),KFLD(100),KFL(100),ITRY(100),ISI(100),ISL(100),DP(100),
0028      &DPT(5,4),KSH(0:140),KCII(2),NIIS(2),IIIS(2,2),THEIIS(2,2),
0029      &PHIIIS(2,2),ISII(2),ISSET(2),ISCOL(0:140),ISCHG(0:140),
0030      &IREF(1000)
0031  
0032 C...Check that QMAX not too low.
0033       IF(MSTJ(41).LE.0) THEN
0034         RETURN
0035       ELSEIF(MSTJ(41).EQ.1.OR.MSTJ(41).EQ.11) THEN
0036         IF(QMAX.LE.PARJ(82).AND.IP2.GE.-80) RETURN
0037       ELSE
0038         IF(QMAX.LE.MIN(PARJ(82),PARJ(83),PARJ(90)).AND.IP2.GE.-80)
0039      &  RETURN
0040       ENDIF
0041  
0042 C...Store positions of shower initiating partons.
0043       MPSPD=0
0044       IF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.EQ.0) THEN
0045         NPA=1
0046         IPA(1)=IP1
0047       ELSEIF(MIN(IP1,IP2).GT.0.AND.MAX(IP1,IP2).LE.MIN(N,MSTU(4)-
0048      &  MSTU(32))) THEN
0049         NPA=2
0050         IPA(1)=IP1
0051         IPA(2)=IP2
0052       ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.IP2.LT.0
0053      &  .AND.IP2.GE.-80) THEN
0054         NPA=IABS(IP2)
0055         DO 100 I=1,NPA
0056           IPA(I)=IP1+I-1
0057   100   CONTINUE
0058       ELSEIF(IP1.GT.0.AND.IP1.LE.MIN(N,MSTU(4)-MSTU(32)).AND.
0059      &IP2.EQ.-100) THEN
0060         MPSPD=1
0061         NPA=2
0062         IPA(1)=IP1+6
0063         IPA(2)=IP1+7
0064       ELSE
0065         CALL PYERRM(12,
0066      &  '(PYSHOW:) failed to reconstruct showering system')
0067         IF(MSTU(21).GE.1) RETURN
0068       ENDIF
0069  
0070 C...Send off to PYPTFS for pT-ordered evolution if requested,
0071 C...if at least 2 partons, and without predefined shower branchings.
0072       IF((MSTJ(41).EQ.11.OR.MSTJ(41).EQ.12).AND.NPA.GE.2.AND.
0073      &MPSPD.EQ.0) THEN
0074         NPART=NPA
0075         DO 110 II=1,NPART
0076           IPART(II)=IPA(II)
0077           PTPART(II)=0.5D0*QMAX
0078   110   CONTINUE
0079         CALL PYPTFS(2,0.5D0*QMAX,0D0,PTGEN)
0080         RETURN
0081       ENDIF
0082  
0083 C...Initialization of cutoff masses etc.
0084       DO 120 IFL=0,40
0085         ISCOL(IFL)=0
0086         ISCHG(IFL)=0
0087         KSH(IFL)=0
0088   120 CONTINUE
0089       ISCOL(21)=1
0090       KSH(21)=1
0091       PMTH(1,21)=PYMASS(21)
0092       PMTH(2,21)=SQRT(PMTH(1,21)**2+0.25D0*PARJ(82)**2)
0093       PMTH(3,21)=2D0*PMTH(2,21)
0094       PMTH(4,21)=PMTH(3,21)
0095       PMTH(5,21)=PMTH(3,21)
0096       PMTH(1,22)=PYMASS(22)
0097       PMTH(2,22)=SQRT(PMTH(1,22)**2+0.25D0*PARJ(83)**2)
0098       PMTH(3,22)=2D0*PMTH(2,22)
0099       PMTH(4,22)=PMTH(3,22)
0100       PMTH(5,22)=PMTH(3,22)
0101       PMQTH1=PARJ(82)
0102       IF(MSTJ(41).GE.2) PMQTH1=MIN(PARJ(82),PARJ(83))
0103       PMQT1E=MIN(PMQTH1,PARJ(90))
0104       PMQTH2=PMTH(2,21)
0105       IF(MSTJ(41).GE.2) PMQTH2=MIN(PMTH(2,21),PMTH(2,22))
0106       PMQT2E=MIN(PMQTH2,0.5D0*PARJ(90))
0107       DO 130 IFL=1,5
0108         ISCOL(IFL)=1
0109         IF(MSTJ(41).GE.2) ISCHG(IFL)=1
0110         KSH(IFL)=1
0111         PMTH(1,IFL)=PYMASS(IFL)
0112         PMTH(2,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PMQTH1**2)
0113         PMTH(3,IFL)=PMTH(2,IFL)+PMQTH2
0114         PMTH(4,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PARJ(82)**2)+PMTH(2,21)
0115         PMTH(5,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PARJ(83)**2)+PMTH(2,22)
0116   130 CONTINUE
0117       DO 140 IFL=11,15,2
0118         IF(MSTJ(41).EQ.2.OR.MSTJ(41).GE.4) ISCHG(IFL)=1
0119         IF(MSTJ(41).EQ.2.OR.MSTJ(41).GE.4) KSH(IFL)=1
0120         PMTH(1,IFL)=PYMASS(IFL)
0121         PMTH(2,IFL)=SQRT(PMTH(1,IFL)**2+0.25D0*PARJ(90)**2)
0122         PMTH(3,IFL)=PMTH(2,IFL)+0.5D0*PARJ(90)
0123         PMTH(4,IFL)=PMTH(3,IFL)
0124         PMTH(5,IFL)=PMTH(3,IFL)
0125   140 CONTINUE
0126       PT2MIN=MAX(0.5D0*PARJ(82),1.1D0*PARJ(81))**2
0127       ALAMS=PARJ(81)**2
0128       ALFM=LOG(PT2MIN/ALAMS)
0129  
0130 C...Check on phase space available for emission.
0131       IREJ=0
0132       DO 150 J=1,5
0133         PS(J)=0D0
0134   150 CONTINUE
0135       PM=0D0
0136       KFLA(2)=0
0137       DO 170 I=1,NPA
0138         KFLA(I)=IABS(K(IPA(I),2))
0139         PMA(I)=P(IPA(I),5)
0140 C...Special cutoff masses for initial partons (may be a heavy quark,
0141 C...squark, ..., and need not be on the mass shell).
0142         IR=30+I
0143         IF(NPA.LE.1) IREF(I)=IR
0144         IF(NPA.GE.2) IREF(I+1)=IR
0145         ISCOL(IR)=0
0146         ISCHG(IR)=0
0147         KSH(IR)=0
0148         IF(KFLA(I).LE.8) THEN
0149           ISCOL(IR)=1
0150           IF(MSTJ(41).GE.2) ISCHG(IR)=1
0151         ELSEIF(KFLA(I).EQ.11.OR.KFLA(I).EQ.13.OR.KFLA(I).EQ.15.OR.
0152      &  KFLA(I).EQ.17) THEN
0153           IF(MSTJ(41).EQ.2.OR.MSTJ(41).GE.4) ISCHG(IR)=1
0154         ELSEIF(KFLA(I).EQ.21) THEN
0155           ISCOL(IR)=1
0156         ELSEIF((KFLA(I).GE.KSUSY1+1.AND.KFLA(I).LE.KSUSY1+8).OR.
0157      &  (KFLA(I).GE.KSUSY2+1.AND.KFLA(I).LE.KSUSY2+8)) THEN
0158           ISCOL(IR)=1
0159         ELSEIF(KFLA(I).EQ.KSUSY1+21) THEN
0160           ISCOL(IR)=1
0161 C...QUARKONIA+++
0162 C...same for QQ~[3S18]
0163         ELSEIF(MSTP(148).GE.1.AND.(KFLA(I).EQ.9900443.OR.
0164      &  KFLA(I).EQ.9900553)) THEN
0165           ISCOL(IR)=1
0166 C...QUARKONIA---
0167         ENDIF
0168         IF(ISCOL(IR).EQ.1.OR.ISCHG(IR).EQ.1) KSH(IR)=1
0169         PMTH(1,IR)=PMA(I)
0170         IF(ISCOL(IR).EQ.1.AND.ISCHG(IR).EQ.1) THEN
0171           PMTH(2,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PMQTH1**2)
0172           PMTH(3,IR)=PMTH(2,IR)+PMQTH2
0173           PMTH(4,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(82)**2)+PMTH(2,21)
0174           PMTH(5,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(83)**2)+PMTH(2,22)
0175         ELSEIF(ISCOL(IR).EQ.1) THEN
0176           PMTH(2,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(82)**2)
0177           PMTH(3,IR)=PMTH(2,IR)+0.5D0*PARJ(82)
0178           PMTH(4,IR)=PMTH(3,IR)
0179           PMTH(5,IR)=PMTH(3,IR)
0180         ELSEIF(ISCHG(IR).EQ.1) THEN
0181           PMTH(2,IR)=SQRT(PMTH(1,IR)**2+0.25D0*PARJ(90)**2)
0182           PMTH(3,IR)=PMTH(2,IR)+0.5D0*PARJ(90)
0183           PMTH(4,IR)=PMTH(3,IR)
0184           PMTH(5,IR)=PMTH(3,IR)
0185         ENDIF
0186         IF(KSH(IR).EQ.1) PMA(I)=PMTH(3,IR)
0187         PM=PM+PMA(I)
0188         IF(KSH(IR).EQ.0.OR.PMA(I).GT.10D0*QMAX) IREJ=IREJ+1
0189         DO 160 J=1,4
0190           PS(J)=PS(J)+P(IPA(I),J)
0191   160   CONTINUE
0192   170 CONTINUE
0193       IF(IREJ.EQ.NPA.AND.IP2.GE.-7) RETURN
0194       PS(5)=SQRT(MAX(0D0,PS(4)**2-PS(1)**2-PS(2)**2-PS(3)**2))
0195       IF(NPA.EQ.1) PS(5)=PS(4)
0196       IF(PS(5).LE.PM+PMQT1E) RETURN
0197  
0198 C...Identify source: q(1), ~q(2), V(3), S(4), chi(5), ~g(6), unknown(0).
0199       KFSRCE=0
0200       IF(IP2.LE.0) THEN
0201       ELSEIF(K(IP1,3).EQ.K(IP2,3).AND.K(IP1,3).GT.0) THEN
0202         KFSRCE=IABS(K(K(IP1,3),2))
0203       ELSE
0204         IPAR1=MAX(1,K(IP1,3))
0205         IPAR2=MAX(1,K(IP2,3))
0206         IF(K(IPAR1,3).EQ.K(IPAR2,3).AND.K(IPAR1,3).GT.0)
0207      &       KFSRCE=IABS(K(K(IPAR1,3),2))
0208       ENDIF
0209       ITYPES=0
0210       IF(KFSRCE.GE.1.AND.KFSRCE.LE.8) ITYPES=1
0211       IF(KFSRCE.GE.KSUSY1+1.AND.KFSRCE.LE.KSUSY1+8) ITYPES=2
0212       IF(KFSRCE.GE.KSUSY2+1.AND.KFSRCE.LE.KSUSY2+8) ITYPES=2
0213       IF(KFSRCE.GE.21.AND.KFSRCE.LE.24) ITYPES=3
0214       IF(KFSRCE.GE.32.AND.KFSRCE.LE.34) ITYPES=3
0215       IF(KFSRCE.EQ.25.OR.(KFSRCE.GE.35.AND.KFSRCE.LE.37)) ITYPES=4
0216       IF(KFSRCE.GE.KSUSY1+22.AND.KFSRCE.LE.KSUSY1+37) ITYPES=5
0217       IF(KFSRCE.EQ.KSUSY1+21) ITYPES=6
0218  
0219 C...Identify two primary showerers.
0220       ITYPE1=0
0221       IF(KFLA(1).GE.1.AND.KFLA(1).LE.8) ITYPE1=1
0222       IF(KFLA(1).GE.KSUSY1+1.AND.KFLA(1).LE.KSUSY1+8) ITYPE1=2
0223       IF(KFLA(1).GE.KSUSY2+1.AND.KFLA(1).LE.KSUSY2+8) ITYPE1=2
0224       IF(KFLA(1).GE.21.AND.KFLA(1).LE.24) ITYPE1=3
0225       IF(KFLA(1).GE.32.AND.KFLA(1).LE.34) ITYPE1=3
0226       IF(KFLA(1).EQ.25.OR.(KFLA(1).GE.35.AND.KFLA(1).LE.37)) ITYPE1=4
0227       IF(KFLA(1).GE.KSUSY1+22.AND.KFLA(1).LE.KSUSY1+37) ITYPE1=5
0228       IF(KFLA(1).EQ.KSUSY1+21) ITYPE1=6
0229       ITYPE2=0
0230       IF(KFLA(2).GE.1.AND.KFLA(2).LE.8) ITYPE2=1
0231       IF(KFLA(2).GE.KSUSY1+1.AND.KFLA(2).LE.KSUSY1+8) ITYPE2=2
0232       IF(KFLA(2).GE.KSUSY2+1.AND.KFLA(2).LE.KSUSY2+8) ITYPE2=2
0233       IF(KFLA(2).GE.21.AND.KFLA(2).LE.24) ITYPE2=3
0234       IF(KFLA(2).GE.32.AND.KFLA(2).LE.34) ITYPE2=3
0235       IF(KFLA(2).EQ.25.OR.(KFLA(2).GE.35.AND.KFLA(2).LE.37)) ITYPE2=4
0236       IF(KFLA(2).GE.KSUSY1+22.AND.KFLA(2).LE.KSUSY1+37) ITYPE2=5
0237       IF(KFLA(2).EQ.KSUSY1+21) ITYPE2=6
0238  
0239 C...Order of showerers. Presence of gluino.
0240       ITYPMN=MIN(ITYPE1,ITYPE2)
0241       ITYPMX=MAX(ITYPE1,ITYPE2)
0242       IORD=1
0243       IF(ITYPE1.GT.ITYPE2) IORD=2
0244       IGLUI=0
0245       IF(ITYPE1.EQ.6.OR.ITYPE2.EQ.6) IGLUI=1
0246  
0247 C...Check if 3-jet matrix elements to be used.
0248       M3JC=0
0249       ALPHA=0.5D0
0250       IF(NPA.EQ.2.AND.MSTJ(47).GE.1.AND.MPSPD.EQ.0) THEN
0251         IF(MSTJ(38).NE.0) THEN
0252           M3JC=MSTJ(38)
0253           ALPHA=PARJ(80)
0254           MSTJ(38)=0
0255         ELSEIF(MSTJ(47).GE.6) THEN
0256           M3JC=MSTJ(47)
0257         ELSE
0258           ICLASS=1
0259           ICOMBI=4
0260  
0261 C...Vector/axial vector -> q + qbar; q -> q + V.
0262           IF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.(ITYPES.EQ.0.OR.
0263      &    ITYPES.EQ.3)) THEN
0264             ICLASS=2
0265             IF(KFSRCE.EQ.21.OR.KFSRCE.EQ.22) THEN
0266               ICOMBI=1
0267             ELSEIF(KFSRCE.EQ.23.OR.(KFSRCE.EQ.0.AND.
0268      &      K(IPA(1),2)+K(IPA(2),2).EQ.0)) THEN
0269 C...gamma*/Z0: assume e+e- initial state if unknown.
0270               EI=-1D0
0271               IF(KFSRCE.EQ.23) THEN
0272                 IANNFL=K(K(IP1,3),3)
0273                 IF(IANNFL.NE.0) THEN
0274                   KANNFL=IABS(K(IANNFL,2))
0275                   IF(KANNFL.GE.1.AND.KANNFL.LE.18) EI=KCHG(KANNFL,1)/3D0
0276                 ENDIF
0277               ENDIF
0278               AI=SIGN(1D0,EI+0.1D0)
0279               VI=AI-4D0*EI*PARU(102)
0280               EF=KCHG(KFLA(1),1)/3D0
0281               AF=SIGN(1D0,EF+0.1D0)
0282               VF=AF-4D0*EF*PARU(102)
0283               XWC=1D0/(16D0*PARU(102)*(1D0-PARU(102)))
0284               SH=PS(5)**2
0285               SQMZ=PMAS(23,1)**2
0286               SQWZ=PS(5)*PMAS(23,2)
0287               SBWZ=1D0/((SH-SQMZ)**2+SQWZ**2)
0288               VECT=EI**2*EF**2+2D0*EI*VI*EF*VF*XWC*SH*(SH-SQMZ)*SBWZ+
0289      &        (VI**2+AI**2)*VF**2*XWC**2*SH**2*SBWZ
0290               AXIV=(VI**2+AI**2)*AF**2*XWC**2*SH**2*SBWZ
0291               ICOMBI=3
0292               ALPHA=VECT/(VECT+AXIV)
0293             ELSEIF(KFSRCE.EQ.24.OR.KFSRCE.EQ.0) THEN
0294               ICOMBI=4
0295             ENDIF
0296 C...For chi -> chi q qbar, use V/A -> q qbar as first approximation.
0297           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.ITYPES.EQ.5) THEN
0298             ICLASS=2
0299           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.3.AND.(ITYPES.EQ.0.OR.
0300      &    ITYPES.EQ.1)) THEN
0301             ICLASS=3
0302  
0303 C...Scalar/pseudoscalar -> q + qbar; q -> q + S.
0304           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.1.AND.ITYPES.EQ.4) THEN
0305             ICLASS=4
0306             IF(KFSRCE.EQ.25.OR.KFSRCE.EQ.35.OR.KFSRCE.EQ.37) THEN
0307               ICOMBI=1
0308             ELSEIF(KFSRCE.EQ.36) THEN
0309               ICOMBI=2
0310             ENDIF
0311           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.4.AND.(ITYPES.EQ.0.OR.
0312      &    ITYPES.EQ.1)) THEN
0313             ICLASS=5
0314  
0315 C...V -> ~q + ~qbar; ~q -> ~q + V; S -> ~q + ~qbar; ~q -> ~q + S.
0316           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.2.AND.(ITYPES.EQ.0.OR.
0317      &    ITYPES.EQ.3)) THEN
0318             ICLASS=6
0319           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.3.AND.(ITYPES.EQ.0.OR.
0320      &    ITYPES.EQ.2)) THEN
0321             ICLASS=7
0322           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.2.AND.ITYPES.EQ.4) THEN
0323             ICLASS=8
0324           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.4.AND.(ITYPES.EQ.0.OR.
0325      &    ITYPES.EQ.2)) THEN
0326             ICLASS=9
0327  
0328 C...chi -> q + ~qbar; ~q -> q + chi; q -> ~q + chi.
0329           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.2.AND.(ITYPES.EQ.0.OR.
0330      &    ITYPES.EQ.5)) THEN
0331             ICLASS=10
0332           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.5.AND.(ITYPES.EQ.0.OR.
0333      &    ITYPES.EQ.2)) THEN
0334             ICLASS=11
0335           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.5.AND.(ITYPES.EQ.0.OR.
0336      &    ITYPES.EQ.1)) THEN
0337             ICLASS=12
0338  
0339 C...~g -> q + ~qbar; ~q -> q + ~g; q -> ~q + ~g.
0340           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.2.AND.ITYPES.EQ.6) THEN
0341             ICLASS=13
0342           ELSEIF(ITYPMN.EQ.1.AND.ITYPMX.EQ.6.AND.(ITYPES.EQ.0.OR.
0343      &    ITYPES.EQ.2)) THEN
0344             ICLASS=14
0345           ELSEIF(ITYPMN.EQ.2.AND.ITYPMX.EQ.6.AND.(ITYPES.EQ.0.OR.
0346      &    ITYPES.EQ.1)) THEN
0347             ICLASS=15
0348  
0349 C...g -> ~g + ~g (eikonal approximation).
0350           ELSEIF(ITYPMN.EQ.6.AND.ITYPMX.EQ.6.AND.ITYPES.EQ.0) THEN
0351             ICLASS=16
0352           ENDIF
0353           M3JC=5*ICLASS+ICOMBI
0354         ENDIF
0355       ENDIF
0356  
0357 C...Find if interference with initial state partons.
0358       MIIS=0
0359       IF(MSTJ(50).GE.1.AND.MSTJ(50).LE.3.AND.NPA.EQ.2.AND.KFSRCE.EQ.0
0360      &.AND.MPSPD.EQ.0) MIIS=MSTJ(50)
0361       IF(MSTJ(50).GE.4.AND.MSTJ(50).LE.6.AND.NPA.EQ.2.AND.MPSPD.EQ.0)
0362      &MIIS=MSTJ(50)-3
0363       IF(MIIS.NE.0) THEN
0364         DO 190 I=1,2
0365           KCII(I)=0
0366           KCA=PYCOMP(KFLA(I))
0367           IF(KCA.NE.0) KCII(I)=KCHG(KCA,2)*ISIGN(1,K(IPA(I),2))
0368           NIIS(I)=0
0369           IF(KCII(I).NE.0) THEN
0370             DO 180 J=1,2
0371               ICSI=MOD(K(IPA(I),3+J)/MSTU(5),MSTU(5))
0372               IF(ICSI.GT.0.AND.ICSI.NE.IPA(1).AND.ICSI.NE.IPA(2).AND.
0373      &        (KCII(I).EQ.(-1)**(J+1).OR.KCII(I).EQ.2)) THEN
0374                 NIIS(I)=NIIS(I)+1
0375                 IIIS(I,NIIS(I))=ICSI
0376               ENDIF
0377   180       CONTINUE
0378           ENDIF
0379   190   CONTINUE
0380         IF(NIIS(1)+NIIS(2).EQ.0) MIIS=0
0381       ENDIF
0382  
0383 C...Boost interfering initial partons to rest frame
0384 C...and reconstruct their polar and azimuthal angles.
0385       IF(MIIS.NE.0) THEN
0386         DO 210 I=1,2
0387           DO 200 J=1,5
0388             K(N+I,J)=K(IPA(I),J)
0389             P(N+I,J)=P(IPA(I),J)
0390             V(N+I,J)=0D0
0391   200     CONTINUE
0392   210   CONTINUE
0393         DO 230 I=3,2+NIIS(1)
0394           DO 220 J=1,5
0395             K(N+I,J)=K(IIIS(1,I-2),J)
0396             P(N+I,J)=P(IIIS(1,I-2),J)
0397             V(N+I,J)=0D0
0398   220     CONTINUE
0399   230   CONTINUE
0400         DO 250 I=3+NIIS(1),2+NIIS(1)+NIIS(2)
0401           DO 240 J=1,5
0402             K(N+I,J)=K(IIIS(2,I-2-NIIS(1)),J)
0403             P(N+I,J)=P(IIIS(2,I-2-NIIS(1)),J)
0404             V(N+I,J)=0D0
0405   240     CONTINUE
0406   250   CONTINUE
0407         CALL PYROBO(N+1,N+2+NIIS(1)+NIIS(2),0D0,0D0,-PS(1)/PS(4),
0408      &  -PS(2)/PS(4),-PS(3)/PS(4))
0409         PHI=PYANGL(P(N+1,1),P(N+1,2))
0410         CALL PYROBO(N+1,N+2+NIIS(1)+NIIS(2),0D0,-PHI,0D0,0D0,0D0)
0411         THE=PYANGL(P(N+1,3),P(N+1,1))
0412         CALL PYROBO(N+1,N+2+NIIS(1)+NIIS(2),-THE,0D0,0D0,0D0,0D0)
0413         DO 260 I=3,2+NIIS(1)
0414           THEIIS(1,I-2)=PYANGL(P(N+I,3),SQRT(P(N+I,1)**2+P(N+I,2)**2))
0415           PHIIIS(1,I-2)=PYANGL(P(N+I,1),P(N+I,2))
0416   260   CONTINUE
0417         DO 270 I=3+NIIS(1),2+NIIS(1)+NIIS(2)
0418           THEIIS(2,I-2-NIIS(1))=PARU(1)-PYANGL(P(N+I,3),
0419      &    SQRT(P(N+I,1)**2+P(N+I,2)**2))
0420           PHIIIS(2,I-2-NIIS(1))=PYANGL(P(N+I,1),P(N+I,2))
0421   270   CONTINUE
0422       ENDIF
0423  
0424 C...Boost 3 or more partons to their rest frame.
0425       IF(NPA.GE.3) CALL PYROBO(IPA(1),IPA(NPA),0D0,0D0,-PS(1)/PS(4),
0426      &-PS(2)/PS(4),-PS(3)/PS(4))
0427  
0428 C...Define imagined single initiator of shower for parton system.
0429       NS=N
0430       IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
0431         CALL PYERRM(11,'(PYSHOW:) no more memory left in PYJETS')
0432         IF(MSTU(21).GE.1) RETURN
0433       ENDIF
0434   280 N=NS
0435       IF(NPA.GE.2) THEN
0436         K(N+1,1)=11
0437         K(N+1,2)=21
0438         K(N+1,3)=0
0439         K(N+1,4)=0
0440         K(N+1,5)=0
0441         P(N+1,1)=0D0
0442         P(N+1,2)=0D0
0443         P(N+1,3)=0D0
0444         P(N+1,4)=PS(5)
0445         P(N+1,5)=PS(5)
0446         V(N+1,5)=PS(5)**2
0447         N=N+1
0448         IREF(1)=21
0449       ENDIF
0450  
0451 C...Loop over partons that may branch.
0452       NEP=NPA
0453       IM=NS
0454       IF(NPA.EQ.1) IM=NS-1
0455   290 IM=IM+1
0456       IF(N.GT.NS) THEN
0457         IF(IM.GT.N) GOTO 600
0458         KFLM=IABS(K(IM,2))
0459         IR=IREF(IM-NS)
0460         IF(KSH(IR).EQ.0) GOTO 290
0461         IF(P(IM,5).LT.PMTH(2,IR)) GOTO 290
0462         IGM=K(IM,3)
0463       ELSE
0464         IGM=-1
0465       ENDIF
0466       IF(N+NEP.GT.MSTU(4)-MSTU(32)-10) THEN
0467         CALL PYERRM(11,'(PYSHOW:) no more memory left in PYJETS')
0468         IF(MSTU(21).GE.1) RETURN
0469       ENDIF
0470  
0471 C...Position of aunt (sister to branching parton).
0472 C...Origin and flavour of daughters.
0473       IAU=0
0474       IF(IGM.GT.0) THEN
0475         IF(K(IM-1,3).EQ.IGM) IAU=IM-1
0476         IF(N.GE.IM+1.AND.K(IM+1,3).EQ.IGM) IAU=IM+1
0477       ENDIF
0478       IF(IGM.GE.0) THEN
0479         K(IM,4)=N+1
0480         DO 300 I=1,NEP
0481           K(N+I,3)=IM
0482   300   CONTINUE
0483       ELSE
0484         K(N+1,3)=IPA(1)
0485       ENDIF
0486       IF(IGM.LE.0) THEN
0487         DO 310 I=1,NEP
0488           K(N+I,2)=K(IPA(I),2)
0489   310   CONTINUE
0490       ELSEIF(KFLM.NE.21) THEN
0491         K(N+1,2)=K(IM,2)
0492         K(N+2,2)=K(IM,5)
0493         IREF(N+1-NS)=IREF(IM-NS)
0494         IREF(N+2-NS)=IABS(K(N+2,2))
0495       ELSEIF(K(IM,5).EQ.21) THEN
0496         K(N+1,2)=21
0497         K(N+2,2)=21
0498         IREF(N+1-NS)=21
0499         IREF(N+2-NS)=21
0500       ELSE
0501         K(N+1,2)=K(IM,5)
0502         K(N+2,2)=-K(IM,5)
0503         IREF(N+1-NS)=IABS(K(N+1,2))
0504         IREF(N+2-NS)=IABS(K(N+2,2))
0505       ENDIF
0506  
0507 C...Reset flags on daughters and tries made.
0508       DO 320 IP=1,NEP
0509         K(N+IP,1)=3
0510         K(N+IP,4)=0
0511         K(N+IP,5)=0
0512         KFLD(IP)=IABS(K(N+IP,2))
0513         IF(KCHG(PYCOMP(KFLD(IP)),2).EQ.0) K(N+IP,1)=1
0514         ITRY(IP)=0
0515         ISL(IP)=0
0516         ISI(IP)=0
0517         IF(KSH(IREF(N+IP-NS)).EQ.1) ISI(IP)=1
0518   320 CONTINUE
0519       ISLM=0
0520  
0521 C...Maximum virtuality of daughters.
0522       IF(IGM.LE.0) THEN
0523         DO 330 I=1,NPA
0524           IF(NPA.GE.3) P(N+I,4)=P(IPA(I),4)
0525           P(N+I,5)=MIN(QMAX,PS(5))
0526           IR=IREF(N+I-NS)
0527           IF(IP2.LE.-8) P(N+I,5)=MAX(P(N+I,5),2D0*PMTH(3,IR))
0528           IF(ISI(I).EQ.0) P(N+I,5)=P(IPA(I),5)
0529   330   CONTINUE
0530       ELSE
0531         IF(MSTJ(43).LE.2) PEM=V(IM,2)
0532         IF(MSTJ(43).GE.3) PEM=P(IM,4)
0533         P(N+1,5)=MIN(P(IM,5),V(IM,1)*PEM)
0534         P(N+2,5)=MIN(P(IM,5),(1D0-V(IM,1))*PEM)
0535         IF(K(N+2,2).EQ.22) P(N+2,5)=PMTH(1,22)
0536       ENDIF
0537       DO 340 I=1,NEP
0538         PMSD(I)=P(N+I,5)
0539         IF(ISI(I).EQ.1) THEN
0540           IR=IREF(N+I-NS)
0541           IF(P(N+I,5).LE.PMTH(3,IR)) P(N+I,5)=PMTH(1,IR)
0542         ENDIF
0543         V(N+I,5)=P(N+I,5)**2
0544   340 CONTINUE
0545  
0546 C...Choose one of the daughters for evolution.
0547   350 INUM=0
0548       IF(NEP.EQ.1) INUM=1
0549       DO 360 I=1,NEP
0550         IF(INUM.EQ.0.AND.ISL(I).EQ.1) INUM=I
0551   360 CONTINUE
0552       DO 370 I=1,NEP
0553         IF(INUM.EQ.0.AND.ITRY(I).EQ.0.AND.ISI(I).EQ.1) THEN
0554           IR=IREF(N+I-NS)
0555           IF(P(N+I,5).GE.PMTH(2,IR)) INUM=I
0556         ENDIF
0557   370 CONTINUE
0558       IF(INUM.EQ.0) THEN
0559         RMAX=0D0
0560         DO 380 I=1,NEP
0561           IF(ISI(I).EQ.1.AND.PMSD(I).GE.PMQT2E) THEN
0562             RPM=P(N+I,5)/PMSD(I)
0563             IR=IREF(N+I-NS)
0564             IF(RPM.GT.RMAX.AND.P(N+I,5).GE.PMTH(2,IR)) THEN
0565               RMAX=RPM
0566               INUM=I
0567             ENDIF
0568           ENDIF
0569   380   CONTINUE
0570       ENDIF
0571  
0572 C...Cancel choice of predetermined daughter already treated.
0573       INUM=MAX(1,INUM)
0574       INUMT=INUM
0575       IF(MPSPD.EQ.1.AND.IGM.EQ.0.AND.ITRY(INUMT).GE.1) THEN
0576         IF(K(IP1-1+INUM,4).GT.0) INUM=3-INUM
0577       ELSEIF(MPSPD.EQ.1.AND.IM.EQ.NS+2.AND.ITRY(INUMT).GE.1) THEN
0578         IF(KFLD(INUMT).NE.21.AND.K(IP1+2,4).GT.0) INUM=3-INUM
0579         IF(KFLD(INUMT).EQ.21.AND.K(IP1+3,4).GT.0) INUM=3-INUM
0580       ENDIF
0581  
0582 C...Store information on choice of evolving daughter.
0583       IEP(1)=N+INUM
0584       DO 390 I=2,NEP
0585         IEP(I)=IEP(I-1)+1
0586         IF(IEP(I).GT.N+NEP) IEP(I)=N+1
0587   390 CONTINUE
0588       DO 400 I=1,NEP
0589         KFL(I)=IABS(K(IEP(I),2))
0590   400 CONTINUE
0591       ITRY(INUM)=ITRY(INUM)+1
0592       IF(ITRY(INUM).GT.200) THEN
0593         CALL PYERRM(14,'(PYSHOW:) caught in infinite loop')
0594         IF(MSTU(21).GE.1) RETURN
0595       ENDIF
0596       Z=0.5D0
0597       IR=IREF(IEP(1)-NS)
0598       IF(KSH(IR).EQ.0) GOTO 450
0599       IF(P(IEP(1),5).LT.PMTH(2,IR)) GOTO 450
0600  
0601 C...Check if evolution already predetermined for daughter.
0602       IPSPD=0
0603       IF(MPSPD.EQ.1.AND.IGM.EQ.0) THEN
0604         IF(K(IP1-1+INUM,4).GT.0) IPSPD=IP1-1+INUM
0605       ELSEIF(MPSPD.EQ.1.AND.IM.EQ.NS+2) THEN
0606         IF(KFL(1).NE.21.AND.K(IP1+2,4).GT.0) IPSPD=IP1+2
0607         IF(KFL(1).EQ.21.AND.K(IP1+3,4).GT.0) IPSPD=IP1+3
0608       ENDIF
0609       IF(INUM.EQ.1.OR.INUM.EQ.2) THEN
0610         ISSET(INUM)=0
0611         IF(IPSPD.NE.0) ISSET(INUM)=1
0612       ENDIF
0613  
0614 C...Select side for interference with initial state partons.
0615       IF(MIIS.GE.1.AND.IEP(1).LE.NS+3) THEN
0616         III=IEP(1)-NS-1
0617         ISII(III)=0
0618         IF(IABS(KCII(III)).EQ.1.AND.NIIS(III).EQ.1) THEN
0619           ISII(III)=1
0620         ELSEIF(KCII(III).EQ.2.AND.NIIS(III).EQ.1) THEN
0621           IF(PYR(0).GT.0.5D0) ISII(III)=1
0622         ELSEIF(KCII(III).EQ.2.AND.NIIS(III).EQ.2) THEN
0623           ISII(III)=1
0624           IF(PYR(0).GT.0.5D0) ISII(III)=2
0625         ENDIF
0626       ENDIF
0627  
0628 C...Calculate allowed z range.
0629       IF(NEP.EQ.1) THEN
0630         PMED=PS(4)
0631       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
0632         PMED=P(IM,5)
0633       ELSE
0634         IF(INUM.EQ.1) PMED=V(IM,1)*PEM
0635         IF(INUM.EQ.2) PMED=(1D0-V(IM,1))*PEM
0636       ENDIF
0637       IF(MOD(MSTJ(43),2).EQ.1) THEN
0638         ZC=PMTH(2,21)/PMED
0639         ZCE=PMTH(2,22)/PMED
0640         IF(ISCOL(IR).EQ.0) ZCE=0.5D0*PARJ(90)/PMED
0641       ELSE
0642         ZC=0.5D0*(1D0-SQRT(MAX(0D0,1D0-(2D0*PMTH(2,21)/PMED)**2)))
0643         IF(ZC.LT.1D-6) ZC=(PMTH(2,21)/PMED)**2
0644         PMTMPE=PMTH(2,22)
0645         IF(ISCOL(IR).EQ.0) PMTMPE=0.5D0*PARJ(90)
0646         ZCE=0.5D0*(1D0-SQRT(MAX(0D0,1D0-(2D0*PMTMPE/PMED)**2)))
0647         IF(ZCE.LT.1D-6) ZCE=(PMTMPE/PMED)**2
0648       ENDIF
0649       ZC=MIN(ZC,0.491D0)
0650       ZCE=MIN(ZCE,0.49991D0)
0651       IF(((MSTJ(41).EQ.1.AND.ZC.GT.0.49D0).OR.(MSTJ(41).GE.2.AND.
0652      &MIN(ZC,ZCE).GT.0.4999D0)).AND.IPSPD.EQ.0) THEN
0653         P(IEP(1),5)=PMTH(1,IR)
0654         V(IEP(1),5)=P(IEP(1),5)**2
0655         GOTO 450
0656       ENDIF
0657  
0658 C...Integral of Altarelli-Parisi z kernel for QCD.
0659 C...(Includes squark and gluino; with factor N_C/C_F extra for latter).
0660       IF(MSTJ(49).EQ.0.AND.KFL(1).EQ.21) THEN
0661         FBR=6D0*LOG((1D0-ZC)/ZC)+MSTJ(45)*0.5D0
0662 C...QUARKONIA+++
0663 C...Evolution of QQ~[3S18] state if MSTP(148)=1.
0664       ELSEIF(MSTJ(49).EQ.0.AND.MSTP(149).GE.0.AND.
0665      &       (KFL(1).EQ.9900443.OR.KFL(1).EQ.9900553)) THEN
0666         FBR=6D0*LOG((1D0-ZC)/ZC)
0667 C...QUARKONIA---
0668       ELSEIF(MSTJ(49).EQ.0) THEN
0669         FBR=(8D0/3D0)*LOG((1D0-ZC)/ZC)
0670         IF(IGLUI.EQ.1.AND.IR.GE.31) FBR=FBR*(9D0/4D0)
0671  
0672 C...Integral of Altarelli-Parisi z kernel for scalar gluon.
0673       ELSEIF(MSTJ(49).EQ.1.AND.KFL(1).EQ.21) THEN
0674         FBR=(PARJ(87)+MSTJ(45)*PARJ(88))*(1D0-2D0*ZC)
0675       ELSEIF(MSTJ(49).EQ.1) THEN
0676         FBR=(1D0-2D0*ZC)/3D0
0677         IF(IGM.EQ.0.AND.M3JC.GE.1) FBR=4D0*FBR
0678  
0679 C...Integral of Altarelli-Parisi z kernel for Abelian vector gluon.
0680       ELSEIF(KFL(1).EQ.21) THEN
0681         FBR=6D0*MSTJ(45)*(0.5D0-ZC)
0682       ELSE
0683         FBR=2D0*LOG((1D0-ZC)/ZC)
0684       ENDIF
0685  
0686 C...Reset QCD probability for colourless.
0687       IF(ISCOL(IR).EQ.0) FBR=0D0
0688  
0689 C...Integral of Altarelli-Parisi kernel for photon emission.
0690       FBRE=0D0
0691       IF(MSTJ(41).GE.2.AND.ISCHG(IR).EQ.1) THEN
0692         IF(KFL(1).LE.18) THEN
0693           FBRE=(KCHG(KFL(1),1)/3D0)**2*2D0*LOG((1D0-ZCE)/ZCE)
0694         ENDIF
0695         IF(MSTJ(41).EQ.10) FBRE=PARJ(84)*FBRE
0696       ENDIF
0697  
0698 C...Inner veto algorithm starts. Find maximum mass for evolution.
0699   410 PMS=V(IEP(1),5)
0700       IF(IGM.GE.0) THEN
0701         PM2=0D0
0702         DO 420 I=2,NEP
0703           PM=P(IEP(I),5)
0704           IRI=IREF(IEP(I)-NS)
0705           IF(KSH(IRI).EQ.1) PM=PMTH(2,IRI)
0706           PM2=PM2+PM
0707   420   CONTINUE
0708         PMS=MIN(PMS,(P(IM,5)-PM2)**2)
0709       ENDIF
0710  
0711 C...Select mass for daughter in QCD evolution.
0712       B0=27D0/6D0
0713       DO 430 IFF=4,MSTJ(45)
0714         IF(PMS.GT.4D0*PMTH(2,IFF)**2) B0=(33D0-2D0*IFF)/6D0
0715   430 CONTINUE
0716 C...Shift m^2 for evolution in Q^2 = m^2 - m(onshell)^2.
0717       PMSC=MAX(0.5D0*PARJ(82),PMS-PMTH(1,IR)**2)
0718 C...Already predetermined choice.
0719       IF(IPSPD.NE.0) THEN
0720         PMSQCD=P(IPSPD,5)**2
0721       ELSEIF(FBR.LT.1D-3) THEN
0722         PMSQCD=0D0
0723       ELSEIF(MSTJ(44).LE.0) THEN
0724         PMSQCD=PMSC*EXP(MAX(-50D0,LOG(PYR(0))*PARU(2)/(PARU(111)*FBR)))
0725       ELSEIF(MSTJ(44).EQ.1) THEN
0726         PMSQCD=4D0*ALAMS*(0.25D0*PMSC/ALAMS)**(PYR(0)**(B0/FBR))
0727       ELSE
0728         PMSQCD=PMSC*EXP(MAX(-50D0,ALFM*B0*LOG(PYR(0))/FBR))
0729       ENDIF
0730 C...Shift back m^2 from evolution in Q^2 = m^2 - m(onshell)^2.
0731       IF(IPSPD.EQ.0) PMSQCD=PMSQCD+PMTH(1,IR)**2
0732       IF(ZC.GT.0.49D0.OR.PMSQCD.LE.PMTH(4,IR)**2) PMSQCD=PMTH(2,IR)**2
0733       V(IEP(1),5)=PMSQCD
0734       MCE=1
0735  
0736 C...Select mass for daughter in QED evolution.
0737       IF(MSTJ(41).GE.2.AND.ISCHG(IR).EQ.1.AND.IPSPD.EQ.0) THEN
0738 C...Shift m^2 for evolution in Q^2 = m^2 - m(onshell)^2.
0739         PMSE=MAX(0.5D0*PARJ(83),PMS-PMTH(1,IR)**2)
0740         IF(FBRE.LT.1D-3) THEN
0741           PMSQED=0D0
0742         ELSE
0743           PMSQED=PMSE*EXP(MAX(-50D0,LOG(PYR(0))*PARU(2)/
0744      &    (PARU(101)*FBRE)))
0745         ENDIF
0746 C...Shift back m^2 from evolution in Q^2 = m^2 - m(onshell)^2.
0747         PMSQED=PMSQED+PMTH(1,IR)**2
0748         IF(ZCE.GT.0.4999D0.OR.PMSQED.LE.PMTH(5,IR)**2) PMSQED=
0749      &  PMTH(2,IR)**2
0750         IF(PMSQED.GT.PMSQCD) THEN
0751           V(IEP(1),5)=PMSQED
0752           MCE=2
0753         ENDIF
0754       ENDIF
0755  
0756 C...Check whether daughter mass below cutoff.
0757       P(IEP(1),5)=SQRT(V(IEP(1),5))
0758       IF(P(IEP(1),5).LE.PMTH(3,IR)) THEN
0759         P(IEP(1),5)=PMTH(1,IR)
0760         V(IEP(1),5)=P(IEP(1),5)**2
0761         GOTO 450
0762       ENDIF
0763  
0764 C...Already predetermined choice of z, and flavour in g -> qqbar.
0765       IF(IPSPD.NE.0) THEN
0766         IPSGD1=K(IPSPD,4)
0767         IPSGD2=K(IPSPD,5)
0768         PMSGD1=P(IPSGD1,5)**2
0769         PMSGD2=P(IPSGD2,5)**2
0770         ALAMPS=SQRT(MAX(1D-10,(PMSQCD-PMSGD1-PMSGD2)**2-
0771      &  4D0*PMSGD1*PMSGD2))
0772         Z=0.5D0*(PMSQCD*(2D0*P(IPSGD1,4)/P(IPSPD,4)-1D0)+ALAMPS-
0773      &  PMSGD1+PMSGD2)/ALAMPS
0774         Z=MAX(0.00001D0,MIN(0.99999D0,Z))
0775         IF(KFL(1).NE.21) THEN
0776           K(IEP(1),5)=21
0777         ELSE
0778           K(IEP(1),5)=IABS(K(IPSGD1,2))
0779         ENDIF
0780  
0781 C...Select z value of branching: q -> qgamma.
0782       ELSEIF(MCE.EQ.2) THEN
0783         Z=1D0-(1D0-ZCE)*(ZCE/(1D0-ZCE))**PYR(0)
0784         IF(1D0+Z**2.LT.2D0*PYR(0)) GOTO 410
0785         K(IEP(1),5)=22
0786  
0787 C...QUARKONIA+++
0788 C...Select z value of branching: QQ~[3S18] -> QQ~[3S18]g.
0789       ELSEIF(MSTJ(49).EQ.0.AND.
0790      &       (KFL(1).EQ.9900443.OR.KFL(1).EQ.9900553)) THEN
0791         Z=(1D0-ZC)*(ZC/(1D0-ZC))**PYR(0)
0792 C...Select always the harder 'gluon' if the switch MSTP(149)<=0.
0793         IF(MSTP(149).LE.0.OR.PYR(0).GT.0.5D0) Z=1D0-Z
0794         IF((1D0-Z*(1D0-Z))**2.LT.PYR(0)) GOTO 410
0795         K(IEP(1),5)=21
0796 C...QUARKONIA---
0797  
0798 C...Select z value of branching: q -> qg, g -> gg, g -> qqbar.
0799       ELSEIF(MSTJ(49).NE.1.AND.KFL(1).NE.21) THEN
0800         Z=1D0-(1D0-ZC)*(ZC/(1D0-ZC))**PYR(0)
0801 C...Only do z weighting when no ME correction afterwards.
0802         IF(M3JC.EQ.0.AND.1D0+Z**2.LT.2D0*PYR(0)) GOTO 410
0803         K(IEP(1),5)=21
0804       ELSEIF(MSTJ(49).EQ.0.AND.MSTJ(45)*0.5D0.LT.PYR(0)*FBR) THEN
0805         Z=(1D0-ZC)*(ZC/(1D0-ZC))**PYR(0)
0806         IF(PYR(0).GT.0.5D0) Z=1D0-Z
0807         IF((1D0-Z*(1D0-Z))**2.LT.PYR(0)) GOTO 410
0808         K(IEP(1),5)=21
0809       ELSEIF(MSTJ(49).NE.1) THEN
0810         Z=PYR(0)
0811         IF(Z**2+(1D0-Z)**2.LT.PYR(0)) GOTO 410
0812         KFLB=1+INT(MSTJ(45)*PYR(0))
0813         PMQ=4D0*PMTH(2,KFLB)**2/V(IEP(1),5)
0814         IF(PMQ.GE.1D0) GOTO 410
0815         IF(MSTJ(44).LE.2.OR.MSTJ(44).EQ.4) THEN
0816           IF(Z.LT.ZC.OR.Z.GT.1D0-ZC) GOTO 410
0817           PMQ0=4D0*PMTH(2,21)**2/V(IEP(1),5)
0818           IF(MOD(MSTJ(43),2).EQ.0.AND.(1D0+0.5D0*PMQ)*SQRT(1D0-PMQ)
0819      &    .LT.PYR(0)*(1D0+0.5D0*PMQ0)*SQRT(1D0-PMQ0)) GOTO 410
0820         ELSE
0821           IF((1D0+0.5D0*PMQ)*SQRT(1D0-PMQ).LT.PYR(0)) GOTO 410
0822         ENDIF
0823         K(IEP(1),5)=KFLB
0824  
0825 C...Ditto for scalar gluon model.
0826       ELSEIF(KFL(1).NE.21) THEN
0827         Z=1D0-SQRT(ZC**2+PYR(0)*(1D0-2D0*ZC))
0828         K(IEP(1),5)=21
0829       ELSEIF(PYR(0)*(PARJ(87)+MSTJ(45)*PARJ(88)).LE.PARJ(87)) THEN
0830         Z=ZC+(1D0-2D0*ZC)*PYR(0)
0831         K(IEP(1),5)=21
0832       ELSE
0833         Z=ZC+(1D0-2D0*ZC)*PYR(0)
0834         KFLB=1+INT(MSTJ(45)*PYR(0))
0835         PMQ=4D0*PMTH(2,KFLB)**2/V(IEP(1),5)
0836         IF(PMQ.GE.1D0) GOTO 410
0837         K(IEP(1),5)=KFLB
0838       ENDIF
0839  
0840 C...Correct to alpha_s(pT^2) (optionally m^2/4 for g -> q qbar).
0841       IF(MCE.EQ.1.AND.MSTJ(44).GE.2.AND.IPSPD.EQ.0) THEN
0842         IF(KFL(1).EQ.21.AND.K(IEP(1),5).LT.10.AND.
0843      &  (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
0844           IF(ALFM/LOG(V(IEP(1),5)*0.25D0/ALAMS).LT.PYR(0)) GOTO 410
0845         ELSE
0846           PT2APP=Z*(1D0-Z)*V(IEP(1),5)
0847           IF(MSTJ(44).GE.4) PT2APP=PT2APP*
0848      &    (1D0-PMTH(1,IR)**2/V(IEP(1),5))**2
0849           IF(PT2APP.LT.PT2MIN) GOTO 410
0850           IF(ALFM/LOG(PT2APP/ALAMS).LT.PYR(0)) GOTO 410
0851         ENDIF
0852       ENDIF
0853  
0854 C...Check if z consistent with chosen m.
0855       IF(KFL(1).EQ.21) THEN
0856         IRGD1=IABS(K(IEP(1),5))
0857         IRGD2=IRGD1
0858       ELSE
0859         IRGD1=IR
0860         IRGD2=IABS(K(IEP(1),5))
0861       ENDIF
0862       IF(NEP.EQ.1) THEN
0863         PED=PS(4)
0864       ELSEIF(NEP.GE.3) THEN
0865         PED=P(IEP(1),4)
0866       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
0867         PED=0.5D0*(V(IM,5)+V(IEP(1),5)-PM2**2)/P(IM,5)
0868       ELSE
0869         IF(IEP(1).EQ.N+1) PED=V(IM,1)*PEM
0870         IF(IEP(1).EQ.N+2) PED=(1D0-V(IM,1))*PEM
0871       ENDIF
0872       IF(MOD(MSTJ(43),2).EQ.1) THEN
0873         PMQTH3=0.5D0*PARJ(82)
0874         IF(IRGD2.EQ.22) PMQTH3=0.5D0*PARJ(83)
0875         IF(IRGD2.EQ.22.AND.ISCOL(IR).EQ.0) PMQTH3=0.5D0*PARJ(90)
0876         PMQ1=(PMTH(1,IRGD1)**2+PMQTH3**2)/V(IEP(1),5)
0877         PMQ2=(PMTH(1,IRGD2)**2+PMQTH3**2)/V(IEP(1),5)
0878         ZD=SQRT(MAX(0D0,(1D0-V(IEP(1),5)/PED**2)*((1D0-PMQ1-PMQ2)**2-
0879      &  4D0*PMQ1*PMQ2)))
0880         ZH=1D0+PMQ1-PMQ2
0881       ELSE
0882         ZD=SQRT(MAX(0D0,1D0-V(IEP(1),5)/PED**2))
0883         ZH=1D0
0884       ENDIF
0885       IF(KFL(1).EQ.21.AND.K(IEP(1),5).LT.10.AND.
0886      &(MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
0887       ELSEIF(IPSPD.NE.0) THEN
0888       ELSE
0889         ZL=0.5D0*(ZH-ZD)
0890         ZU=0.5D0*(ZH+ZD)
0891         IF(Z.LT.ZL.OR.Z.GT.ZU) GOTO 410
0892       ENDIF
0893       IF(KFL(1).EQ.21) V(IEP(1),3)=LOG(ZU*(1D0-ZL)/MAX(1D-20,ZL*
0894      &(1D0-ZU)))
0895       IF(KFL(1).NE.21) V(IEP(1),3)=LOG((1D0-ZL)/MAX(1D-10,1D0-ZU))
0896  
0897 C...Width suppression for q -> q + g.
0898       IF(MSTJ(40).NE.0.AND.KFL(1).NE.21.AND.IPSPD.EQ.0) THEN
0899         IF(IGM.EQ.0) THEN
0900           EGLU=0.5D0*PS(5)*(1D0-Z)*(1D0+V(IEP(1),5)/V(NS+1,5))
0901         ELSE
0902           EGLU=PMED*(1D0-Z)
0903         ENDIF
0904         CHI=PARJ(89)**2/(PARJ(89)**2+EGLU**2)
0905         IF(MSTJ(40).EQ.1) THEN
0906           IF(CHI.LT.PYR(0)) GOTO 410
0907         ELSEIF(MSTJ(40).EQ.2) THEN
0908           IF(1D0-CHI.LT.PYR(0)) GOTO 410
0909         ENDIF
0910       ENDIF
0911  
0912 C...Three-jet matrix element correction.
0913       IF(M3JC.GE.1) THEN
0914         WME=1D0
0915         WSHOW=1D0
0916  
0917 C...QED matrix elements: only for massless case so far.
0918         IF(MCE.EQ.2.AND.IGM.EQ.0) THEN
0919           X1=Z*(1D0+V(IEP(1),5)/V(NS+1,5))
0920           X2=1D0-V(IEP(1),5)/V(NS+1,5)
0921           X3=(1D0-X1)+(1D0-X2)
0922           KI1=K(IPA(INUM),2)
0923           KI2=K(IPA(3-INUM),2)
0924           QF1=KCHG(PYCOMP(KI1),1)*ISIGN(1,KI1)/3D0
0925           QF2=KCHG(PYCOMP(KI2),1)*ISIGN(1,KI2)/3D0
0926           WSHOW=QF1**2*(1D0-X1)/X3*(1D0+(X1/(2D0-X2))**2)+
0927      &    QF2**2*(1D0-X2)/X3*(1D0+(X2/(2D0-X1))**2)
0928           WME=(QF1*(1D0-X1)/X3-QF2*(1D0-X2)/X3)**2*(X1**2+X2**2)
0929         ELSEIF(MCE.EQ.2) THEN
0930  
0931 C...QCD matrix elements, including mass effects.
0932         ELSEIF(MSTJ(49).NE.1.AND.K(IEP(1),2).NE.21) THEN
0933           PS1ME=V(IEP(1),5)
0934           PM1ME=PMTH(1,IR)
0935           M3JCC=M3JC
0936           IF(IR.GE.31.AND.IGM.EQ.0) THEN
0937 C...QCD ME: original parton, first branching.
0938             PM2ME=PMTH(1,63-IR)
0939             ECMME=PS(5)
0940           ELSEIF(IR.GE.31) THEN
0941 C...QCD ME: original parton, subsequent branchings.
0942             PM2ME=PMTH(1,63-IR)
0943             PEDME=PEM*(V(IM,1)+(1D0-V(IM,1))*PS1ME/V(IM,5))
0944             ECMME=PEDME+SQRT(MAX(0D0,PEDME**2-PS1ME+PM2ME**2))
0945           ELSEIF(K(IM,2).EQ.21) THEN
0946 C...QCD ME: secondary partons, first branching.
0947             PM2ME=PM1ME
0948             ZMME=V(IM,1)
0949             IF(IEP(1).GT.IEP(2)) ZMME=1D0-ZMME
0950             PMLME=SQRT(MAX(0D0,(V(IM,5)-PS1ME-PM2ME**2)**2-
0951      &      4D0*PS1ME*PM2ME**2))
0952             PEDME=PEM*(0.5D0*(V(IM,5)-PMLME+PS1ME-PM2ME**2)+PMLME*ZMME)/
0953      &      V(IM,5)
0954             ECMME=PEDME+SQRT(MAX(0D0,PEDME**2-PS1ME+PM2ME**2))
0955             M3JCC=66
0956           ELSE
0957 C...QCD ME: secondary partons, subsequent branchings.
0958             PM2ME=PM1ME
0959             PEDME=PEM*(V(IM,1)+(1D0-V(IM,1))*PS1ME/V(IM,5))
0960             ECMME=PEDME+SQRT(MAX(0D0,PEDME**2-PS1ME+PM2ME**2))
0961             M3JCC=66
0962           ENDIF
0963 C...Construct ME variables.
0964           R1ME=PM1ME/ECMME
0965           R2ME=PM2ME/ECMME
0966           X1=(1D0+PS1ME/ECMME**2-R2ME**2)*(Z+(1D0-Z)*PM1ME**2/PS1ME)
0967           X2=1D0+R2ME**2-PS1ME/ECMME**2
0968 C...Call ME, with right order important for two inequivalent showerers.
0969           IF(IR.EQ.IORD+30) THEN
0970             WME=PYMAEL(M3JCC,X1,X2,R1ME,R2ME,ALPHA)
0971           ELSE
0972             WME=PYMAEL(M3JCC,X2,X1,R2ME,R1ME,ALPHA)
0973           ENDIF
0974 C...Split up total ME when two radiating partons.
0975           ISPRAD=1
0976           IF((M3JCC.GE.16.AND.M3JCC.LE.19).OR.
0977      &    (M3JCC.GE.26.AND.M3JCC.LE.29).OR.
0978      &    (M3JCC.GE.36.AND.M3JCC.LE.39).OR.
0979      &    (M3JCC.GE.46.AND.M3JCC.LE.49).OR.
0980      &    (M3JCC.GE.56.AND.M3JCC.LE.64)) ISPRAD=0
0981           IF(ISPRAD.EQ.1) WME=WME*MAX(1D-10,1D0+R1ME**2-R2ME**2-X1)/
0982      &    MAX(1D-10,2D0-X1-X2)
0983 C...Evaluate shower rate to be compared with.
0984           WSHOW=2D0/(MAX(1D-10,2D0-X1-X2)*
0985      &    MAX(1D-10,1D0+R2ME**2-R1ME**2-X2))
0986           IF(IGLUI.EQ.1.AND.IR.GE.31) WSHOW=(9D0/4D0)*WSHOW
0987         ELSEIF(MSTJ(49).NE.1) THEN
0988  
0989 C...Toy model scalar theory matrix elements; no mass effects.
0990         ELSE
0991           X1=Z*(1D0+V(IEP(1),5)/V(NS+1,5))
0992           X2=1D0-V(IEP(1),5)/V(NS+1,5)
0993           X3=(1D0-X1)+(1D0-X2)
0994           WSHOW=4D0*X3*((1D0-X1)/(2D0-X2)**2+(1D0-X2)/(2D0-X1)**2)
0995           WME=X3**2
0996           IF(MSTJ(102).GE.2) WME=X3**2-2D0*(1D0+X3)*(1D0-X1)*(1D0-X2)*
0997      &    PARJ(171)
0998         ENDIF
0999  
1000         IF(WME.LT.PYR(0)*WSHOW) GOTO 410
1001       ENDIF
1002  
1003 C...Impose angular ordering by rejection of nonordered emission.
1004       IF(MCE.EQ.1.AND.IGM.GT.0.AND.MSTJ(42).GE.2.AND.IPSPD.EQ.0) THEN
1005         PEMAO=V(IM,1)*P(IM,4)
1006         IF(IEP(1).EQ.N+2) PEMAO=(1D0-V(IM,1))*P(IM,4)
1007         IF(IR.GE.31.AND.MSTJ(42).GE.5) THEN
1008           MAOD=0
1009         ELSEIF(KFL(1).EQ.21.AND.K(IEP(1),5).LE.10.AND.(MSTJ(42).EQ.4
1010      &  .OR.MSTJ(42).EQ.7)) THEN
1011           MAOD=0
1012         ELSEIF(KFL(1).EQ.21.AND.K(IEP(1),5).LE.10.AND.(MSTJ(42).EQ.3
1013      &  .OR.MSTJ(42).EQ.6)) THEN
1014           MAOD=1
1015           PMDAO=PMTH(2,K(IEP(1),5))
1016           THE2ID=Z*(1D0-Z)*PEMAO**2/(V(IEP(1),5)-4D0*PMDAO**2)
1017         ELSE
1018           MAOD=1
1019           THE2ID=Z*(1D0-Z)*PEMAO**2/V(IEP(1),5)
1020           IF(MSTJ(42).GE.3.AND.MSTJ(42).NE.5) THE2ID=THE2ID*
1021      &    (1D0+PMTH(1,IR)**2*(1D0-Z)/(V(IEP(1),5)*Z))**2
1022         ENDIF
1023         MAOM=1
1024         IAOM=IM
1025   440   IF(K(IAOM,5).EQ.22) THEN
1026           IAOM=K(IAOM,3)
1027           IF(K(IAOM,3).LE.NS) MAOM=0
1028           IF(MAOM.EQ.1) GOTO 440
1029         ENDIF
1030         IF(MAOM.EQ.1.AND.MAOD.EQ.1) THEN
1031           THE2IM=V(IAOM,1)*(1D0-V(IAOM,1))*P(IAOM,4)**2/V(IAOM,5)
1032           IF(THE2ID.LT.THE2IM) GOTO 410
1033         ENDIF
1034       ENDIF
1035  
1036 C...Impose user-defined maximum angle at first branching.
1037       IF(MSTJ(48).EQ.1.AND.IPSPD.EQ.0) THEN
1038         IF(NEP.EQ.1.AND.IM.EQ.NS) THEN
1039           THE2ID=Z*(1D0-Z)*PS(4)**2/V(IEP(1),5)
1040           IF(PARJ(85)**2*THE2ID.LT.1D0) GOTO 410
1041         ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+2) THEN
1042           THE2ID=Z*(1D0-Z)*(0.5D0*P(IM,4))**2/V(IEP(1),5)
1043           IF(PARJ(85)**2*THE2ID.LT.1D0) GOTO 410
1044         ELSEIF(NEP.EQ.2.AND.IEP(1).EQ.NS+3) THEN
1045           THE2ID=Z*(1D0-Z)*(0.5D0*P(IM,4))**2/V(IEP(1),5)
1046           IF(PARJ(86)**2*THE2ID.LT.1D0) GOTO 410
1047         ENDIF
1048       ENDIF
1049  
1050 C...Impose angular constraint in first branching from interference
1051 C...with initial state partons.
1052       IF(MIIS.GE.2.AND.IEP(1).LE.NS+3) THEN
1053         THE2D=MAX((1D0-Z)/Z,Z/(1D0-Z))*V(IEP(1),5)/(0.5D0*P(IM,4))**2
1054         IF(IEP(1).EQ.NS+2.AND.ISII(1).GE.1) THEN
1055           IF(THE2D.GT.THEIIS(1,ISII(1))**2) GOTO 410
1056         ELSEIF(IEP(1).EQ.NS+3.AND.ISII(2).GE.1) THEN
1057           IF(THE2D.GT.THEIIS(2,ISII(2))**2) GOTO 410
1058         ENDIF
1059       ENDIF
1060  
1061 C...End of inner veto algorithm. Check if only one leg evolved so far.
1062   450 V(IEP(1),1)=Z
1063       ISL(1)=0
1064       ISL(2)=0
1065       IF(NEP.EQ.1) GOTO 490
1066       IF(NEP.EQ.2.AND.P(IEP(1),5)+P(IEP(2),5).GE.P(IM,5)) GOTO 350
1067       DO 460 I=1,NEP
1068         IR=IREF(N+I-NS)
1069         IF(ITRY(I).EQ.0.AND.KSH(IR).EQ.1) THEN
1070           IF(P(N+I,5).GE.PMTH(2,IR)) GOTO 350
1071         ENDIF
1072   460 CONTINUE
1073  
1074 C...Check if chosen multiplet m1,m2,z1,z2 is physical.
1075       IF(NEP.GE.3) THEN
1076         PMSUM=0D0
1077         DO 470 I=1,NEP
1078           PMSUM=PMSUM+P(N+I,5)
1079   470   CONTINUE
1080         IF(PMSUM.GE.PS(5)) GOTO 350
1081       ELSEIF(IGM.EQ.0.OR.MSTJ(43).LE.2.OR.MOD(MSTJ(43),2).EQ.0) THEN
1082         DO 480 I1=N+1,N+2
1083           IRDA=IREF(I1-NS)
1084           IF(KSH(IRDA).EQ.0) GOTO 480
1085           IF(P(I1,5).LT.PMTH(2,IRDA)) GOTO 480
1086           IF(IRDA.EQ.21) THEN
1087             IRGD1=IABS(K(I1,5))
1088             IRGD2=IRGD1
1089           ELSE
1090             IRGD1=IRDA
1091             IRGD2=IABS(K(I1,5))
1092           ENDIF
1093           I2=2*N+3-I1
1094           IF(IGM.EQ.0.OR.MSTJ(43).LE.2) THEN
1095             PED=0.5D0*(V(IM,5)+V(I1,5)-V(I2,5))/P(IM,5)
1096           ELSE
1097             IF(I1.EQ.N+1) ZM=V(IM,1)
1098             IF(I1.EQ.N+2) ZM=1D0-V(IM,1)
1099             PML=SQRT((V(IM,5)-V(N+1,5)-V(N+2,5))**2-
1100      &      4D0*V(N+1,5)*V(N+2,5))
1101             PED=PEM*(0.5D0*(V(IM,5)-PML+V(I1,5)-V(I2,5))+PML*ZM)/
1102      &      V(IM,5)
1103           ENDIF
1104           IF(MOD(MSTJ(43),2).EQ.1) THEN
1105             PMQTH3=0.5D0*PARJ(82)
1106             IF(IRGD2.EQ.22) PMQTH3=0.5D0*PARJ(83)
1107             IF(IRGD2.EQ.22.AND.ISCOL(IRDA).EQ.0) PMQTH3=0.5D0*PARJ(90)
1108             PMQ1=(PMTH(1,IRGD1)**2+PMQTH3**2)/V(I1,5)
1109             PMQ2=(PMTH(1,IRGD2)**2+PMQTH3**2)/V(I1,5)
1110             ZD=SQRT(MAX(0D0,(1D0-V(I1,5)/PED**2)*((1D0-PMQ1-PMQ2)**2-
1111      &      4D0*PMQ1*PMQ2)))
1112             ZH=1D0+PMQ1-PMQ2
1113           ELSE
1114             ZD=SQRT(MAX(0D0,1D0-V(I1,5)/PED**2))
1115             ZH=1D0
1116           ENDIF
1117           IF(IRDA.EQ.21.AND.IRGD1.LT.10.AND.
1118      &    (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
1119           ELSE
1120             ZL=0.5D0*(ZH-ZD)
1121             ZU=0.5D0*(ZH+ZD)
1122             IF(I1.EQ.N+1.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU).AND.
1123      &      ISSET(1).EQ.0) THEN
1124               ISL(1)=1
1125             ELSEIF(I1.EQ.N+2.AND.(V(I1,1).LT.ZL.OR.V(I1,1).GT.ZU).AND.
1126      &      ISSET(2).EQ.0) THEN
1127               ISL(2)=1
1128             ENDIF
1129           ENDIF
1130           IF(IRDA.EQ.21) V(I1,4)=LOG(ZU*(1D0-ZL)/MAX(1D-20,
1131      &    ZL*(1D0-ZU)))
1132           IF(IRDA.NE.21) V(I1,4)=LOG((1D0-ZL)/MAX(1D-10,1D0-ZU))
1133   480   CONTINUE
1134         IF(ISL(1).EQ.1.AND.ISL(2).EQ.1.AND.ISLM.NE.0) THEN
1135           ISL(3-ISLM)=0
1136           ISLM=3-ISLM
1137         ELSEIF(ISL(1).EQ.1.AND.ISL(2).EQ.1) THEN
1138           ZDR1=MAX(0D0,V(N+1,3)/MAX(1D-6,V(N+1,4))-1D0)
1139           ZDR2=MAX(0D0,V(N+2,3)/MAX(1D-6,V(N+2,4))-1D0)
1140           IF(ZDR2.GT.PYR(0)*(ZDR1+ZDR2)) ISL(1)=0
1141           IF(ISL(1).EQ.1) ISL(2)=0
1142           IF(ISL(1).EQ.0) ISLM=1
1143           IF(ISL(2).EQ.0) ISLM=2
1144         ENDIF
1145         IF(ISL(1).EQ.1.OR.ISL(2).EQ.1) GOTO 350
1146       ENDIF
1147       IRD1=IREF(N+1-NS)
1148       IRD2=IREF(N+2-NS)
1149       IF(IGM.GT.0) THEN
1150         IF(MOD(MSTJ(43),2).EQ.1.AND.(P(N+1,5).GE.
1151      &  PMTH(2,IRD1).OR.P(N+2,5).GE.PMTH(2,IRD2))) THEN
1152           PMQ1=V(N+1,5)/V(IM,5)
1153           PMQ2=V(N+2,5)/V(IM,5)
1154           ZD=SQRT(MAX(0D0,(1D0-V(IM,5)/PEM**2)*((1D0-PMQ1-PMQ2)**2-
1155      &    4D0*PMQ1*PMQ2)))
1156           ZH=1D0+PMQ1-PMQ2
1157           ZL=0.5D0*(ZH-ZD)
1158           ZU=0.5D0*(ZH+ZD)
1159           IF(V(IM,1).LT.ZL.OR.V(IM,1).GT.ZU) GOTO 350
1160         ENDIF
1161       ENDIF
1162  
1163 C...Accepted branch. Construct four-momentum for initial partons.
1164   490 MAZIP=0
1165       MAZIC=0
1166       IF(NEP.EQ.1) THEN
1167         P(N+1,1)=0D0
1168         P(N+1,2)=0D0
1169         P(N+1,3)=SQRT(MAX(0D0,(P(IPA(1),4)+P(N+1,5))*(P(IPA(1),4)-
1170      &  P(N+1,5))))
1171         P(N+1,4)=P(IPA(1),4)
1172         V(N+1,2)=P(N+1,4)
1173       ELSEIF(IGM.EQ.0.AND.NEP.EQ.2) THEN
1174         PED1=0.5D0*(V(IM,5)+V(N+1,5)-V(N+2,5))/P(IM,5)
1175         P(N+1,1)=0D0
1176         P(N+1,2)=0D0
1177         P(N+1,3)=SQRT(MAX(0D0,(PED1+P(N+1,5))*(PED1-P(N+1,5))))
1178         P(N+1,4)=PED1
1179         P(N+2,1)=0D0
1180         P(N+2,2)=0D0
1181         P(N+2,3)=-P(N+1,3)
1182         P(N+2,4)=P(IM,5)-PED1
1183         V(N+1,2)=P(N+1,4)
1184         V(N+2,2)=P(N+2,4)
1185       ELSEIF(NEP.GE.3) THEN
1186 C...Rescale all momenta for energy conservation.
1187         LOOP=0
1188         PES=0D0
1189         PQS=0D0
1190         DO 510 I=1,NEP
1191           DO 500 J=1,4
1192             P(N+I,J)=P(IPA(I),J)
1193   500     CONTINUE
1194           PES=PES+P(N+I,4)
1195           PQS=PQS+P(N+I,5)**2/P(N+I,4)
1196   510   CONTINUE
1197   520   LOOP=LOOP+1
1198         FAC=(PS(5)-PQS)/(PES-PQS)
1199         PES=0D0
1200         PQS=0D0
1201         DO 540 I=1,NEP
1202           DO 530 J=1,3
1203             P(N+I,J)=FAC*P(N+I,J)
1204   530     CONTINUE
1205           P(N+I,4)=SQRT(P(N+I,5)**2+P(N+I,1)**2+P(N+I,2)**2+P(N+I,3)**2)
1206           V(N+I,2)=P(N+I,4)
1207           PES=PES+P(N+I,4)
1208           PQS=PQS+P(N+I,5)**2/P(N+I,4)
1209   540   CONTINUE
1210         IF(LOOP.LT.10.AND.ABS(PES-PS(5)).GT.1D-12*PS(5)) GOTO 520
1211  
1212 C...Construct transverse momentum for ordinary branching in shower.
1213       ELSE
1214         ZM=V(IM,1)
1215         LOOPPT=0
1216   550   LOOPPT=LOOPPT+1
1217         PZM=SQRT(MAX(0D0,(PEM+P(IM,5))*(PEM-P(IM,5))))
1218         PMLS=(V(IM,5)-V(N+1,5)-V(N+2,5))**2-4D0*V(N+1,5)*V(N+2,5)
1219         IF(PZM.LE.0D0) THEN
1220           PTS=0D0
1221         ELSEIF(K(IM,2).EQ.21.AND.IABS(K(N+1,2)).LE.10.AND.
1222      &  (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
1223           PTS=PMLS*ZM*(1D0-ZM)/V(IM,5)
1224         ELSEIF(MOD(MSTJ(43),2).EQ.1) THEN
1225           PTS=(PEM**2*(ZM*(1D0-ZM)*V(IM,5)-(1D0-ZM)*V(N+1,5)-
1226      &    ZM*V(N+2,5))-0.25D0*PMLS)/PZM**2
1227         ELSE
1228           PTS=PMLS*(ZM*(1D0-ZM)*PEM**2/V(IM,5)-0.25D0)/PZM**2
1229         ENDIF
1230         IF(PTS.LT.0D0.AND.LOOPPT.LT.10) THEN
1231           ZM=0.05D0+0.9D0*ZM
1232           GOTO 550
1233         ELSEIF(PTS.LT.0D0) THEN
1234           GOTO 280
1235         ENDIF
1236         PT=SQRT(MAX(0D0,PTS))
1237  
1238 C...Global statistics.
1239         MINT(353)=MINT(353)+1
1240         VINT(353)=VINT(353)+PT
1241         IF (MINT(353).EQ.1) VINT(358)=PT
1242  
1243 C...Find coefficient of azimuthal asymmetry due to gluon polarization.
1244         HAZIP=0D0
1245         IF(MSTJ(49).NE.1.AND.MOD(MSTJ(46),2).EQ.1.AND.K(IM,2).EQ.21
1246      &  .AND.IAU.NE.0) THEN
1247           IF(K(IGM,3).NE.0) MAZIP=1
1248           ZAU=V(IGM,1)
1249           IF(IAU.EQ.IM+1) ZAU=1D0-V(IGM,1)
1250           IF(MAZIP.EQ.0) ZAU=0D0
1251           IF(K(IGM,2).NE.21) THEN
1252             HAZIP=2D0*ZAU/(1D0+ZAU**2)
1253           ELSE
1254             HAZIP=(ZAU/(1D0-ZAU*(1D0-ZAU)))**2
1255           ENDIF
1256           IF(K(N+1,2).NE.21) THEN
1257             HAZIP=HAZIP*(-2D0*ZM*(1D0-ZM))/(1D0-2D0*ZM*(1D0-ZM))
1258           ELSE
1259             HAZIP=HAZIP*(ZM*(1D0-ZM)/(1D0-ZM*(1D0-ZM)))**2
1260           ENDIF
1261         ENDIF
1262  
1263 C...Find coefficient of azimuthal asymmetry due to soft gluon
1264 C...interference.
1265         HAZIC=0D0
1266         IF(MSTJ(49).NE.2.AND.MSTJ(46).GE.2.AND.(K(N+1,2).EQ.21.OR.
1267      &  K(N+2,2).EQ.21).AND.IAU.NE.0) THEN
1268           IF(K(IGM,3).NE.0) MAZIC=N+1
1269           IF(K(IGM,3).NE.0.AND.K(N+1,2).NE.21) MAZIC=N+2
1270           IF(K(IGM,3).NE.0.AND.K(N+1,2).EQ.21.AND.K(N+2,2).EQ.21.AND.
1271      &    ZM.GT.0.5D0) MAZIC=N+2
1272           IF(K(IAU,2).EQ.22) MAZIC=0
1273           ZS=ZM
1274           IF(MAZIC.EQ.N+2) ZS=1D0-ZM
1275           ZGM=V(IGM,1)
1276           IF(IAU.EQ.IM-1) ZGM=1D0-V(IGM,1)
1277           IF(MAZIC.EQ.0) ZGM=1D0
1278           IF(MAZIC.NE.0) HAZIC=(P(IM,5)/P(IGM,5))*
1279      &    SQRT((1D0-ZS)*(1D0-ZGM)/(ZS*ZGM))
1280           HAZIC=MIN(0.95D0,HAZIC)
1281         ENDIF
1282       ENDIF
1283  
1284 C...Construct energies for ordinary branching in shower.
1285   560 IF(NEP.EQ.2.AND.IGM.GT.0) THEN
1286         IF(K(IM,2).EQ.21.AND.IABS(K(N+1,2)).LE.10.AND.
1287      &  (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
1288           P(N+1,4)=0.5D0*(PEM*(V(IM,5)+V(N+1,5)-V(N+2,5))+
1289      &    PZM*SQRT(MAX(0D0,PMLS))*(2D0*ZM-1D0))/V(IM,5)
1290         ELSEIF(MOD(MSTJ(43),2).EQ.1) THEN
1291           P(N+1,4)=PEM*V(IM,1)
1292         ELSE
1293           P(N+1,4)=PEM*(0.5D0*(V(IM,5)-SQRT(PMLS)+V(N+1,5)-V(N+2,5))+
1294      &    SQRT(PMLS)*ZM)/V(IM,5)
1295         ENDIF
1296  
1297 C...Already predetermined choice of phi angle or not
1298         PHI=PARU(2)*PYR(0)
1299         IF(MPSPD.EQ.1.AND.IGM.EQ.NS+1) THEN
1300           IPSPD=IP1+IM-NS-2
1301           IF(K(IPSPD,4).GT.0) THEN
1302             IPSGD1=K(IPSPD,4)
1303             IF(IM.EQ.NS+2) THEN
1304               PHI=PYANGL(P(IPSGD1,1),P(IPSGD1,2))
1305             ELSE
1306               PHI=PYANGL(-P(IPSGD1,1),P(IPSGD1,2))
1307             ENDIF
1308           ENDIF
1309         ELSEIF(MPSPD.EQ.1.AND.IGM.EQ.NS+2) THEN
1310           IPSPD=IP1+IM-NS-2
1311           IF(K(IPSPD,4).GT.0) THEN
1312             IPSGD1=K(IPSPD,4)
1313             PHIPSM=PYANGL(P(IPSPD,1),P(IPSPD,2))
1314             THEPSM=PYANGL(P(IPSPD,3),SQRT(P(IPSPD,1)**2+P(IPSPD,2)**2))
1315             CALL PYROBO(IPSGD1,IPSGD1,0D0,-PHIPSM,0D0,0D0,0D0)
1316             CALL PYROBO(IPSGD1,IPSGD1,-THEPSM,0D0,0D0,0D0,0D0)
1317             PHI=PYANGL(P(IPSGD1,1),P(IPSGD1,2))
1318             CALL PYROBO(IPSGD1,IPSGD1,THEPSM,PHIPSM,0D0,0D0,0D0)
1319           ENDIF
1320         ENDIF
1321  
1322 C...Construct momenta for ordinary branching in shower.
1323         P(N+1,1)=PT*COS(PHI)
1324         P(N+1,2)=PT*SIN(PHI)
1325         IF(K(IM,2).EQ.21.AND.IABS(K(N+1,2)).LE.10.AND.
1326      &  (MSTJ(44).EQ.3.OR.MSTJ(44).EQ.5)) THEN
1327           P(N+1,3)=0.5D0*(PZM*(V(IM,5)+V(N+1,5)-V(N+2,5))+
1328      &    PEM*SQRT(MAX(0D0,PMLS))*(2D0*ZM-1D0))/V(IM,5)
1329         ELSEIF(PZM.GT.0D0) THEN
1330           P(N+1,3)=0.5D0*(V(N+2,5)-V(N+1,5)-V(IM,5)+
1331      &    2D0*PEM*P(N+1,4))/PZM
1332         ELSE
1333           P(N+1,3)=0D0
1334         ENDIF
1335         P(N+2,1)=-P(N+1,1)
1336         P(N+2,2)=-P(N+1,2)
1337         P(N+2,3)=PZM-P(N+1,3)
1338         P(N+2,4)=PEM-P(N+1,4)
1339         IF(MSTJ(43).LE.2) THEN
1340           V(N+1,2)=(PEM*P(N+1,4)-PZM*P(N+1,3))/P(IM,5)
1341           V(N+2,2)=(PEM*P(N+2,4)-PZM*P(N+2,3))/P(IM,5)
1342         ENDIF
1343       ENDIF
1344  
1345 C...Rotate and boost daughters.
1346       IF(IGM.GT.0) THEN
1347         IF(MSTJ(43).LE.2) THEN
1348           BEX=P(IGM,1)/P(IGM,4)
1349           BEY=P(IGM,2)/P(IGM,4)
1350           BEZ=P(IGM,3)/P(IGM,4)
1351           GA=P(IGM,4)/P(IGM,5)
1352           GABEP=GA*(GA*(BEX*P(IM,1)+BEY*P(IM,2)+BEZ*P(IM,3))/(1D0+GA)-
1353      &    P(IM,4))
1354         ELSE
1355           BEX=0D0
1356           BEY=0D0
1357           BEZ=0D0
1358           GA=1D0
1359           GABEP=0D0
1360         ENDIF
1361         PTIMB=SQRT((P(IM,1)+GABEP*BEX)**2+(P(IM,2)+GABEP*BEY)**2)
1362         THE=PYANGL(P(IM,3)+GABEP*BEZ,PTIMB)
1363         IF(PTIMB.GT.1D-4) THEN
1364           PHI=PYANGL(P(IM,1)+GABEP*BEX,P(IM,2)+GABEP*BEY)
1365         ELSE
1366           PHI=0D0
1367         ENDIF
1368         DO 570 I=N+1,N+2
1369           DP(1)=COS(THE)*COS(PHI)*P(I,1)-SIN(PHI)*P(I,2)+
1370      &    SIN(THE)*COS(PHI)*P(I,3)
1371           DP(2)=COS(THE)*SIN(PHI)*P(I,1)+COS(PHI)*P(I,2)+
1372      &    SIN(THE)*SIN(PHI)*P(I,3)
1373           DP(3)=-SIN(THE)*P(I,1)+COS(THE)*P(I,3)
1374           DP(4)=P(I,4)
1375           DBP=BEX*DP(1)+BEY*DP(2)+BEZ*DP(3)
1376           DGABP=GA*(GA*DBP/(1D0+GA)+DP(4))
1377           P(I,1)=DP(1)+DGABP*BEX
1378           P(I,2)=DP(2)+DGABP*BEY
1379           P(I,3)=DP(3)+DGABP*BEZ
1380           P(I,4)=GA*(DP(4)+DBP)
1381   570   CONTINUE
1382       ENDIF
1383  
1384 C...Weight with azimuthal distribution, if required.
1385       IF(MAZIP.NE.0.OR.MAZIC.NE.0) THEN
1386         DO 580 J=1,3
1387           DPT(1,J)=P(IM,J)
1388           DPT(2,J)=P(IAU,J)
1389           DPT(3,J)=P(N+1,J)
1390   580   CONTINUE
1391         DPMA=DPT(1,1)*DPT(2,1)+DPT(1,2)*DPT(2,2)+DPT(1,3)*DPT(2,3)
1392         DPMD=DPT(1,1)*DPT(3,1)+DPT(1,2)*DPT(3,2)+DPT(1,3)*DPT(3,3)
1393         DPMM=DPT(1,1)**2+DPT(1,2)**2+DPT(1,3)**2
1394         DO 590 J=1,3
1395           DPT(4,J)=DPT(2,J)-DPMA*DPT(1,J)/MAX(1D-10,DPMM)
1396           DPT(5,J)=DPT(3,J)-DPMD*DPT(1,J)/MAX(1D-10,DPMM)
1397   590   CONTINUE
1398         DPT(4,4)=SQRT(DPT(4,1)**2+DPT(4,2)**2+DPT(4,3)**2)
1399         DPT(5,4)=SQRT(DPT(5,1)**2+DPT(5,2)**2+DPT(5,3)**2)
1400         IF(MIN(DPT(4,4),DPT(5,4)).GT.0.1D0*PARJ(82)) THEN
1401           CAD=(DPT(4,1)*DPT(5,1)+DPT(4,2)*DPT(5,2)+
1402      &    DPT(4,3)*DPT(5,3))/(DPT(4,4)*DPT(5,4))
1403           IF(MAZIP.NE.0) THEN
1404             IF(1D0+HAZIP*(2D0*CAD**2-1D0).LT.PYR(0)*(1D0+ABS(HAZIP)))
1405      &      GOTO 560
1406           ENDIF
1407           IF(MAZIC.NE.0) THEN
1408             IF(MAZIC.EQ.N+2) CAD=-CAD
1409             IF((1D0-HAZIC)*(1D0-HAZIC*CAD)/(1D0+HAZIC**2-2D0*HAZIC*CAD)
1410      &      .LT.PYR(0)) GOTO 560
1411           ENDIF
1412         ENDIF
1413       ENDIF
1414  
1415 C...Azimuthal anisotropy due to interference with initial state partons.
1416       IF(MOD(MIIS,2).EQ.1.AND.IGM.EQ.NS+1.AND.(K(N+1,2).EQ.21.OR.
1417      &K(N+2,2).EQ.21)) THEN
1418         III=IM-NS-1
1419         IF(ISII(III).GE.1) THEN
1420           IAZIID=N+1
1421           IF(K(N+1,2).NE.21) IAZIID=N+2
1422           IF(K(N+1,2).EQ.21.AND.K(N+2,2).EQ.21.AND.
1423      &    P(N+1,4).GT.P(N+2,4)) IAZIID=N+2
1424           THEIID=PYANGL(P(IAZIID,3),SQRT(P(IAZIID,1)**2+P(IAZIID,2)**2))
1425           IF(III.EQ.2) THEIID=PARU(1)-THEIID
1426           PHIIID=PYANGL(P(IAZIID,1),P(IAZIID,2))
1427           HAZII=MIN(0.95D0,THEIID/THEIIS(III,ISII(III)))
1428           CAD=COS(PHIIID-PHIIIS(III,ISII(III)))
1429           PHIREL=ABS(PHIIID-PHIIIS(III,ISII(III)))
1430           IF(PHIREL.GT.PARU(1)) PHIREL=PARU(2)-PHIREL
1431           IF((1D0-HAZII)*(1D0-HAZII*CAD)/(1D0+HAZII**2-2D0*HAZII*CAD)
1432      &    .LT.PYR(0)) GOTO 560
1433         ENDIF
1434       ENDIF
1435  
1436 C...Continue loop over partons that may branch, until none left.
1437       IF(IGM.GE.0) K(IM,1)=14
1438       N=N+NEP
1439       NEP=2
1440       IF(N.GT.MSTU(4)-MSTU(32)-10) THEN
1441         CALL PYERRM(11,'(PYSHOW:) no more memory left in PYJETS')
1442         IF(MSTU(21).GE.1) N=NS
1443         IF(MSTU(21).GE.1) RETURN
1444       ENDIF
1445       GOTO 290
1446  
1447 C...Set information on imagined shower initiator.
1448   600 IF(NPA.GE.2) THEN
1449         K(NS+1,1)=11
1450         K(NS+1,2)=94
1451         K(NS+1,3)=IP1
1452         IF(IP2.GT.0.AND.IP2.LT.IP1) K(NS+1,3)=IP2
1453         K(NS+1,4)=NS+2
1454         K(NS+1,5)=NS+1+NPA
1455         IIM=1
1456       ELSE
1457         IIM=0
1458       ENDIF
1459  
1460 C...Reconstruct string drawing information.
1461       DO 610 I=NS+1+IIM,N
1462         KQ=KCHG(PYCOMP(K(I,2)),2)
1463         IF(K(I,1).LE.10.AND.K(I,2).EQ.22) THEN
1464           K(I,1)=1
1465         ELSEIF(K(I,1).LE.10.AND.IABS(K(I,2)).GE.11.AND.
1466      &    IABS(K(I,2)).LE.18) THEN
1467           K(I,1)=1
1468         ELSEIF(K(I,1).LE.10) THEN
1469           K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))
1470           K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))
1471         ELSEIF(K(MOD(K(I,4),MSTU(5))+1,2).NE.22) THEN
1472           ID1=MOD(K(I,4),MSTU(5))
1473           IF(KQ.EQ.1.AND.K(I,2).GT.0) ID1=MOD(K(I,4),MSTU(5))+1
1474           IF(KQ.EQ.2.AND.(K(ID1,2).EQ.21.OR.K(ID1+1,2).EQ.21).AND.
1475      &    PYR(0).GT.0.5D0) ID1=MOD(K(I,4),MSTU(5))+1
1476           ID2=2*MOD(K(I,4),MSTU(5))+1-ID1
1477           K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1
1478           K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID2
1479           K(ID1,4)=K(ID1,4)+MSTU(5)*I
1480           K(ID1,5)=K(ID1,5)+MSTU(5)*ID2
1481           K(ID2,4)=K(ID2,4)+MSTU(5)*ID1
1482           K(ID2,5)=K(ID2,5)+MSTU(5)*I
1483         ELSE
1484           ID1=MOD(K(I,4),MSTU(5))
1485           ID2=ID1+1
1486           K(I,4)=MSTU(5)*(K(I,4)/MSTU(5))+ID1
1487           K(I,5)=MSTU(5)*(K(I,5)/MSTU(5))+ID1
1488           IF(KQ.EQ.1.OR.K(ID1,1).GE.11) THEN
1489             K(ID1,4)=K(ID1,4)+MSTU(5)*I
1490             K(ID1,5)=K(ID1,5)+MSTU(5)*I
1491           ELSE
1492             K(ID1,4)=0
1493             K(ID1,5)=0
1494           ENDIF
1495           K(ID2,4)=0
1496           K(ID2,5)=0
1497         ENDIF
1498   610 CONTINUE
1499  
1500 C...Transformation from CM frame.
1501       IF(NPA.EQ.1) THEN
1502         THE=PYANGL(P(IPA(1),3),SQRT(P(IPA(1),1)**2+P(IPA(1),2)**2))
1503         PHI=PYANGL(P(IPA(1),1),P(IPA(1),2))
1504         MSTU(33)=1
1505         CALL PYROBO(NS+1,N,THE,PHI,0D0,0D0,0D0)
1506       ELSEIF(NPA.EQ.2) THEN
1507         BEX=PS(1)/PS(4)
1508         BEY=PS(2)/PS(4)
1509         BEZ=PS(3)/PS(4)
1510         GA=PS(4)/PS(5)
1511         GABEP=GA*(GA*(BEX*P(IPA(1),1)+BEY*P(IPA(1),2)+BEZ*P(IPA(1),3))
1512      &  /(1D0+GA)-P(IPA(1),4))
1513         THE=PYANGL(P(IPA(1),3)+GABEP*BEZ,SQRT((P(IPA(1),1)
1514      &  +GABEP*BEX)**2+(P(IPA(1),2)+GABEP*BEY)**2))
1515         PHI=PYANGL(P(IPA(1),1)+GABEP*BEX,P(IPA(1),2)+GABEP*BEY)
1516         MSTU(33)=1
1517         CALL PYROBO(NS+1,N,THE,PHI,BEX,BEY,BEZ)
1518       ELSE
1519         CALL PYROBO(IPA(1),IPA(NPA),0D0,0D0,PS(1)/PS(4),PS(2)/PS(4),
1520      &  PS(3)/PS(4))
1521         MSTU(33)=1
1522         CALL PYROBO(NS+1,N,0D0,0D0,PS(1)/PS(4),PS(2)/PS(4),PS(3)/PS(4))
1523       ENDIF
1524  
1525 C...Decay vertex of shower.
1526       DO 630 I=NS+1,N
1527         DO 620 J=1,5
1528           V(I,J)=V(IP1,J)
1529   620   CONTINUE
1530   630 CONTINUE
1531  
1532 C...Delete trivial shower, else connect initiators.
1533       IF(N.LE.NS+NPA+IIM) THEN
1534         N=NS
1535       ELSE
1536         DO 640 IP=1,NPA
1537           K(IPA(IP),1)=14
1538           K(IPA(IP),4)=K(IPA(IP),4)+NS+IIM+IP
1539           K(IPA(IP),5)=K(IPA(IP),5)+NS+IIM+IP
1540           K(NS+IIM+IP,3)=IPA(IP)
1541           IF(IIM.EQ.1.AND.MSTU(16).NE.2) K(NS+IIM+IP,3)=NS+1
1542           IF(K(NS+IIM+IP,1).NE.1) THEN
1543             K(NS+IIM+IP,4)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,4)
1544             K(NS+IIM+IP,5)=MSTU(5)*IPA(IP)+K(NS+IIM+IP,5)
1545           ENDIF
1546   640   CONTINUE
1547       ENDIF
1548  
1549       RETURN
1550       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!