ISAJET/code/qcdjet.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE QCDJET(NJMIN)
   3 C
   4 C          Carry out final state QCD jet evolution using the algorithm
   5 C          of Fox and Wolfram.  Evolve each parton in T with fixed ZC
   6 C          and iterate as follows--
   7 C
   8 C          (0) Evolve initial partons.
   9 C          (1) Pick I and find matching J>I.
  10 C          (2) Solve kinematics.
  11 C          (3) For K=I,J, generate Z(K) and evolve T(K1), T(K2).  If no
  12 C              good, evolve T(K). Otherwise, add K1 and K2 to /JETSET/.
  13 C          (4) If I or J no good, then (2).
  14 C          (5) Then (1).
  15 C
  16 C          Use Z=(E+P)/(E0+P0) and a large TCUT.
  17 C          JMATCH(J1)=J2 if J1 and J2 match.
  18 C          JMATCH(J)=JPACK*J1+J2 if J1,...,J2 match. Used for multiple
  19 C          initial partons.
  20 C          JMATCH(J)=0 for initial jet partons
  21 C
  22 C          Include W+- and Z0 radiation.
  23 C
  24 #if defined(CERNLIB_IMPNONE)
  25       IMPLICIT NONE
  26 #endif
  27 #include "isajet/itapes.inc"
  28 #include "isajet/partcl.inc"
  29 #include "isajet/qcdpar.inc"
  30 #include "isajet/jetset.inc"
  31 #include "isajet/jwork.inc"
  32 #include "isajet/const.inc"
  33 C
  34       INTEGER J,NJMIN,JPRNT,JI1,JI2,NJI,JI,NJ1,NJ2,L,K,NPTLV1,IFAIL,J0
  35       REAL AM0,AM1,AM2,RANF,AMSUM,PCM2,POLD2,RATIO,PSUM,P12CM,E0,P0,Z1,
  36      $E1MAX,P1MAX,ZMAX,E1MIN,P1MIN,ZMIN,ZEP,E1,P1,CTHCMZ,Z2,E2MAX,P2MAX,
  37      $E2MIN,P2MIN,P2,E2,CTHCM,STHCM,PHICM,CPHICM,SPHICM,PT0,CTH0,STH0,
  38      $CPHI0,SPHI0,SGN,BP,ZLIM,ZLIM1
  39       DIMENSION PSUM(5)
  40       DATA PSUM/5*0./
  41 C
  42 C          (0) Evolve initial parton masses.
  43 C
  44       DO 100 J=NJMIN,NJSET
  45         J1=J
  46         J2=JMATCH(J)
  47         IF(J2.GT.JPACK) GO TO 150
  48         IF(J2.LE.J1) GO TO 100
  49 C          Two partons
  50         IF(JDCAY(J1).EQ.-1) CALL QCDT(J1)
  51         IF(JDCAY(J2).EQ.-1) CALL QCDT(J2)
  52         JPRNT=MOD(JORIG(J),JPACK)
  53         IF(JPRNT.EQ.0) THEN
  54           AM0=PJSET(4,J1)+PJSET(4,J2)
  55         ELSE
  56           AM0=PJSET(5,JPRNT)
  57         ENDIF
  58 110     AM1=PJSET(5,J1)
  59         AM2=PJSET(5,J2)
  60         IF(AM0.LE.AM1+AM2) THEN
  61           J3=J1
  62           IF(RANF().GT..5) J3=J2
  63           IF(JDCAY(J3).EQ.-1) CALL QCDT(J3)
  64           GO TO 110
  65         ENDIF
  66         GO TO 100
  67 C          More than two partons
  68 150     JI1=JMATCH(J)/JPACK
  69         IF(J.NE.JI1) GO TO 100
  70         JI2=JMATCH(J)-JPACK*JI1
  71         NJI=JI2-JI1+1
  72         AM0=0.
  73         AMSUM=0.
  74         DO 160 JI=JI1,JI2
  75           IF(JDCAY(JI).EQ.-1) CALL QCDT(JI)
  76           AM0=AM0+PJSET(4,JI)
  77           AMSUM=AMSUM+PJSET(5,JI)
  78 160     CONTINUE
  79 170     IF(AM0.LT.AMSUM) THEN
  80           J3=NJI*RANF()+JI1
  81           AMSUM=AMSUM-PJSET(5,J3)
  82           IF(JDCAY(J3).EQ.-1) CALL QCDT(J3)
  83           AMSUM=AMSUM+PJSET(5,J3)
  84           GO TO 170
  85         ENDIF
  86 100   CONTINUE
  87 C
  88 C          (1) Loop over active partons
  89 C
  90       NJ1=NJMIN
  91 1     NJ2=NJSET
  92       DO 200 J=NJ1,NJ2
  93         J1=J
  94         J2=JMATCH(J1)
  95         NJI=2
  96         IF(J2.LE.J1) GO TO 200
  97 C
  98 C          (2) Solve kinematics.
  99 C
 100 C          Initial partons--keep directions fixed.
 101 210     IF(MOD(JORIG(J),JPACK).NE.0) GO TO 230
 102         IF(JMATCH(J).GT.JPACK) GO TO 400
 103         AM0=PJSET(4,J1)+PJSET(4,J2)
 104         AM1=PJSET(5,J1)
 105         AM2=PJSET(5,J2)
 106         PJSET(4,J1)=(AM0**2+AM1**2-AM2**2)/(2*AM0)
 107         PJSET(4,J2)=(AM0**2+AM2**2-AM1**2)/(2*AM0)
 108         PCM2=((AM0**2-AM1**2-AM2**2)**2-(2*AM1*AM2)**2)/(4*AM0**2)
 109         DO 220 L=1,2
 110           POLD2=PJSET(1,JJ(L))**2+PJSET(2,JJ(L))**2+PJSET(3,JJ(L))**2
 111           RATIO=SQRT(PCM2/POLD2)
 112           DO 225 K=1,3
 113 225       PJSET(K,JJ(L))=RATIO*PJSET(K,JJ(L))
 114 220     CONTINUE
 115         GO TO 300
 116 C
 117 C          NJI.LE.5 initial partons
 118 400     CONTINUE
 119         JI1=JMATCH(J)/JPACK
 120         IF(J.NE.JI1) GO TO 200
 121         JI2=JMATCH(J)-JPACK*JI1
 122         NJI=JI2-JI1+1
 123         AM0=0.
 124         DO 410 JI=JI1,JI2
 125           AM0=AM0+PJSET(4,JI)
 126           JJ(JI-JI1+1)=JI
 127           PJSET(4,JI)=SQRT(PJSET(1,JI)**2+PJSET(2,JI)**2+PJSET(3,JI)**2
 128      1    +PJSET(5,JI)**2)
 129           DO 420 K=1,5
 130 420       PPTCL(K,NPTCL+JI-JI1+1)=PJSET(K,JI)
 131 410     CONTINUE
 132         PSUM(4)=AM0
 133         PSUM(5)=PSUM(4)
 134         NPTLV1=NPTCL
 135         CALL RESCAL(NPTLV1+1,NPTLV1+NJI,PSUM,IFAIL)
 136         DO 430 JI=JI1,JI2
 137         DO 430 K=1,5
 138           PJSET(K,JI)=PPTCL(K,NPTCL+JI-JI1+1)
 139 430     CONTINUE
 140         GO TO 300
 141 C
 142 C          Solve kinematics for general partons.
 143 C
 144 230     J0=MOD(JORIG(J),JPACK)
 145         AM0=PJSET(5,J0)
 146         AM1=PJSET(5,J1)
 147         AM2=PJSET(5,J2)
 148         E1CM=(AM0**2+AM1**2-AM2**2)/(2*AM0)
 149         E2CM=(AM0**2+AM2**2-AM1**2)/(2*AM0)
 150         P12CM=SQRT((AM0**2-AM1**2-AM2**2)**2-(2*AM1*AM2)**2)/(2*AM0)
 151         NJI=2
 152 C          Determine E1, P1, and COS(THCM) from Z(J0).
 153 C          Occasionally COS(TH)>1.  If so then reset Z.
 154         E0=PJSET(4,J0)
 155         P0=SQRT(PJSET(1,J0)**2+PJSET(2,J0)**2+PJSET(3,J0)**2)
 156         Z1=ZZC(J0)
 157         IF(Z1.GT.0.5) THEN
 158           E1MAX=(E0*E1CM+P0*P12CM)/AM0
 159           P1MAX=(P0*E1CM+E0*P12CM)/AM0
 160           ZMAX=(E1MAX+P1MAX)/(E0+P0)
 161           E1MIN=(E0*E1CM-P0*P12CM)/AM0
 162           P1MIN=(P0*E1CM-E0*P12CM)/AM0
 163           P1MIN=ABS(P1MIN)
 164           ZMIN=(E1MIN+P1MIN)/(E0+P0)
 165           IF(Z1.LT.ZMIN.OR.Z1.GT.ZMAX) Z1=ZMIN+Z1*(ZMAX-ZMIN)
 166           ZZC(J0)=Z1
 167           ZEP=Z1*(E0+P0)
 168           P1=(ZEP**2-AM1**2)/(2.*ZEP)
 169           E1=(ZEP**2+AM1**2)/(2.*ZEP)
 170           CTHCM=(E1*AM0-E0*E1CM)/(P0*P12CM)
 171         ELSE
 172           Z2=1.-Z1
 173           E2MAX=(E0*E2CM+P0*P12CM)/AM0
 174           P2MAX=(P0*E2CM+E0*P12CM)/AM0
 175           ZMAX=(E2MAX+P2MAX)/(E0+P0)
 176           E2MIN=(E0*E2CM-P0*P12CM)/AM0
 177           P2MIN=(P0*E2CM-E0*P12CM)/AM0
 178           P2MIN=ABS(P2MIN)
 179           ZMIN=(E2MIN+P2MIN)/(E0+P0)
 180           IF(Z2.LT.ZMIN.OR.Z2.GT.ZMAX) Z2=ZMIN+Z2*(ZMAX-ZMIN)
 181           ZZC(J0)=Z2
 182           ZEP=Z2*(E0+P0)
 183           P2=(ZEP**2-AM2**2)/(2.*ZEP)
 184           E2=(ZEP**2+AM2**2)/(2.*ZEP)
 185           CTHCM=-(E2*AM0-E0*E2CM)/(P0*P12CM)
 186         ENDIF
 187 C          Avoid disaster
 188         IF(ABS(CTHCM).GT.1.) CTHCM=SIGN(RANF(),CTHCM)
 189         STHCM=SQRT(1.-CTHCM**2)
 190         PHICM=2*PI*RANF()
 191         CPHICM=COS(PHICM)
 192         SPHICM=SIN(PHICM)
 193 C
 194 C          Construct cm momenta.
 195         PT0=SQRT(PJSET(1,J0)**2+PJSET(2,J0)**2)
 196         CTH0=PJSET(3,J0)/P0
 197         STH0=PT0/P0
 198         CPHI0=PJSET(1,J0)/PT0
 199         SPHI0=PJSET(2,J0)/PT0
 200         P1CM(1)=P12CM*(CPHI0*(CTH0*CPHICM*STHCM+STH0*CTHCM)
 201      1  -SPHI0*SPHICM*STHCM)
 202         P1CM(2)=P12CM*(SPHI0*(CTH0*CPHICM*STHCM+STH0*CTHCM)
 203      1  +CPHI0*SPHICM*STHCM)
 204         P1CM(3)=P12CM*(-STH0*CPHICM*STHCM+CTH0*CTHCM)
 205 C          Boost with P0 to get lab momenta
 206         DO 240 L=1,2
 207           SGN=3-2*L
 208           BP=0
 209           DO 241 K=1,3
 210 241       BP=BP+PJSET(K,J0)*SGN*P1CM(K)
 211           BP=BP/AM0
 212           PJSET(4,JJ(L))=PJSET(4,J0)*EE(L)/PJSET(5,J0)+BP
 213           DO 242 K=1,3
 214 242       PJSET(K,JJ(L))=SGN*P1CM(K)+PJSET(K,J0)*EE(L)/PJSET(5,J0)
 215      1    +PJSET(K,J0)*BP/(PJSET(4,J0)+PJSET(5,J0))
 216 240     CONTINUE
 217 C
 218 C          (3) Pick Z and decay partons. Check.
 219 C
 220 300     CONTINUE
 221         TNEW=.FALSE.
 222         DO 310 L=1,NJI
 223           IF(JDCAY(JJ(L)).GE.0) GO TO 310
 224           IF(NJSET+2.GT.MXJSET) GO TO 9999
 225           CALL QCDZ(JJ(L))
 226           CALL QCDT(NJSET+1)
 227           CALL QCDT(NJSET+2)
 228 C
 229 C          Check whether masses allowed.
 230           AM0=PJSET(5,JJ(L))
 231           AM1=PJSET(5,NJSET+1)
 232           AM2=PJSET(5,NJSET+2)
 233           IF(AM1+AM2.GE.AM0) GO TO 320
 234 C
 235 C          Check whether Z allowed.
 236           E1CM=(AM0**2+AM1**2-AM2**2)/(2*AM0)
 237           E2CM=(AM0**2+AM2**2-AM1**2)/(2.*AM0)
 238           P12CM=SQRT((AM0**2-AM1**2-AM2**2)**2-(2*AM1*AM2)**2)/(2*AM0)
 239           E0=PJSET(4,JJ(L))
 240           P0=SQRT(PJSET(1,JJ(L))**2+PJSET(2,JJ(L))**2+PJSET(3,JJ(L))**2)
 241           IF(ZZC(JJ(L)).GT.0.5) THEN
 242             ZEP=ZZC(JJ(L))*(E0+P0)
 243             P1=(ZEP**2-AM1**2)/(2.*ZEP)
 244             E1=(ZEP**2+AM1**2)/(2.*ZEP)
 245             CTHCM=(E1*AM0-E0*E1CM)/(P0*P12CM)
 246             IF((ABS(CTHCM).GE.1..OR.P1.LE.0.).AND.IABS(JTYPE(JJ(L)))
 247      $      .LT.80) GO TO 320
 248           ELSE
 249             ZEP=(1.-ZZC(JJ(L)))*(E0+P0)
 250             P2=(ZEP**2-AM2**2)/(2.*ZEP)
 251             E2=(ZEP**2+AM2**2)/(2.*ZEP)
 252             CTHCM=-(E2*AM0-E0*E2CM)/(P0*P12CM)
 253             IF((ABS(CTHCM).GE.1..OR.P2.LE.0.).AND.IABS(JTYPE(JJ(L)))
 254      $      .LT.80) GO TO 320
 255           ENDIF
 256 C
 257 C          Require Z and 1-Z within kinematic limits.
 258 C
 259           ZLIM=(AM0/(E0+P0))**2
 260           ZLIM1=CUTJET/(E0+P0)
 261           ZLIM=AMAX1(ZLIM,ZLIM1)
 262           IF((ZZC(JJ(L)).GT.ZLIM.AND.ZZC(JJ(L)).LT.(1.-ZLIM)).OR.
 263      $    IABS(JTYPE(JJ(L))).GE.80) THEN
 264 C          Add new partons to /JETSET/.
 265             JDCAY(JJ(L))=JPACK*(NJSET+1)+(NJSET+2)
 266             NJSET=NJSET+2
 267             GO TO 310
 268           ENDIF
 269 C          Discard partons and evolve JJ(L) again.
 270 320       TNEW=.TRUE.
 271           CALL QCDT(JJ(L))
 272 310     CONTINUE
 273 C
 274 C          (4) Resolve kinematics if any parton mass is changed.
 275 C
 276         IF(TNEW) GO TO 210
 277 200   CONTINUE
 278 C
 279 C          (5) Iterate entire proceedure.
 280 C
 281       NJ1=NJ2+1
 282       IF(NJ1.LE.NJSET) GO TO 1
 283       RETURN
 284 C
 285 C          Error message
 286 C
 287 9999  CALL PRTEVT(0)
 288       WRITE(ITLIS,10) NJSET
 289 10    FORMAT(//' ERROR IN QCDJET...NJSET > ',I4)
 290       RETURN
 291       END