]>
Commit | Line | Data |
---|---|---|
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 |