]>
Commit | Line | Data |
---|---|---|
e74335a4 | 1 | * $Id$ |
2 | ||
3 | C********************************************************************* | |
4 | ||
5 | SUBROUTINE LUDECY_HIJING(IP) | |
6 | ||
7 | C...Purpose: to handle the decay of unstable particles. | |
8 | #include "lujets_hijing.inc" | |
9 | #include "ludat1_hijing.inc" | |
10 | #include "ludat2_hijing.inc" | |
11 | #include "ludat3_hijing.inc" | |
12 | DIMENSION VDCY(4),KFLO(4),KFL1(4),PV(10,5),RORD(10),UE(3),BE(3), | |
13 | &WTCOR(10) | |
14 | DATA WTCOR/2.,5.,15.,60.,250.,1500.,1.2E4,1.2E5,150.,16./ | |
15 | ||
16 | C...Functions: momentum in two-particle decays, four-product and | |
17 | C...matrix element times phase space in weak decays. | |
18 | PAWT(A,B,C)=SQRT((A**2-(B+C)**2)*(A**2-(B-C)**2))/(2.*A) | |
19 | FOUR(I,J)=P(I,4)*P(J,4)-P(I,1)*P(J,1)-P(I,2)*P(J,2)-P(I,3)*P(J,3) | |
20 | HMEPS(HA)=((1.-HRQ-HA)**2+3.*HA*(1.+HRQ-HA))* | |
21 | &SQRT((1.-HRQ-HA)**2-4.*HRQ*HA) | |
22 | ||
23 | C...Initial values. | |
24 | NTRY=0 | |
25 | NSAV=N | |
26 | KFA=IABS(K(IP,2)) | |
27 | KFS=ISIGN(1,K(IP,2)) | |
28 | KC=LUCOMP_HIJING(KFA) | |
29 | MSTJ(92)=0 | |
30 | ||
31 | C...Choose lifetime and determine decay vertex. | |
32 | IF(K(IP,1).EQ.5) THEN | |
33 | V(IP,5)=0. | |
34 | ELSEIF(K(IP,1).NE.4) THEN | |
35 | V(IP,5)=-PMAS(KC,4)*LOG(RLU_HIJING(0)) | |
36 | ENDIF | |
37 | DO 100 J=1,4 | |
38 | 100 VDCY(J)=V(IP,J)+V(IP,5)*P(IP,J)/P(IP,5) | |
39 | ||
40 | C...Determine whether decay allowed or not. | |
41 | MOUT=0 | |
42 | IF(MSTJ(22).EQ.2) THEN | |
43 | IF(PMAS(KC,4).GT.PARJ(71)) MOUT=1 | |
44 | ELSEIF(MSTJ(22).EQ.3) THEN | |
45 | IF(VDCY(1)**2+VDCY(2)**2+VDCY(3)**2.GT.PARJ(72)**2) MOUT=1 | |
46 | ELSEIF(MSTJ(22).EQ.4) THEN | |
47 | IF(VDCY(1)**2+VDCY(2)**2.GT.PARJ(73)**2) MOUT=1 | |
48 | IF(ABS(VDCY(3)).GT.PARJ(74)) MOUT=1 | |
49 | ENDIF | |
50 | IF(MOUT.EQ.1.AND.K(IP,1).NE.5) THEN | |
51 | K(IP,1)=4 | |
52 | RETURN | |
53 | ENDIF | |
54 | ||
55 | C...Check existence of decay channels. Particle/antiparticle rules. | |
56 | KCA=KC | |
57 | IF(MDCY(KC,2).GT.0) THEN | |
58 | MDMDCY=MDME(MDCY(KC,2),2) | |
59 | IF(MDMDCY.GT.80.AND.MDMDCY.LE.90) KCA=MDMDCY | |
60 | ENDIF | |
61 | IF(MDCY(KCA,2).LE.0.OR.MDCY(KCA,3).LE.0) THEN | |
62 | CALL LUERRM_HIJING(9 | |
63 | $ ,'(LUDECY_HIJING:) no decay channel defined') | |
64 | RETURN | |
65 | ENDIF | |
66 | IF(MOD(KFA/1000,10).EQ.0.AND.(KCA.EQ.85.OR.KCA.EQ.87)) KFS=-KFS | |
67 | IF(KCHG(KC,3).EQ.0) THEN | |
68 | KFSP=1 | |
69 | KFSN=0 | |
70 | IF(RLU_HIJING(0).GT.0.5) KFS=-KFS | |
71 | ELSEIF(KFS.GT.0) THEN | |
72 | KFSP=1 | |
73 | KFSN=0 | |
74 | ELSE | |
75 | KFSP=0 | |
76 | KFSN=1 | |
77 | ENDIF | |
78 | ||
79 | C...Sum branching ratios of allowed decay channels. | |
80 | 110 NOPE=0 | |
81 | BRSU=0. | |
82 | DO 120 IDL=MDCY(KCA,2),MDCY(KCA,2)+MDCY(KCA,3)-1 | |
83 | IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND. | |
84 | &KFSN*MDME(IDL,1).NE.3) GOTO 120 | |
85 | IF(MDME(IDL,2).GT.100) GOTO 120 | |
86 | NOPE=NOPE+1 | |
87 | BRSU=BRSU+BRAT(IDL) | |
88 | 120 CONTINUE | |
89 | IF(NOPE.EQ.0) THEN | |
90 | CALL LUERRM_HIJING(2 | |
91 | $ ,'(LUDECY_HIJING:) all decay channels closed by user') | |
92 | RETURN | |
93 | ENDIF | |
94 | ||
95 | C...Select decay channel among allowed ones. | |
96 | 130 RBR=BRSU*RLU_HIJING(0) | |
97 | IDL=MDCY(KCA,2)-1 | |
98 | 140 IDL=IDL+1 | |
99 | IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND. | |
100 | &KFSN*MDME(IDL,1).NE.3) THEN | |
101 | IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 140 | |
102 | ELSEIF(MDME(IDL,2).GT.100) THEN | |
103 | IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 140 | |
104 | ELSE | |
105 | IDC=IDL | |
106 | RBR=RBR-BRAT(IDL) | |
107 | IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1.AND.RBR.GT.0.) GOTO 140 | |
108 | ENDIF | |
109 | ||
110 | C...Start readout of decay channel: matrix element, reset counters. | |
111 | MMAT=MDME(IDC,2) | |
112 | 150 NTRY=NTRY+1 | |
113 | IF(NTRY.GT.1000) THEN | |
114 | CALL LUERRM_HIJING(14 | |
115 | $ ,'(LUDECY_HIJING:) caught in infinite loop') | |
116 | IF(MSTU(21).GE.1) RETURN | |
117 | ENDIF | |
118 | I=N | |
119 | NP=0 | |
120 | NQ=0 | |
121 | MBST=0 | |
122 | IF(MMAT.GE.11.AND.MMAT.NE.46.AND.P(IP,4).GT.20.*P(IP,5)) MBST=1 | |
123 | DO 160 J=1,4 | |
124 | PV(1,J)=0. | |
125 | 160 IF(MBST.EQ.0) PV(1,J)=P(IP,J) | |
126 | IF(MBST.EQ.1) PV(1,4)=P(IP,5) | |
127 | PV(1,5)=P(IP,5) | |
128 | PS=0. | |
129 | PSQ=0. | |
130 | MREM=0 | |
131 | ||
132 | C...Read out decay products. Convert to standard flavour code. | |
133 | JTMAX=5 | |
134 | IF(MDME(IDC+1,2).EQ.101) JTMAX=10 | |
135 | DO 170 JT=1,JTMAX | |
136 | IF(JT.LE.5) KP=KFDP(IDC,JT) | |
137 | IF(JT.GE.6) KP=KFDP(IDC+1,JT-5) | |
138 | IF(KP.EQ.0) GOTO 170 | |
139 | KPA=IABS(KP) | |
140 | KCP=LUCOMP_HIJING(KPA) | |
141 | IF(KCHG(KCP,3).EQ.0.AND.KPA.NE.81.AND.KPA.NE.82) THEN | |
142 | KFP=KP | |
143 | ELSEIF(KPA.NE.81.AND.KPA.NE.82) THEN | |
144 | KFP=KFS*KP | |
145 | ELSEIF(KPA.EQ.81.AND.MOD(KFA/1000,10).EQ.0) THEN | |
146 | KFP=-KFS*MOD(KFA/10,10) | |
147 | ELSEIF(KPA.EQ.81.AND.MOD(KFA/100,10).GE.MOD(KFA/10,10)) THEN | |
148 | KFP=KFS*(100*MOD(KFA/10,100)+3) | |
149 | ELSEIF(KPA.EQ.81) THEN | |
150 | KFP=KFS*(1000*MOD(KFA/10,10)+100*MOD(KFA/100,10)+1) | |
151 | ELSEIF(KP.EQ.82) THEN | |
152 | CALL LUKFDI_HIJING(-KFS*INT(1.+(2.+PARJ(2))*RLU_HIJING(0)),0 | |
153 | $ ,KFP,KDUMP) | |
154 | IF(KFP.EQ.0) GOTO 150 | |
155 | MSTJ(93)=1 | |
156 | IF(PV(1,5).LT.PARJ(32)+2.*ULMASS_HIJING(KFP)) GOTO 150 | |
157 | ELSEIF(KP.EQ.-82) THEN | |
158 | KFP=-KFP | |
159 | IF(IABS(KFP).GT.10) KFP=KFP+ISIGN(10000,KFP) | |
160 | ENDIF | |
161 | IF(KPA.EQ.81.OR.KPA.EQ.82) KCP=LUCOMP_HIJING(KFP) | |
162 | ||
163 | C...Add decay product to event record or to quark flavour list. | |
164 | KFPA=IABS(KFP) | |
165 | KQP=KCHG(KCP,2) | |
166 | IF(MMAT.GE.11.AND.MMAT.LE.30.AND.KQP.NE.0) THEN | |
167 | NQ=NQ+1 | |
168 | KFLO(NQ)=KFP | |
169 | MSTJ(93)=2 | |
170 | PSQ=PSQ+ULMASS_HIJING(KFLO(NQ)) | |
171 | ELSEIF(MMAT.GE.42.AND.MMAT.LE.43.AND.NP.EQ.3.AND.MOD(NQ,2).EQ.1) | |
172 | &THEN | |
173 | NQ=NQ-1 | |
174 | PS=PS-P(I,5) | |
175 | K(I,1)=1 | |
176 | KFI=K(I,2) | |
177 | CALL LUKFDI_HIJING(KFP,KFI,KFLDMP,K(I,2)) | |
178 | IF(K(I,2).EQ.0) GOTO 150 | |
179 | MSTJ(93)=1 | |
180 | P(I,5)=ULMASS_HIJING(K(I,2)) | |
181 | PS=PS+P(I,5) | |
182 | ELSE | |
183 | I=I+1 | |
184 | NP=NP+1 | |
185 | IF(MMAT.NE.33.AND.KQP.NE.0) NQ=NQ+1 | |
186 | IF(MMAT.EQ.33.AND.KQP.NE.0.AND.KQP.NE.2) NQ=NQ+1 | |
187 | K(I,1)=1+MOD(NQ,2) | |
188 | IF(MMAT.EQ.4.AND.JT.LE.2.AND.KFP.EQ.21) K(I,1)=2 | |
189 | IF(MMAT.EQ.4.AND.JT.EQ.3) K(I,1)=1 | |
190 | K(I,2)=KFP | |
191 | K(I,3)=IP | |
192 | K(I,4)=0 | |
193 | K(I,5)=0 | |
194 | P(I,5)=ULMASS_HIJING(KFP) | |
195 | IF(MMAT.EQ.45.AND.KFPA.EQ.89) P(I,5)=PARJ(32) | |
196 | PS=PS+P(I,5) | |
197 | ENDIF | |
198 | 170 CONTINUE | |
199 | ||
200 | C...Choose decay multiplicity in phase space model. | |
201 | 180 IF(MMAT.GE.11.AND.MMAT.LE.30) THEN | |
202 | PSP=PS | |
203 | CNDE=PARJ(61)*LOG(MAX((PV(1,5)-PS-PSQ)/PARJ(62),1.1)) | |
204 | IF(MMAT.EQ.12) CNDE=CNDE+PARJ(63) | |
205 | 190 NTRY=NTRY+1 | |
206 | IF(NTRY.GT.1000) THEN | |
207 | CALL LUERRM_HIJING(14 | |
208 | $ ,'(LUDECY_HIJING:) caught in infinite loop') | |
209 | IF(MSTU(21).GE.1) RETURN | |
210 | ENDIF | |
211 | IF(MMAT.LE.20) THEN | |
212 | GAUSS=SQRT(-2.*CNDE*LOG(MAX(1E-10,RLU_HIJING(0))))* | |
213 | & SIN(PARU(2)*RLU_HIJING(0)) | |
214 | ND=0.5+0.5*NP+0.25*NQ+CNDE+GAUSS | |
215 | IF(ND.LT.NP+NQ/2.OR.ND.LT.2.OR.ND.GT.10) GOTO 190 | |
216 | IF(MMAT.EQ.13.AND.ND.EQ.2) GOTO 190 | |
217 | IF(MMAT.EQ.14.AND.ND.LE.3) GOTO 190 | |
218 | IF(MMAT.EQ.15.AND.ND.LE.4) GOTO 190 | |
219 | ELSE | |
220 | ND=MMAT-20 | |
221 | ENDIF | |
222 | ||
223 | C...Form hadrons from flavour content. | |
224 | DO 200 JT=1,4 | |
225 | 200 KFL1(JT)=KFLO(JT) | |
226 | IF(ND.EQ.NP+NQ/2) GOTO 220 | |
227 | DO 210 I=N+NP+1,N+ND-NQ/2 | |
228 | JT=1+INT((NQ-1)*RLU_HIJING(0)) | |
229 | CALL LUKFDI_HIJING(KFL1(JT),0,KFL2,K(I,2)) | |
230 | IF(K(I,2).EQ.0) GOTO 190 | |
231 | 210 KFL1(JT)=-KFL2 | |
232 | 220 JT=2 | |
233 | JT2=3 | |
234 | JT3=4 | |
235 | IF(NQ.EQ.4.AND.RLU_HIJING(0).LT.PARJ(66)) JT=4 | |
236 | IF(JT.EQ.4.AND.ISIGN(1,KFL1(1)*(10-IABS(KFL1(1))))* | |
237 | & ISIGN(1,KFL1(JT)*(10-IABS(KFL1(JT)))).GT.0) JT=3 | |
238 | IF(JT.EQ.3) JT2=2 | |
239 | IF(JT.EQ.4) JT3=2 | |
240 | CALL LUKFDI_HIJING(KFL1(1),KFL1(JT),KFLDMP,K(N+ND-NQ/2+1,2)) | |
241 | IF(K(N+ND-NQ/2+1,2).EQ.0) GOTO 190 | |
242 | IF(NQ.EQ.4) CALL LUKFDI_HIJING(KFL1(JT2),KFL1(JT3),KFLDMP,K(N+ND | |
243 | $ ,2)) | |
244 | IF(NQ.EQ.4.AND.K(N+ND,2).EQ.0) GOTO 190 | |
245 | ||
246 | C...Check that sum of decay product masses not too large. | |
247 | PS=PSP | |
248 | DO 230 I=N+NP+1,N+ND | |
249 | K(I,1)=1 | |
250 | K(I,3)=IP | |
251 | K(I,4)=0 | |
252 | K(I,5)=0 | |
253 | P(I,5)=ULMASS_HIJING(K(I,2)) | |
254 | 230 PS=PS+P(I,5) | |
255 | IF(PS+PARJ(64).GT.PV(1,5)) GOTO 190 | |
256 | ||
257 | C...Rescale energy to subtract off spectator quark mass. | |
258 | ELSEIF((MMAT.EQ.31.OR.MMAT.EQ.33.OR.MMAT.EQ.44.OR.MMAT.EQ.45). | |
259 | &AND.NP.GE.3) THEN | |
260 | PS=PS-P(N+NP,5) | |
261 | PQT=(P(N+NP,5)+PARJ(65))/PV(1,5) | |
262 | DO 240 J=1,5 | |
263 | P(N+NP,J)=PQT*PV(1,J) | |
264 | 240 PV(1,J)=(1.-PQT)*PV(1,J) | |
265 | IF(PS+PARJ(64).GT.PV(1,5)) GOTO 150 | |
266 | ND=NP-1 | |
267 | MREM=1 | |
268 | ||
269 | C...Phase space factors imposed in W decay. | |
270 | ELSEIF(MMAT.EQ.46) THEN | |
271 | MSTJ(93)=1 | |
272 | PSMC=ULMASS_HIJING(K(N+1,2)) | |
273 | MSTJ(93)=1 | |
274 | PSMC=PSMC+ULMASS_HIJING(K(N+2,2)) | |
275 | IF(MAX(PS,PSMC)+PARJ(32).GT.PV(1,5)) GOTO 130 | |
276 | HR1=(P(N+1,5)/PV(1,5))**2 | |
277 | HR2=(P(N+2,5)/PV(1,5))**2 | |
278 | IF((1.-HR1-HR2)*(2.+HR1+HR2)*SQRT((1.-HR1-HR2)**2-4.*HR1*HR2). | |
279 | & LT.2.*RLU_HIJING(0)) GOTO 130 | |
280 | ND=NP | |
281 | ||
282 | C...Fully specified final state: check mass broadening effects. | |
283 | ELSE | |
284 | IF(NP.GE.2.AND.PS+PARJ(64).GT.PV(1,5)) GOTO 150 | |
285 | ND=NP | |
286 | ENDIF | |
287 | ||
288 | C...Select W mass in decay Q -> W + q, without W propagator. | |
289 | IF(MMAT.EQ.45.AND.MSTJ(25).LE.0) THEN | |
290 | HLQ=(PARJ(32)/PV(1,5))**2 | |
291 | HUQ=(1.-(P(N+2,5)+PARJ(64))/PV(1,5))**2 | |
292 | HRQ=(P(N+2,5)/PV(1,5))**2 | |
293 | 250 HW=HLQ+RLU_HIJING(0)*(HUQ-HLQ) | |
294 | IF(HMEPS(HW).LT.RLU_HIJING(0)) GOTO 250 | |
295 | P(N+1,5)=PV(1,5)*SQRT(HW) | |
296 | ||
297 | C...Ditto, including W propagator. Divide mass range into three regions. | |
298 | ELSEIF(MMAT.EQ.45) THEN | |
299 | HQW=(PV(1,5)/PMAS(24,1))**2 | |
300 | HLW=(PARJ(32)/PMAS(24,1))**2 | |
301 | HUW=((PV(1,5)-P(N+2,5)-PARJ(64))/PMAS(24,1))**2 | |
302 | HRQ=(P(N+2,5)/PV(1,5))**2 | |
303 | HG=PMAS(24,2)/PMAS(24,1) | |
304 | HATL=ATAN((HLW-1.)/HG) | |
305 | HM=MIN(1.,HUW-0.001) | |
306 | HMV1=HMEPS(HM/HQW)/((HM-1.)**2+HG**2) | |
307 | 260 HM=HM-HG | |
308 | HMV2=HMEPS(HM/HQW)/((HM-1.)**2+HG**2) | |
309 | HSAV1=HMEPS(HM/HQW) | |
310 | HSAV2=1./((HM-1.)**2+HG**2) | |
311 | IF(HMV2.GT.HMV1.AND.HM-HG.GT.HLW) THEN | |
312 | HMV1=HMV2 | |
313 | GOTO 260 | |
314 | ENDIF | |
315 | HMV=MIN(2.*HMV1,HMEPS(HM/HQW)/HG**2) | |
316 | HM1=1.-SQRT(1./HMV-HG**2) | |
317 | IF(HM1.GT.HLW.AND.HM1.LT.HM) THEN | |
318 | HM=HM1 | |
319 | ELSEIF(HMV2.LE.HMV1) THEN | |
320 | HM=MAX(HLW,HM-MIN(0.1,1.-HM)) | |
321 | ENDIF | |
322 | HATM=ATAN((HM-1.)/HG) | |
323 | HWT1=(HATM-HATL)/HG | |
324 | HWT2=HMV*(MIN(1.,HUW)-HM) | |
325 | HWT3=0. | |
326 | IF(HUW.GT.1.) THEN | |
327 | HATU=ATAN((HUW-1.)/HG) | |
328 | HMP1=HMEPS(1./HQW) | |
329 | HWT3=HMP1*HATU/HG | |
330 | ENDIF | |
331 | ||
332 | C...Select mass region and W mass there. Accept according to weight. | |
333 | 270 HREG=RLU_HIJING(0)*(HWT1+HWT2+HWT3) | |
334 | IF(HREG.LE.HWT1) THEN | |
335 | HW=1.+HG*TAN(HATL+RLU_HIJING(0)*(HATM-HATL)) | |
336 | HACC=HMEPS(HW/HQW) | |
337 | ELSEIF(HREG.LE.HWT1+HWT2) THEN | |
338 | HW=HM+RLU_HIJING(0)*(MIN(1.,HUW)-HM) | |
339 | HACC=HMEPS(HW/HQW)/((HW-1.)**2+HG**2)/HMV | |
340 | ELSE | |
341 | HW=1.+HG*TAN(RLU_HIJING(0)*HATU) | |
342 | HACC=HMEPS(HW/HQW)/HMP1 | |
343 | ENDIF | |
344 | IF(HACC.LT.RLU_HIJING(0)) GOTO 270 | |
345 | P(N+1,5)=PMAS(24,1)*SQRT(HW) | |
346 | ENDIF | |
347 | ||
348 | C...Determine position of grandmother, number of sisters, Q -> W sign. | |
349 | NM=0 | |
350 | MSGN=0 | |
351 | IF(MMAT.EQ.3.OR.MMAT.EQ.46) THEN | |
352 | IM=K(IP,3) | |
353 | IF(IM.LT.0.OR.IM.GE.IP) IM=0 | |
354 | IF(IM.NE.0) KFAM=IABS(K(IM,2)) | |
355 | IF(IM.NE.0.AND.MMAT.EQ.3) THEN | |
356 | DO 280 IL=MAX(IP-2,IM+1),MIN(IP+2,N) | |
357 | 280 IF(K(IL,3).EQ.IM) NM=NM+1 | |
358 | IF(NM.NE.2.OR.KFAM.LE.100.OR.MOD(KFAM,10).NE.1.OR. | |
359 | & MOD(KFAM/1000,10).NE.0) NM=0 | |
360 | ELSEIF(IM.NE.0.AND.MMAT.EQ.46) THEN | |
361 | MSGN=ISIGN(1,K(IM,2)*K(IP,2)) | |
362 | IF(KFAM.GT.100.AND.MOD(KFAM/1000,10).EQ.0) MSGN= | |
363 | & MSGN*(-1)**MOD(KFAM/100,10) | |
364 | ENDIF | |
365 | ENDIF | |
366 | ||
367 | C...Kinematics of one-particle decays. | |
368 | IF(ND.EQ.1) THEN | |
369 | DO 290 J=1,4 | |
370 | 290 P(N+1,J)=P(IP,J) | |
371 | GOTO 510 | |
372 | ENDIF | |
373 | ||
374 | C...Calculate maximum weight ND-particle decay. | |
375 | PV(ND,5)=P(N+ND,5) | |
376 | IF(ND.GE.3) THEN | |
377 | WTMAX=1./WTCOR(ND-2) | |
378 | PMAX=PV(1,5)-PS+P(N+ND,5) | |
379 | PMIN=0. | |
380 | DO 300 IL=ND-1,1,-1 | |
381 | PMAX=PMAX+P(N+IL,5) | |
382 | PMIN=PMIN+P(N+IL+1,5) | |
383 | 300 WTMAX=WTMAX*PAWT(PMAX,PMIN,P(N+IL,5)) | |
384 | ENDIF | |
385 | ||
386 | C...Find virtual gamma mass in Dalitz decay. | |
387 | 310 IF(ND.EQ.2) THEN | |
388 | ELSEIF(MMAT.EQ.2) THEN | |
389 | PMES=4.*PMAS(11,1)**2 | |
390 | PMRHO2=PMAS(131,1)**2 | |
391 | PGRHO2=PMAS(131,2)**2 | |
392 | 320 PMST=PMES*(P(IP,5)**2/PMES)**RLU_HIJING(0) | |
393 | WT=(1+0.5*PMES/PMST)*SQRT(MAX(0.,1.-PMES/PMST))* | |
394 | & (1.-PMST/P(IP,5)**2)**3*(1.+PGRHO2/PMRHO2)/ | |
395 | & ((1.-PMST/PMRHO2)**2+PGRHO2/PMRHO2) | |
396 | IF(WT.LT.RLU_HIJING(0)) GOTO 320 | |
397 | PV(2,5)=MAX(2.00001*PMAS(11,1),SQRT(PMST)) | |
398 | ||
399 | C...M-generator gives weight. If rejected, try again. | |
400 | ELSE | |
401 | 330 RORD(1)=1. | |
402 | DO 350 IL1=2,ND-1 | |
403 | RSAV=RLU_HIJING(0) | |
404 | DO 340 IL2=IL1-1,1,-1 | |
405 | IF(RSAV.LE.RORD(IL2)) GOTO 350 | |
406 | 340 RORD(IL2+1)=RORD(IL2) | |
407 | 350 RORD(IL2+1)=RSAV | |
408 | RORD(ND)=0. | |
409 | WT=1. | |
410 | DO 360 IL=ND-1,1,-1 | |
411 | PV(IL,5)=PV(IL+1,5)+P(N+IL,5)+(RORD(IL)-RORD(IL+1))*(PV(1,5)-PS) | |
412 | 360 WT=WT*PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5)) | |
413 | IF(WT.LT.RLU_HIJING(0)*WTMAX) GOTO 330 | |
414 | ENDIF | |
415 | ||
416 | C...Perform two-particle decays in respective CM frame. | |
417 | 370 DO 390 IL=1,ND-1 | |
418 | PA=PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5)) | |
419 | UE(3)=2.*RLU_HIJING(0)-1. | |
420 | PHI=PARU(2)*RLU_HIJING(0) | |
421 | UE(1)=SQRT(1.-UE(3)**2)*COS(PHI) | |
422 | UE(2)=SQRT(1.-UE(3)**2)*SIN(PHI) | |
423 | DO 380 J=1,3 | |
424 | P(N+IL,J)=PA*UE(J) | |
425 | 380 PV(IL+1,J)=-PA*UE(J) | |
426 | P(N+IL,4)=SQRT(PA**2+P(N+IL,5)**2) | |
427 | 390 PV(IL+1,4)=SQRT(PA**2+PV(IL+1,5)**2) | |
428 | ||
429 | C...Lorentz transform decay products to lab frame. | |
430 | DO 400 J=1,4 | |
431 | 400 P(N+ND,J)=PV(ND,J) | |
432 | DO 430 IL=ND-1,1,-1 | |
433 | DO 410 J=1,3 | |
434 | 410 BE(J)=PV(IL,J)/PV(IL,4) | |
435 | GA=PV(IL,4)/PV(IL,5) | |
436 | DO 430 I=N+IL,N+ND | |
437 | BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3) | |
438 | DO 420 J=1,3 | |
439 | 420 P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J) | |
440 | 430 P(I,4)=GA*(P(I,4)+BEP) | |
441 | ||
442 | C...Matrix elements for omega and phi decays. | |
443 | IF(MMAT.EQ.1) THEN | |
444 | WT=(P(N+1,5)*P(N+2,5)*P(N+3,5))**2-(P(N+1,5)*FOUR(N+2,N+3))**2 | |
445 | & -(P(N+2,5)*FOUR(N+1,N+3))**2-(P(N+3,5)*FOUR(N+1,N+2))**2 | |
446 | & +2.*FOUR(N+1,N+2)*FOUR(N+1,N+3)*FOUR(N+2,N+3) | |
447 | IF(MAX(WT*WTCOR(9)/P(IP,5)**6,0.001).LT.RLU_HIJING(0)) GOTO 310 | |
448 | ||
449 | C...Matrix elements for pi0 or eta Dalitz decay to gamma e+ e-. | |
450 | ELSEIF(MMAT.EQ.2) THEN | |
451 | FOUR12=FOUR(N+1,N+2) | |
452 | FOUR13=FOUR(N+1,N+3) | |
453 | FOUR23=0.5*PMST-0.25*PMES | |
454 | WT=(PMST-0.5*PMES)*(FOUR12**2+FOUR13**2)+ | |
455 | & PMES*(FOUR12*FOUR13+FOUR12**2+FOUR13**2) | |
456 | IF(WT.LT.RLU_HIJING(0)*0.25*PMST*(P(IP,5)**2-PMST)**2) GOTO 370 | |
457 | ||
458 | C...Matrix element for S0 -> S1 + V1 -> S1 + S2 + S3 (S scalar, | |
459 | C...V vector), of form cos**2(theta02) in V1 rest frame. | |
460 | ELSEIF(MMAT.EQ.3.AND.NM.EQ.2) THEN | |
461 | IF((P(IP,5)**2*FOUR(IM,N+1)-FOUR(IP,IM)*FOUR(IP,N+1))**2.LE. | |
462 | & RLU_HIJING(0)*(FOUR(IP,IM)**2-(P(IP,5)*P(IM,5))**2) | |
463 | $ *(FOUR(IP,N+1)**2-(P(IP,5)*P(N+1,5))**2)) GOTO 370 | |
464 | ||
465 | C...Matrix element for "onium" -> g + g + g or gamma + g + g. | |
466 | ELSEIF(MMAT.EQ.4) THEN | |
467 | HX1=2.*FOUR(IP,N+1)/P(IP,5)**2 | |
468 | HX2=2.*FOUR(IP,N+2)/P(IP,5)**2 | |
469 | HX3=2.*FOUR(IP,N+3)/P(IP,5)**2 | |
470 | WT=((1.-HX1)/(HX2*HX3))**2+((1.-HX2)/(HX1*HX3))**2+ | |
471 | & ((1.-HX3)/(HX1*HX2))**2 | |
472 | IF(WT.LT.2.*RLU_HIJING(0)) GOTO 310 | |
473 | IF(K(IP+1,2).EQ.22.AND.(1.-HX1)*P(IP,5)**2.LT.4.*PARJ(32)**2) | |
474 | & GOTO 310 | |
475 | ||
476 | C...Effective matrix element for nu spectrum in tau -> nu + hadrons. | |
477 | ELSEIF(MMAT.EQ.41) THEN | |
478 | HX1=2.*FOUR(IP,N+1)/P(IP,5)**2 | |
479 | IF(8.*HX1*(3.-2.*HX1)/9..LT.RLU_HIJING(0)) GOTO 310 | |
480 | ||
481 | C...Matrix elements for weak decays (only semileptonic for c and b) | |
482 | ELSEIF(MMAT.GE.42.AND.MMAT.LE.44.AND.ND.EQ.3) THEN | |
483 | IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+3) | |
484 | IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+3) | |
485 | IF(WT.LT.RLU_HIJING(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 310 | |
486 | ELSEIF(MMAT.GE.42.AND.MMAT.LE.44) THEN | |
487 | DO 440 J=1,4 | |
488 | P(N+NP+1,J)=0. | |
489 | DO 440 IS=N+3,N+NP | |
490 | 440 P(N+NP+1,J)=P(N+NP+1,J)+P(IS,J) | |
491 | IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+NP+1) | |
492 | IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+NP+1) | |
493 | IF(WT.LT.RLU_HIJING(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 310 | |
494 | ||
495 | C...Angular distribution in W decay. | |
496 | ELSEIF(MMAT.EQ.46.AND.MSGN.NE.0) THEN | |
497 | IF(MSGN.GT.0) WT=FOUR(IM,N+1)*FOUR(N+2,IP+1) | |
498 | IF(MSGN.LT.0) WT=FOUR(IM,N+2)*FOUR(N+1,IP+1) | |
499 | IF(WT.LT.RLU_HIJING(0)*P(IM,5)**4/WTCOR(10)) GOTO 370 | |
500 | ENDIF | |
501 | ||
502 | C...Scale back energy and reattach spectator. | |
503 | IF(MREM.EQ.1) THEN | |
504 | DO 450 J=1,5 | |
505 | 450 PV(1,J)=PV(1,J)/(1.-PQT) | |
506 | ND=ND+1 | |
507 | MREM=0 | |
508 | ENDIF | |
509 | ||
510 | C...Low invariant mass for system with spectator quark gives particle, | |
511 | C...not two jets. Readjust momenta accordingly. | |
512 | IF((MMAT.EQ.31.OR.MMAT.EQ.45).AND.ND.EQ.3) THEN | |
513 | MSTJ(93)=1 | |
514 | PM2=ULMASS_HIJING(K(N+2,2)) | |
515 | MSTJ(93)=1 | |
516 | PM3=ULMASS_HIJING(K(N+3,2)) | |
517 | IF(P(N+2,5)**2+P(N+3,5)**2+2.*FOUR(N+2,N+3).GE. | |
518 | & (PARJ(32)+PM2+PM3)**2) GOTO 510 | |
519 | K(N+2,1)=1 | |
520 | KFTEMP=K(N+2,2) | |
521 | CALL LUKFDI_HIJING(KFTEMP,K(N+3,2),KFLDMP,K(N+2,2)) | |
522 | IF(K(N+2,2).EQ.0) GOTO 150 | |
523 | P(N+2,5)=ULMASS_HIJING(K(N+2,2)) | |
524 | PS=P(N+1,5)+P(N+2,5) | |
525 | PV(2,5)=P(N+2,5) | |
526 | MMAT=0 | |
527 | ND=2 | |
528 | GOTO 370 | |
529 | ELSEIF(MMAT.EQ.44) THEN | |
530 | MSTJ(93)=1 | |
531 | PM3=ULMASS_HIJING(K(N+3,2)) | |
532 | MSTJ(93)=1 | |
533 | PM4=ULMASS_HIJING(K(N+4,2)) | |
534 | IF(P(N+3,5)**2+P(N+4,5)**2+2.*FOUR(N+3,N+4).GE. | |
535 | & (PARJ(32)+PM3+PM4)**2) GOTO 480 | |
536 | K(N+3,1)=1 | |
537 | KFTEMP=K(N+3,2) | |
538 | CALL LUKFDI_HIJING(KFTEMP,K(N+4,2),KFLDMP,K(N+3,2)) | |
539 | IF(K(N+3,2).EQ.0) GOTO 150 | |
540 | P(N+3,5)=ULMASS_HIJING(K(N+3,2)) | |
541 | DO 460 J=1,3 | |
542 | 460 P(N+3,J)=P(N+3,J)+P(N+4,J) | |
543 | P(N+3,4)=SQRT(P(N+3,1)**2+P(N+3,2)**2+P(N+3,3)**2+P(N+3,5)**2) | |
544 | HA=P(N+1,4)**2-P(N+2,4)**2 | |
545 | HB=HA-(P(N+1,5)**2-P(N+2,5)**2) | |
546 | HC=(P(N+1,1)-P(N+2,1))**2+(P(N+1,2)-P(N+2,2))**2+ | |
547 | & (P(N+1,3)-P(N+2,3))**2 | |
548 | HD=(PV(1,4)-P(N+3,4))**2 | |
549 | HE=HA**2-2.*HD*(P(N+1,4)**2+P(N+2,4)**2)+HD**2 | |
550 | HF=HD*HC-HB**2 | |
551 | HG=HD*HC-HA*HB | |
552 | HH=(SQRT(HG**2+HE*HF)-HG)/(2.*HF) | |
553 | DO 470 J=1,3 | |
554 | PCOR=HH*(P(N+1,J)-P(N+2,J)) | |
555 | P(N+1,J)=P(N+1,J)+PCOR | |
556 | 470 P(N+2,J)=P(N+2,J)-PCOR | |
557 | P(N+1,4)=SQRT(P(N+1,1)**2+P(N+1,2)**2+P(N+1,3)**2+P(N+1,5)**2) | |
558 | P(N+2,4)=SQRT(P(N+2,1)**2+P(N+2,2)**2+P(N+2,3)**2+P(N+2,5)**2) | |
559 | ND=ND-1 | |
560 | ENDIF | |
561 | ||
562 | C...Check invariant mass of W jets. May give one particle or start over. | |
563 | 480 IF(MMAT.GE.42.AND.MMAT.LE.44.AND.IABS(K(N+1,2)).LT.10) THEN | |
564 | PMR=SQRT(MAX(0.,P(N+1,5)**2+P(N+2,5)**2+2.*FOUR(N+1,N+2))) | |
565 | MSTJ(93)=1 | |
566 | PM1=ULMASS_HIJING(K(N+1,2)) | |
567 | MSTJ(93)=1 | |
568 | PM2=ULMASS_HIJING(K(N+2,2)) | |
569 | IF(PMR.GT.PARJ(32)+PM1+PM2) GOTO 490 | |
570 | KFLDUM=INT(1.5+RLU_HIJING(0)) | |
571 | CALL LUKFDI_HIJING(K(N+1,2),-ISIGN(KFLDUM,K(N+1,2)),KFLDMP,KF1) | |
572 | CALL LUKFDI_HIJING(K(N+2,2),-ISIGN(KFLDUM,K(N+2,2)),KFLDMP,KF2) | |
573 | IF(KF1.EQ.0.OR.KF2.EQ.0) GOTO 150 | |
574 | PSM=ULMASS_HIJING(KF1)+ULMASS_HIJING(KF2) | |
575 | IF(MMAT.EQ.42.AND.PMR.GT.PARJ(64)+PSM) GOTO 490 | |
576 | IF(MMAT.GE.43.AND.PMR.GT.0.2*PARJ(32)+PSM) GOTO 490 | |
577 | IF(ND.EQ.4.OR.KFA.EQ.15) GOTO 150 | |
578 | K(N+1,1)=1 | |
579 | KFTEMP=K(N+1,2) | |
580 | CALL LUKFDI_HIJING(KFTEMP,K(N+2,2),KFLDMP,K(N+1,2)) | |
581 | IF(K(N+1,2).EQ.0) GOTO 150 | |
582 | P(N+1,5)=ULMASS_HIJING(K(N+1,2)) | |
583 | K(N+2,2)=K(N+3,2) | |
584 | P(N+2,5)=P(N+3,5) | |
585 | PS=P(N+1,5)+P(N+2,5) | |
586 | PV(2,5)=P(N+3,5) | |
587 | MMAT=0 | |
588 | ND=2 | |
589 | GOTO 370 | |
590 | ENDIF | |
591 | ||
592 | C...Phase space decay of partons from W decay. | |
593 | 490 IF(MMAT.EQ.42.AND.IABS(K(N+1,2)).LT.10) THEN | |
594 | KFLO(1)=K(N+1,2) | |
595 | KFLO(2)=K(N+2,2) | |
596 | K(N+1,1)=K(N+3,1) | |
597 | K(N+1,2)=K(N+3,2) | |
598 | DO 500 J=1,5 | |
599 | PV(1,J)=P(N+1,J)+P(N+2,J) | |
600 | 500 P(N+1,J)=P(N+3,J) | |
601 | PV(1,5)=PMR | |
602 | N=N+1 | |
603 | NP=0 | |
604 | NQ=2 | |
605 | PS=0. | |
606 | MSTJ(93)=2 | |
607 | PSQ=ULMASS_HIJING(KFLO(1)) | |
608 | MSTJ(93)=2 | |
609 | PSQ=PSQ+ULMASS_HIJING(KFLO(2)) | |
610 | MMAT=11 | |
611 | GOTO 180 | |
612 | ENDIF | |
613 | ||
614 | C...Boost back for rapidly moving particle. | |
615 | 510 N=N+ND | |
616 | IF(MBST.EQ.1) THEN | |
617 | DO 520 J=1,3 | |
618 | 520 BE(J)=P(IP,J)/P(IP,4) | |
619 | GA=P(IP,4)/P(IP,5) | |
620 | DO 540 I=NSAV+1,N | |
621 | BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3) | |
622 | DO 530 J=1,3 | |
623 | 530 P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J) | |
624 | 540 P(I,4)=GA*(P(I,4)+BEP) | |
625 | ENDIF | |
626 | ||
627 | C...Fill in position of decay vertex. | |
628 | DO 560 I=NSAV+1,N | |
629 | DO 550 J=1,4 | |
630 | 550 V(I,J)=VDCY(J) | |
631 | 560 V(I,5)=0. | |
632 | ||
633 | C...Set up for parton shower evolution from jets. | |
634 | IF(MSTJ(23).GE.1.AND.MMAT.EQ.4.AND.K(NSAV+1,2).EQ.21) THEN | |
635 | K(NSAV+1,1)=3 | |
636 | K(NSAV+2,1)=3 | |
637 | K(NSAV+3,1)=3 | |
638 | K(NSAV+1,4)=MSTU(5)*(NSAV+2) | |
639 | K(NSAV+1,5)=MSTU(5)*(NSAV+3) | |
640 | K(NSAV+2,4)=MSTU(5)*(NSAV+3) | |
641 | K(NSAV+2,5)=MSTU(5)*(NSAV+1) | |
642 | K(NSAV+3,4)=MSTU(5)*(NSAV+1) | |
643 | K(NSAV+3,5)=MSTU(5)*(NSAV+2) | |
644 | MSTJ(92)=-(NSAV+1) | |
645 | ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.4) THEN | |
646 | K(NSAV+2,1)=3 | |
647 | K(NSAV+3,1)=3 | |
648 | K(NSAV+2,4)=MSTU(5)*(NSAV+3) | |
649 | K(NSAV+2,5)=MSTU(5)*(NSAV+3) | |
650 | K(NSAV+3,4)=MSTU(5)*(NSAV+2) | |
651 | K(NSAV+3,5)=MSTU(5)*(NSAV+2) | |
652 | MSTJ(92)=NSAV+2 | |
653 | ELSEIF(MSTJ(23).GE.1.AND.(MMAT.EQ.32.OR.MMAT.EQ.44.OR.MMAT.EQ.46). | |
654 | &AND.IABS(K(NSAV+1,2)).LE.10.AND.IABS(K(NSAV+2,2)).LE.10) THEN | |
655 | K(NSAV+1,1)=3 | |
656 | K(NSAV+2,1)=3 | |
657 | K(NSAV+1,4)=MSTU(5)*(NSAV+2) | |
658 | K(NSAV+1,5)=MSTU(5)*(NSAV+2) | |
659 | K(NSAV+2,4)=MSTU(5)*(NSAV+1) | |
660 | K(NSAV+2,5)=MSTU(5)*(NSAV+1) | |
661 | MSTJ(92)=NSAV+1 | |
662 | ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33.AND.IABS(K(NSAV+2,2)).EQ.21) | |
663 | &THEN | |
664 | K(NSAV+1,1)=3 | |
665 | K(NSAV+2,1)=3 | |
666 | K(NSAV+3,1)=3 | |
667 | KCP=LUCOMP_HIJING(K(NSAV+1,2)) | |
668 | KQP=KCHG(KCP,2)*ISIGN(1,K(NSAV+1,2)) | |
669 | JCON=4 | |
670 | IF(KQP.LT.0) JCON=5 | |
671 | K(NSAV+1,JCON)=MSTU(5)*(NSAV+2) | |
672 | K(NSAV+2,9-JCON)=MSTU(5)*(NSAV+1) | |
673 | K(NSAV+2,JCON)=MSTU(5)*(NSAV+3) | |
674 | K(NSAV+3,9-JCON)=MSTU(5)*(NSAV+2) | |
675 | MSTJ(92)=NSAV+1 | |
676 | ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33) THEN | |
677 | K(NSAV+1,1)=3 | |
678 | K(NSAV+3,1)=3 | |
679 | K(NSAV+1,4)=MSTU(5)*(NSAV+3) | |
680 | K(NSAV+1,5)=MSTU(5)*(NSAV+3) | |
681 | K(NSAV+3,4)=MSTU(5)*(NSAV+1) | |
682 | K(NSAV+3,5)=MSTU(5)*(NSAV+1) | |
683 | MSTJ(92)=NSAV+1 | |
684 | ENDIF | |
685 | ||
686 | C...Mark decayed particle. | |
687 | IF(K(IP,1).EQ.5) K(IP,1)=15 | |
688 | IF(K(IP,1).LE.10) K(IP,1)=11 | |
689 | K(IP,4)=NSAV+1 | |
690 | K(IP,5)=N | |
691 | ||
692 | RETURN | |
693 | END |