]>
Commit | Line | Data |
---|---|---|
0795afa3 | 1 | #include "isajet/pilot.h" |
2 | LOGICAL FUNCTION DECJET(IP,NADD,IDABS,PREST,WDECAY,BETA,GAMMA) | |
3 | C | |
4 | C Auxiliary routine for DECAY. Evolve and hadronize partons. | |
5 | C Check conservation laws. Return TRUE if OK, FALSE otherwise. | |
6 | C | |
7 | C IP = particle to be decayed. | |
8 | C NADD = number of products (NPTCL+1, ..., NPTCL+NADD). | |
9 | C IDABS = absolute values of decay IDENT's. | |
10 | C PREST = 4-momenta in rest frame. | |
11 | C WDECAY = logical flag for real W decay. | |
12 | C BETA,GAMMA = boost parameters. | |
13 | C | |
14 | #if defined(CERNLIB_IMPNONE) | |
15 | IMPLICIT NONE | |
16 | #endif | |
17 | #include "isajet/itapes.inc" | |
18 | #include "isajet/wcon.inc" | |
19 | #include "isajet/partcl.inc" | |
20 | #include "isajet/dkytab.inc" | |
21 | #include "isajet/jetset.inc" | |
22 | #include "isajet/jwork.inc" | |
23 | #include "isajet/const.inc" | |
24 | C | |
25 | REAL PGEN(5,5),RND(5),U(3),BETA(3),IDQK(3),ROT(3,3),PSAVE(3) | |
26 | 1,REDUCE(5),WPROP,Z,TRY,RANF,AMASS,TWOME,CHARGE | |
27 | REAL PSUM(5),POLD(4),PNEW(4),SUM,WTMAX,SUM1,SUM2 | |
28 | REAL PREST(4,6),PWREST(5),BETAW(3),DOT,PCM | |
29 | REAL AMEE,REE,WTEE,SWAP,RNEW,WT,QCM,PHI,S12,S12MAX,GAMMAW,BP | |
30 | REAL PJET,CTHQK,STHQK,CPHIQK,SPHIQK,SUMQ,A,B,C,GAMMA | |
31 | REAL CHARGW | |
32 | LOGICAL WDECAY | |
33 | INTEGER IDLV1,IFL1,IFL2,IFL3,JSPIN,INDEX,IPOINT,ID1,I1,I2,I3 | |
34 | INTEGER NADD,NSTART,NEW,NADD1,J,IP,I,IDABS(5),NEXT | |
35 | INTEGER JJ1,II,K1,K,NJSAVE,NJSAV1,NJSAV2,NJ1,NPRTN,KK,NHDRN1 | |
36 | INTEGER IFAIL,JSAVE,JETIP,JET,NJADD,NPTLV1,IDANTI,NPJET(5) | |
37 | INTEGER NHDRN,NPJET3,NPTCLW,NPBEG(5) | |
38 | C | |
39 | C Copy decay products into /JETSET/ and do QCD evolution. | |
40 | C | |
41 | IF(NJSET+NADD.GT.MXJSET) GO TO 9998 | |
42 | NJSAVE=NJSET | |
43 | NSTART=NPTCL-NADD+1 | |
44 | NPTCL=NSTART-1 | |
45 | DO 100 I=1,NADD | |
46 | NJSET=NJSET+1 | |
47 | DO 110 K=1,4 | |
48 | 110 PJSET(K,NJSET)=PREST(K,I) | |
49 | PJSET(5,NJSET)=PPTCL(5,NPTCL+I) | |
50 | JORIG(NJSET)=JPACK*I | |
51 | JTYPE(NJSET)=IDENT(NPTCL+I) | |
52 | JDCAY(NJSET)=0 | |
53 | JMATCH(NJSET)=JPACK*(NJSAVE+1)+NJSAVE+NADD | |
54 | 100 CONTINUE | |
55 | C | |
56 | C For heavy quarks match 1+2 and 3+(1+2). Boost 1+2 to rest. | |
57 | C | |
58 | IF(WDECAY) THEN | |
59 | JMATCH(NJSAVE+1)=NJSAVE+2 | |
60 | JMATCH(NJSAVE+2)=NJSAVE+1 | |
61 | NJSET=NJSET+1 | |
62 | DO 120 K=1,4 | |
63 | PWREST(K)=PJSET(K,NJSAVE+1)+PJSET(K,NJSAVE+2) | |
64 | PJSET(K,NJSET)=PWREST(K) | |
65 | 120 CONTINUE | |
66 | PWREST(5)=SQRT(PWREST(4)**2-PWREST(1)**2-PWREST(2)**2 | |
67 | $ -PWREST(3)**2) | |
68 | PJSET(5,NJSET)=PWREST(5) | |
69 | JMATCH(NJSAVE+3)=NJSAVE+4 | |
70 | JMATCH(NJSAVE+4)=NJSAVE+3 | |
71 | JORIG(NJSAVE+4)=-1 | |
72 | IDLV1=JTYPE(NJSAVE+1) | |
73 | CHARGW=CHARGE(IDLV1) | |
74 | IDLV1=JTYPE(NJSAVE+2) | |
75 | CHARGW=CHARGW+CHARGE(IDLV1) | |
76 | JTYPE(NJSAVE+4)=80*SIGN(1.,CHARGW) | |
77 | JDCAY(NJSAVE+4)=0 | |
78 | C Boost W vectors to rest. | |
79 | DO 130 K=1,3 | |
80 | 130 BETAW(K)=PWREST(K)/PWREST(4) | |
81 | GAMMAW=PWREST(4)/PWREST(5) | |
82 | NJSAV1=NJSAVE+1 | |
83 | NJSAV2=NJSAVE+2 | |
84 | DO 140 J=NJSAV1,NJSAV2 | |
85 | BP=BETAW(1)*PJSET(1,J)+BETAW(2)*PJSET(2,J)+BETAW(3)*PJSET(3,J) | |
86 | DO 141 K=1,3 | |
87 | 141 PJSET(K,J)=PJSET(K,J)-GAMMAW*BETAW(K)*(PJSET(4,J) | |
88 | $ -BP*GAMMAW/(GAMMAW+1.)) | |
89 | PJSET(4,J)=GAMMAW*(PJSET(4,J)-BP) | |
90 | 140 CONTINUE | |
91 | ENDIF | |
92 | C | |
93 | C Do evolution and save new W momentum. Start from parent | |
94 | C mass or NADD*energy. | |
95 | NJSAV1=NJSAVE+1 | |
96 | DO 150 J=NJSAV1,NJSET | |
97 | IF(IABS(JTYPE(J)).LT.10.OR.MOD(JTYPE(J),100).EQ.0) THEN | |
98 | JDCAY(J)=-1 | |
99 | PJSET(5,J)=AMIN1(PPTCL(5,IP),NADD*PJSET(4,J)) | |
100 | ENDIF | |
101 | 150 CONTINUE | |
102 | C | |
103 | CALL QCDJET(NJSAVE+1) | |
104 | C | |
105 | IF(WDECAY) THEN | |
106 | PWREST(4)=PJSET(4,NJSAVE+4) | |
107 | GAMMAW=PWREST(4)/PWREST(5) | |
108 | DO 200 K=1,3 | |
109 | PWREST(K)=PJSET(K,NJSAVE+4) | |
110 | BETAW(K)=PWREST(K)/PWREST(4) | |
111 | 200 CONTINUE | |
112 | ENDIF | |
113 | C | |
114 | C Put final partons in particle table - temporary IORIG. | |
115 | C Also include virtual or real W momentum for quark decays. | |
116 | C | |
117 | NJ1=NJSAVE+1 | |
118 | IF(WDECAY) THEN | |
119 | C Real or virtual W. | |
120 | NPTCL=NPTCL+1 | |
121 | NPTCLW=NPTCL | |
122 | DO 210 K=1,5 | |
123 | 210 PPTCL(K,NPTCL)=PJSET(K,NJSAVE+4) | |
124 | IORIG(NPTCL)=IP | |
125 | IDENT(NPTCL)=JTYPE(NJSAVE+4) | |
126 | IDCAY(NPTCL)=0 | |
127 | C Jet 3 | |
128 | NPBEG(3)=NPTCL+1 | |
129 | DO 220 J=NJ1,NJSET | |
130 | IF(JDCAY(J).NE.0) GO TO 220 | |
131 | IF(JORIG(J)/JPACK.NE.3) GO TO 220 | |
132 | NPTCL=NPTCL+1 | |
133 | DO 221 K=1,5 | |
134 | 221 PPTCL(K,NPTCL)=PJSET(K,J) | |
135 | IORIG(NPTCL)=3*IPACK+IP | |
136 | IDENT(NPTCL)=JTYPE(J) | |
137 | IDCAY(NPTCL)=0 | |
138 | 220 CONTINUE | |
139 | C Jets 1 and 2 | |
140 | NPJET3=NPTCL | |
141 | DO 230 JET=1,2 | |
142 | NPBEG(JET)=NPTCL+1 | |
143 | DO 240 J=NJ1,NJSET | |
144 | IF(JDCAY(J).NE.0) GO TO 240 | |
145 | IF(JORIG(J)/JPACK.NE.JET) GO TO 240 | |
146 | NPTCL=NPTCL+1 | |
147 | BP=BETAW(1)*PJSET(1,J)+BETAW(2)*PJSET(2,J) | |
148 | $ +BETAW(3)*PJSET(3,J) | |
149 | DO 241 K=1,3 | |
150 | 241 PPTCL(K,NPTCL)=PJSET(K,J)+GAMMAW*BETAW(K)*(PJSET(4,J) | |
151 | $ +BP*GAMMAW/(GAMMAW+1.)) | |
152 | PPTCL(4,NPTCL)=GAMMAW*(PJSET(4,J)+BP) | |
153 | PPTCL(5,NPTCL)=PJSET(5,J) | |
154 | IORIG(NPTCL)=IPACK*(JORIG(J)/JPACK)+NPTCLW | |
155 | IDENT(NPTCL)=JTYPE(J) | |
156 | IDCAY(NPTCL)=0 | |
157 | 240 CONTINUE | |
158 | 230 CONTINUE | |
159 | C Quark decays to W plus jet 3; then W decays. | |
160 | IDCAY(IP)=IPACK*NPTCLW+NPJET3 | |
161 | IDCAY(NPTCLW)=IPACK*(NPJET3+1)+NPTCL | |
162 | ELSE | |
163 | C Not quark decay, so just copy partons. | |
164 | DO 250 JET=1,NADD | |
165 | NPBEG(JET)=NPTCL+1 | |
166 | DO 260 J=NJ1,NJSET | |
167 | IF(JDCAY(J).NE.0) GO TO 260 | |
168 | IF(JORIG(J)/JPACK.NE.JET) GO TO 260 | |
169 | NPTCL=NPTCL+1 | |
170 | DO 261 K=1,5 | |
171 | 261 PPTCL(K,NPTCL)=PJSET(K,J) | |
172 | IORIG(NPTCL)=IPACK*(JORIG(J)/JPACK)+IP | |
173 | IDENT(NPTCL)=JTYPE(J) | |
174 | IDCAY(NPTCL)=0 | |
175 | 260 CONTINUE | |
176 | 250 CONTINUE | |
177 | IDCAY(IP)=NSTART*IPACK+NPTCL | |
178 | ENDIF | |
179 | NHDRN=NPTCL | |
180 | C | |
181 | C Hadronize quarks and rotate to proper angles. | |
182 | C | |
183 | DO 300 JET=1,NADD | |
184 | NPRTN=NPBEG(JET)-1 | |
185 | DO 310 I=NJ1,NJSET | |
186 | IF(JDCAY(I).NE.0) GO TO 310 | |
187 | IF(JORIG(I)/JPACK.NE.JET) GO TO 310 | |
188 | NPRTN=NPRTN+1 | |
189 | IF(IABS(JTYPE(I)).GE.10.AND.MOD(JTYPE(I),100).NE.0) | |
190 | $ GO TO 330 | |
191 | C | |
192 | C Fragment parton: | |
193 | NEXT=NPTCL+1 | |
194 | PJET=SQRT(PJSET(1,I)**2+PJSET(2,I)**2+PJSET(3,I)**2) | |
195 | CTHQK=PJSET(3,I)/PJET | |
196 | STHQK=1.-CTHQK**2 | |
197 | IF(STHQK.LT.1) THEN | |
198 | STHQK=SQRT(STHQK) | |
199 | CPHIQK=PJSET(1,I)/(PJET*STHQK) | |
200 | SPHIQK=PJSET(2,I)/(PJET*STHQK) | |
201 | ELSE | |
202 | STHQK=0 | |
203 | CPHIQK=1 | |
204 | SPHIQK=0 | |
205 | ENDIF | |
206 | CALL JETGEN(I) | |
207 | IF(NEXT.GT.NPTCL) GO TO 310 | |
208 | ROT(1,1)=CPHIQK*CTHQK | |
209 | ROT(2,1)=SPHIQK*CTHQK | |
210 | ROT(3,1)=-STHQK | |
211 | ROT(1,2)=-SPHIQK | |
212 | ROT(2,2)=CPHIQK | |
213 | ROT(3,2)=0. | |
214 | ROT(1,3)=CPHIQK*STHQK | |
215 | ROT(2,3)=SPHIQK*STHQK | |
216 | ROT(3,3)=CTHQK | |
217 | C | |
218 | DO 320 II=NEXT,NPTCL | |
219 | DO 321 K=1,3 | |
220 | PSAVE(K)=PPTCL(K,II) | |
221 | PPTCL(K,II)=0. | |
222 | 321 CONTINUE | |
223 | DO 322 K=1,3 | |
224 | DO 322 KK=1,3 | |
225 | 322 PPTCL(K,II)=PPTCL(K,II)+ROT(K,KK)*PSAVE(KK) | |
226 | IORIG(II)=IPACK*JET+NPRTN | |
227 | IDCAY(II)=0 | |
228 | 320 CONTINUE | |
229 | IDCAY(NPRTN)=NEXT*IPACK+NPTCL | |
230 | GO TO 310 | |
231 | C | |
232 | C or add lepton: | |
233 | 330 NPTCL=NPTCL+1 | |
234 | DO 331 K=1,5 | |
235 | 331 PPTCL(K,NPTCL)=PJSET(K,I) | |
236 | IORIG(NPTCL)=IPACK*JET+NPRTN | |
237 | IDENT(NPTCL)=JTYPE(I) | |
238 | IDCAY(NPTCL)=0 | |
239 | IDCAY(NPRTN)=NPTCL*IPACK+NPTCL | |
240 | 310 CONTINUE | |
241 | NPJET(JET)=NPTCL | |
242 | 300 CONTINUE | |
243 | C | |
244 | C Reset NJSET so decay jets do not appear in /JETSET/ | |
245 | NJADD=NJSET | |
246 | NJSET=NJSAVE | |
247 | C | |
248 | C Check for at least two particles | |
249 | IF(NPTCL.LT.NHDRN+2) THEN | |
250 | NPTCL=NSTART-1 | |
251 | DECJET=.FALSE. | |
252 | RETURN | |
253 | ENDIF | |
254 | C | |
255 | C Conserve charge | |
256 | C | |
257 | SUMQ=0. | |
258 | NHDRN1=NHDRN+1 | |
259 | DO 400 I=NHDRN1,NPTCL | |
260 | IDLV1=IDENT(I) | |
261 | SUMQ=SUMQ+CHARGE(IDLV1) | |
262 | 400 CONTINUE | |
263 | IDLV1=IDENT(IP) | |
264 | SUMQ=SUMQ-CHARGE(IDLV1) | |
265 | C | |
266 | IF(SUMQ.EQ.0.) GO TO 500 | |
267 | C | |
268 | C Charge wrong--fix it by swapping UP and DN quarks. | |
269 | DO 410 I=NHDRN1,NPTCL | |
270 | ID1=IDENT(I) | |
271 | IF(IABS(ID1).GT.1000) GO TO 410 | |
272 | I1=MOD(IABS(ID1)/100,10) | |
273 | I2=MOD(IABS(ID1)/10,10) | |
274 | I3=MOD(IABS(ID1),10) | |
275 | IF(I1.EQ.1.AND.I2.GT.2.AND.SUMQ*ID1.GT.0.) THEN | |
276 | IDENT(I)=ISIGN(200+10*I2+I3,ID1) | |
277 | ELSEIF(I1.EQ.2.AND.I2.GT.2.AND.SUMQ*ID1.LT.0.) THEN | |
278 | IDENT(I)=ISIGN(100+10*I2+I3,ID1) | |
279 | ELSEIF(I1.EQ.1.AND.I2.EQ.2.AND.SUMQ*ID1.GT.0.) THEN | |
280 | IDENT(I)=110+I3 | |
281 | ELSEIF(I1.EQ.1.AND.I2.EQ.1) THEN | |
282 | IDENT(I)=(120+I3)*(-SIGN(1.,SUMQ)) | |
283 | ELSE | |
284 | GO TO 410 | |
285 | ENDIF | |
286 | SUMQ=SIGN(ABS(SUMQ)-1.,SUMQ) | |
287 | IDLV1=IDENT(I) | |
288 | PPTCL(5,I)=AMASS(IDLV1) | |
289 | PPTCL(4,I)=SQRT(PPTCL(1,I)**2+PPTCL(2,I)**2+PPTCL(3,I)**2 | |
290 | $ +PPTCL(5,I)**2) | |
291 | C Sum cannot vanish for fractionally charged initial particle. | |
292 | IF(ABS(SUMQ).LT.1.) GO TO 500 | |
293 | 410 CONTINUE | |
294 | C Failed to conserve charge. | |
295 | NPTCL=NSTART-1 | |
296 | DECJET=.FALSE. | |
297 | RETURN | |
298 | C | |
299 | C Rescale momenta for correct mass | |
300 | C | |
301 | 500 CONTINUE | |
302 | IF(WDECAY) THEN | |
303 | C Quark decay. First rescale jet3 + W | |
304 | DO 510 K=1,5 | |
305 | 510 PPTCL(K,NPTCL+1)=PPTCL(K,NPTCLW) | |
306 | NPTLV1=NPTCL+1 | |
307 | DO 520 K=1,3 | |
308 | 520 PSUM(K)=0. | |
309 | PSUM(4)=PPTCL(5,IP) | |
310 | PSUM(5)=PSUM(4) | |
311 | CALL RESCAL(NPJET(2)+1,NPTLV1,PSUM,IFAIL) | |
312 | IF(IFAIL.NE.0) THEN | |
313 | NPTCL=NSTART-1 | |
314 | DECJET=.FALSE. | |
315 | RETURN | |
316 | ENDIF | |
317 | DO 530 K=1,3 | |
318 | 530 BETAW(K)=PPTCL(K,NPTCL+1)/PPTCL(4,NPTCL+1) | |
319 | GAMMAW=PPTCL(4,NPTCL+1)/PPTCL(5,NPTCL+1) | |
320 | C Then rescale W | |
321 | PSUM(4)=PPTCL(5,NPTCLW) | |
322 | PSUM(5)=PSUM(4) | |
323 | CALL RESCAL(NHDRN1,NPJET(2),PSUM,IFAIL) | |
324 | IF(IFAIL.NE.0) THEN | |
325 | NPTCL=NSTART-1 | |
326 | DECJET=.FALSE. | |
327 | RETURN | |
328 | ENDIF | |
329 | ELSE | |
330 | C General decay with no W. | |
331 | DO 550 K=1,3 | |
332 | 550 PSUM(K)=0. | |
333 | PSUM(4)=PPTCL(5,IP) | |
334 | PSUM(5)=PSUM(4) | |
335 | NPTLV1=NPTCL | |
336 | CALL RESCAL(NHDRN1,NPTLV1,PSUM,IFAIL) | |
337 | IF(IFAIL.NE.0) THEN | |
338 | NPTCL=NSTART-1 | |
339 | DECJET=.FALSE. | |
340 | RETURN | |
341 | ENDIF | |
342 | ENDIF | |
343 | C | |
344 | C Boost back to lab frame. Reset IORIG. | |
345 | C | |
346 | IF(WDECAY) THEN | |
347 | DO 600 I=NHDRN1,NPTCL | |
348 | JET=IORIG(I)/IPACK | |
349 | IF(JET.NE.1.AND.JET.NE.2) GO TO 600 | |
350 | BP=BETAW(1)*PPTCL(1,I)+BETAW(2)*PPTCL(2,I)+BETAW(3)*PPTCL(3,I) | |
351 | DO 610 J=1,3 | |
352 | 610 PPTCL(J,I)=PPTCL(J,I)+GAMMAW*BETAW(J)*(PPTCL(4,I) | |
353 | $ +BP*GAMMAW/(GAMMAW+1.)) | |
354 | PPTCL(4,I)=GAMMAW*(PPTCL(4,I)+BP) | |
355 | 600 CONTINUE | |
356 | ENDIF | |
357 | C | |
358 | DO 620 I=NSTART,NPTCL | |
359 | IORIG(I)=MOD(IORIG(I),IPACK) | |
360 | BP=BETA(1)*PPTCL(1,I)+BETA(2)*PPTCL(2,I)+BETA(3)*PPTCL(3,I) | |
361 | DO 621 J=1,3 | |
362 | PPTCL(J,I)=PPTCL(J,I)+GAMMA*BETA(J)*(PPTCL(4,I) | |
363 | $ +BP*GAMMA/(GAMMA+1.)) | |
364 | 621 CONTINUE | |
365 | PPTCL(4,I)=GAMMA*(PPTCL(4,I)+BP) | |
366 | 620 CONTINUE | |
367 | C | |
368 | C Normal exit | |
369 | C | |
370 | DECJET=.TRUE. | |
371 | RETURN | |
372 | C | |
373 | C Error messages. | |
374 | C | |
375 | 9998 DECJET=.FALSE. | |
376 | CALL PRTEVT(0) | |
377 | WRITE(ITLIS,99980) NJSET | |
378 | 99980 FORMAT(//5X,'ERROR IN DECJET...NJSET > ',I5) | |
379 | RETURN | |
380 | END |