| 0795afa3 |
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 |
git://git.uio.no / u / mrichter / AliRoot.git / blame