1 #include "isajet/pilot.h"
2 LOGICAL FUNCTION DECJET(IP,NADD,IDABS,PREST,WDECAY,BETA,GAMMA)
4 C Auxiliary routine for DECAY. Evolve and hadronize partons.
5 C Check conservation laws. Return TRUE if OK, FALSE otherwise.
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.
14 #if defined(CERNLIB_IMPNONE)
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"
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
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)
39 C Copy decay products into /JETSET/ and do QCD evolution.
41 IF(NJSET+NADD.GT.MXJSET) GO TO 9998
48 110 PJSET(K,NJSET)=PREST(K,I)
49 PJSET(5,NJSET)=PPTCL(5,NPTCL+I)
51 JTYPE(NJSET)=IDENT(NPTCL+I)
53 JMATCH(NJSET)=JPACK*(NJSAVE+1)+NJSAVE+NADD
56 C For heavy quarks match 1+2 and 3+(1+2). Boost 1+2 to rest.
59 JMATCH(NJSAVE+1)=NJSAVE+2
60 JMATCH(NJSAVE+2)=NJSAVE+1
63 PWREST(K)=PJSET(K,NJSAVE+1)+PJSET(K,NJSAVE+2)
64 PJSET(K,NJSET)=PWREST(K)
66 PWREST(5)=SQRT(PWREST(4)**2-PWREST(1)**2-PWREST(2)**2
68 PJSET(5,NJSET)=PWREST(5)
69 JMATCH(NJSAVE+3)=NJSAVE+4
70 JMATCH(NJSAVE+4)=NJSAVE+3
75 CHARGW=CHARGW+CHARGE(IDLV1)
76 JTYPE(NJSAVE+4)=80*SIGN(1.,CHARGW)
78 C Boost W vectors to rest.
80 130 BETAW(K)=PWREST(K)/PWREST(4)
81 GAMMAW=PWREST(4)/PWREST(5)
84 DO 140 J=NJSAV1,NJSAV2
85 BP=BETAW(1)*PJSET(1,J)+BETAW(2)*PJSET(2,J)+BETAW(3)*PJSET(3,J)
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)
93 C Do evolution and save new W momentum. Start from parent
94 C mass or NADD*energy.
97 IF(IABS(JTYPE(J)).LT.10.OR.MOD(JTYPE(J),100).EQ.0) THEN
99 PJSET(5,J)=AMIN1(PPTCL(5,IP),NADD*PJSET(4,J))
103 CALL QCDJET(NJSAVE+1)
106 PWREST(4)=PJSET(4,NJSAVE+4)
107 GAMMAW=PWREST(4)/PWREST(5)
109 PWREST(K)=PJSET(K,NJSAVE+4)
110 BETAW(K)=PWREST(K)/PWREST(4)
114 C Put final partons in particle table - temporary IORIG.
115 C Also include virtual or real W momentum for quark decays.
123 210 PPTCL(K,NPTCL)=PJSET(K,NJSAVE+4)
125 IDENT(NPTCL)=JTYPE(NJSAVE+4)
130 IF(JDCAY(J).NE.0) GO TO 220
131 IF(JORIG(J)/JPACK.NE.3) GO TO 220
134 221 PPTCL(K,NPTCL)=PJSET(K,J)
135 IORIG(NPTCL)=3*IPACK+IP
136 IDENT(NPTCL)=JTYPE(J)
144 IF(JDCAY(J).NE.0) GO TO 240
145 IF(JORIG(J)/JPACK.NE.JET) GO TO 240
147 BP=BETAW(1)*PJSET(1,J)+BETAW(2)*PJSET(2,J)
148 $ +BETAW(3)*PJSET(3,J)
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)
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
163 C Not quark decay, so just copy partons.
167 IF(JDCAY(J).NE.0) GO TO 260
168 IF(JORIG(J)/JPACK.NE.JET) GO TO 260
171 261 PPTCL(K,NPTCL)=PJSET(K,J)
172 IORIG(NPTCL)=IPACK*(JORIG(J)/JPACK)+IP
173 IDENT(NPTCL)=JTYPE(J)
177 IDCAY(IP)=NSTART*IPACK+NPTCL
181 C Hadronize quarks and rotate to proper angles.
186 IF(JDCAY(I).NE.0) GO TO 310
187 IF(JORIG(I)/JPACK.NE.JET) GO TO 310
189 IF(IABS(JTYPE(I)).GE.10.AND.MOD(JTYPE(I),100).NE.0)
194 PJET=SQRT(PJSET(1,I)**2+PJSET(2,I)**2+PJSET(3,I)**2)
195 CTHQK=PJSET(3,I)/PJET
199 CPHIQK=PJSET(1,I)/(PJET*STHQK)
200 SPHIQK=PJSET(2,I)/(PJET*STHQK)
207 IF(NEXT.GT.NPTCL) GO TO 310
208 ROT(1,1)=CPHIQK*CTHQK
209 ROT(2,1)=SPHIQK*CTHQK
214 ROT(1,3)=CPHIQK*STHQK
215 ROT(2,3)=SPHIQK*STHQK
225 322 PPTCL(K,II)=PPTCL(K,II)+ROT(K,KK)*PSAVE(KK)
226 IORIG(II)=IPACK*JET+NPRTN
229 IDCAY(NPRTN)=NEXT*IPACK+NPTCL
235 331 PPTCL(K,NPTCL)=PJSET(K,I)
236 IORIG(NPTCL)=IPACK*JET+NPRTN
237 IDENT(NPTCL)=JTYPE(I)
239 IDCAY(NPRTN)=NPTCL*IPACK+NPTCL
244 C Reset NJSET so decay jets do not appear in /JETSET/
248 C Check for at least two particles
249 IF(NPTCL.LT.NHDRN+2) THEN
259 DO 400 I=NHDRN1,NPTCL
261 SUMQ=SUMQ+CHARGE(IDLV1)
264 SUMQ=SUMQ-CHARGE(IDLV1)
266 IF(SUMQ.EQ.0.) GO TO 500
268 C Charge wrong--fix it by swapping UP and DN quarks.
269 DO 410 I=NHDRN1,NPTCL
271 IF(IABS(ID1).GT.1000) GO TO 410
272 I1=MOD(IABS(ID1)/100,10)
273 I2=MOD(IABS(ID1)/10,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
281 ELSEIF(I1.EQ.1.AND.I2.EQ.1) THEN
282 IDENT(I)=(120+I3)*(-SIGN(1.,SUMQ))
286 SUMQ=SIGN(ABS(SUMQ)-1.,SUMQ)
288 PPTCL(5,I)=AMASS(IDLV1)
289 PPTCL(4,I)=SQRT(PPTCL(1,I)**2+PPTCL(2,I)**2+PPTCL(3,I)**2
291 C Sum cannot vanish for fractionally charged initial particle.
292 IF(ABS(SUMQ).LT.1.) GO TO 500
294 C Failed to conserve charge.
299 C Rescale momenta for correct mass
303 C Quark decay. First rescale jet3 + W
305 510 PPTCL(K,NPTCL+1)=PPTCL(K,NPTCLW)
311 CALL RESCAL(NPJET(2)+1,NPTLV1,PSUM,IFAIL)
318 530 BETAW(K)=PPTCL(K,NPTCL+1)/PPTCL(4,NPTCL+1)
319 GAMMAW=PPTCL(4,NPTCL+1)/PPTCL(5,NPTCL+1)
321 PSUM(4)=PPTCL(5,NPTCLW)
323 CALL RESCAL(NHDRN1,NPJET(2),PSUM,IFAIL)
330 C General decay with no W.
336 CALL RESCAL(NHDRN1,NPTLV1,PSUM,IFAIL)
344 C Boost back to lab frame. Reset IORIG.
347 DO 600 I=NHDRN1,NPTCL
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)
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)
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)
362 PPTCL(J,I)=PPTCL(J,I)+GAMMA*BETA(J)*(PPTCL(4,I)
363 $ +BP*GAMMA/(GAMMA+1.))
365 PPTCL(4,I)=GAMMA*(PPTCL(4,I)+BP)
377 WRITE(ITLIS,99980) NJSET
378 99980 FORMAT(//5X,'ERROR IN DECJET...NJSET > ',I5)