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