]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA/jetset/ludecy.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Moved GetGoodParticles to alimacros
[u/mrichter/AliRoot.git] / PYTHIA / jetset / ludecy.F
CommitLineData
fe4da5cc 1
2C*********************************************************************
3
4 SUBROUTINE LUDECY(IP)
5
6C...Purpose: to handle the decay of unstable particles.
7 COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)
8 COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
9 COMMON/LUDAT2/KCHG(500,3),PMAS(500,4),PARF(2000),VCKM(4,4)
10 COMMON/LUDAT3/MDCY(500,3),MDME(2000,2),BRAT(2000),KFDP(2000,5)
11 SAVE /LUJETS/,/LUDAT1/,/LUDAT2/,/LUDAT3/
12 DIMENSION VDCY(4),KFLO(4),KFL1(4),PV(10,5),RORD(10),UE(3),BE(3),
13 &WTCOR(10),PTAU(4),PCMTAU(4)
14 DOUBLE PRECISION DBETAU(3)
15 DATA WTCOR/2.,5.,15.,60.,250.,1500.,1.2E4,1.2E5,150.,16./
16
17C...Functions: momentum in two-particle decays, four-product and
18C...matrix element times phase space in weak decays.
19 PAWT(A,B,C)=SQRT((A**2-(B+C)**2)*(A**2-(B-C)**2))/(2.*A)
20 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)
21 HMEPS(HA)=((1.-HRQ-HA)**2+3.*HA*(1.+HRQ-HA))*
22 &SQRT((1.-HRQ-HA)**2-4.*HRQ*HA)
23
24C...Initial values.
25 NTRY=0
26 NSAV=N
27 KFA=IABS(K(IP,2))
28 KFS=ISIGN(1,K(IP,2))
29 KC=LUCOMP(KFA)
30 MSTJ(92)=0
31
32C...Choose lifetime and determine decay vertex.
33 IF(K(IP,1).EQ.5) THEN
34 V(IP,5)=0.
35 ELSEIF(K(IP,1).NE.4) THEN
36 V(IP,5)=-PMAS(KC,4)*LOG(RLU(0))
37 ENDIF
38 DO 100 J=1,4
39 VDCY(J)=V(IP,J)+V(IP,5)*P(IP,J)/P(IP,5)
40 100 CONTINUE
41
42C...Determine whether decay allowed or not.
43 MOUT=0
44 IF(MSTJ(22).EQ.2) THEN
45 IF(PMAS(KC,4).GT.PARJ(71)) MOUT=1
46 ELSEIF(MSTJ(22).EQ.3) THEN
47 IF(VDCY(1)**2+VDCY(2)**2+VDCY(3)**2.GT.PARJ(72)**2) MOUT=1
48 ELSEIF(MSTJ(22).EQ.4) THEN
49 IF(VDCY(1)**2+VDCY(2)**2.GT.PARJ(73)**2) MOUT=1
50 IF(ABS(VDCY(3)).GT.PARJ(74)) MOUT=1
51 ENDIF
52 IF(MOUT.EQ.1.AND.K(IP,1).NE.5) THEN
53 K(IP,1)=4
54 RETURN
55 ENDIF
56
57C...Interface to external tau decay library (for tau polarization).
58 IF(KFA.EQ.15.AND.MSTJ(28).GE.1) THEN
59
60C...Starting values for pointers and momenta.
61 ITAU=IP
62 DO 110 J=1,4
63 PTAU(J)=P(ITAU,J)
64 PCMTAU(J)=P(ITAU,J)
65 110 CONTINUE
66
67C...Iterate to find position and code of mother of tau.
68 IMTAU=ITAU
69 120 IMTAU=K(IMTAU,3)
70
71 IF(IMTAU.EQ.0) THEN
72C...If no known origin then impossible to do anything further.
73 KFORIG=0
74 IORIG=0
75
76 ELSEIF(K(IMTAU,2).EQ.K(ITAU,2)) THEN
77C...If tau -> tau + gamma then add gamma energy and loop.
78 IF(K(K(IMTAU,4),2).EQ.22) THEN
79 DO 130 J=1,4
80 PCMTAU(J)=PCMTAU(J)+P(K(IMTAU,4),J)
81 130 CONTINUE
82 ELSEIF(K(K(IMTAU,5),2).EQ.22) THEN
83 DO 140 J=1,4
84 PCMTAU(J)=PCMTAU(J)+P(K(IMTAU,5),J)
85 140 CONTINUE
86 ENDIF
87 GOTO 120
88
89 ELSEIF(IABS(K(IMTAU,2)).GT.100) THEN
90C...If coming from weak decay of hadron then W is not stored in record,
91C...but can be reconstructed by adding neutrino momentum.
92 KFORIG=-ISIGN(24,K(ITAU,2))
93 IORIG=0
94 DO 160 II=K(IMTAU,4),K(IMTAU,5)
95 IF(K(II,2)*ISIGN(1,K(ITAU,2)).EQ.-16) THEN
96 DO 150 J=1,4
97 PCMTAU(J)=PCMTAU(J)+P(II,J)
98 150 CONTINUE
99 ENDIF
100 160 CONTINUE
101
102 ELSE
103C...If coming from resonance decay then find latest copy of this
104C...resonance (may not completely agree).
105 KFORIG=K(IMTAU,2)
106 IORIG=IMTAU
107 DO 170 II=IMTAU+1,IP-1
108 IF(K(II,2).EQ.KFORIG.AND.K(II,3).EQ.IORIG.AND.
109 & ABS(P(II,5)-P(IORIG,5)).LT.1E-5*P(IORIG,5)) IORIG=II
110 170 CONTINUE
111 DO 180 J=1,4
112 PCMTAU(J)=P(IORIG,J)
113 180 CONTINUE
114 ENDIF
115
116C...Boost tau to rest frame of production process (where known)
117C...and rotate it to sit along +z axis.
118 DO 190 J=1,3
119 DBETAU(J)=PCMTAU(J)/PCMTAU(4)
120 190 CONTINUE
121 IF(KFORIG.NE.0) CALL LUDBRB(ITAU,ITAU,0.,0.,-DBETAU(1),
122 & -DBETAU(2),-DBETAU(3))
123 PHITAU=ULANGL(P(ITAU,1),P(ITAU,2))
124 CALL LUDBRB(ITAU,ITAU,0.,-PHITAU,0D0,0D0,0D0)
125 THETAU=ULANGL(P(ITAU,3),P(ITAU,1))
126 CALL LUDBRB(ITAU,ITAU,-THETAU,0.,0D0,0D0,0D0)
127
128C...Call tau decay routine (if meaningful) and fill extra info.
129 IF(KFORIG.NE.0.OR.MSTJ(28).EQ.2) THEN
130 CALL LUTAUD(ITAU,IORIG,KFORIG,NDECAY)
131 DO 200 II=NSAV+1,NSAV+NDECAY
132 K(II,1)=1
133 K(II,3)=IP
134 K(II,4)=0
135 K(II,5)=0
136 200 CONTINUE
137 N=NSAV+NDECAY
138 ENDIF
139
140C...Boost back decay tau and decay products.
141 DO 210 J=1,4
142 P(ITAU,J)=PTAU(J)
143 210 CONTINUE
144 IF(KFORIG.NE.0.OR.MSTJ(28).EQ.2) THEN
145 CALL LUDBRB(NSAV+1,N,THETAU,PHITAU,0D0,0D0,0D0)
146 IF(KFORIG.NE.0) CALL LUDBRB(NSAV+1,N,0.,0.,DBETAU(1),
147 & DBETAU(2),DBETAU(3))
148
149C...Skip past ordinary tau decay treatment.
150 MMAT=0
151 MBST=0
152 ND=0
153 GOTO 660
154 ENDIF
155 ENDIF
156
157C...B-B~ mixing: flip sign of meson appropriately.
158 MMIX=0
159 IF((KFA.EQ.511.OR.KFA.EQ.531).AND.MSTJ(26).GE.1) THEN
160 XBBMIX=PARJ(76)
161 IF(KFA.EQ.531) XBBMIX=PARJ(77)
162 IF(SIN(0.5*XBBMIX*V(IP,5)/PMAS(KC,4))**2.GT.RLU(0)) MMIX=1
163 IF(MMIX.EQ.1) KFS=-KFS
164 ENDIF
165
166C...Check existence of decay channels. Particle/antiparticle rules.
167 KCA=KC
168 IF(MDCY(KC,2).GT.0) THEN
169 MDMDCY=MDME(MDCY(KC,2),2)
170 IF(MDMDCY.GT.80.AND.MDMDCY.LE.90) KCA=MDMDCY
171 ENDIF
172 IF(MDCY(KCA,2).LE.0.OR.MDCY(KCA,3).LE.0) THEN
173 CALL LUERRM(9,'(LUDECY:) no decay channel defined')
174 RETURN
175 ENDIF
176 IF(MOD(KFA/1000,10).EQ.0.AND.(KCA.EQ.85.OR.KCA.EQ.87)) KFS=-KFS
177 IF(KCHG(KC,3).EQ.0) THEN
178 KFSP=1
179 KFSN=0
180 IF(RLU(0).GT.0.5) KFS=-KFS
181 ELSEIF(KFS.GT.0) THEN
182 KFSP=1
183 KFSN=0
184 ELSE
185 KFSP=0
186 KFSN=1
187 ENDIF
188
189C...Sum branching ratios of allowed decay channels.
190 220 NOPE=0
191 BRSU=0.
192 DO 230 IDL=MDCY(KCA,2),MDCY(KCA,2)+MDCY(KCA,3)-1
193 IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND.
194 &KFSN*MDME(IDL,1).NE.3) GOTO 230
195 IF(MDME(IDL,2).GT.100) GOTO 230
196 NOPE=NOPE+1
197 BRSU=BRSU+BRAT(IDL)
198 230 CONTINUE
199 IF(NOPE.EQ.0) THEN
200 CALL LUERRM(2,'(LUDECY:) all decay channels closed by user')
201 RETURN
202 ENDIF
203
204C...Select decay channel among allowed ones.
205 240 RBR=BRSU*RLU(0)
206 IDL=MDCY(KCA,2)-1
207 250 IDL=IDL+1
208 IF(MDME(IDL,1).NE.1.AND.KFSP*MDME(IDL,1).NE.2.AND.
209 &KFSN*MDME(IDL,1).NE.3) THEN
210 IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 250
211 ELSEIF(MDME(IDL,2).GT.100) THEN
212 IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1) GOTO 250
213 ELSE
214 IDC=IDL
215 RBR=RBR-BRAT(IDL)
216 IF(IDL.LT.MDCY(KCA,2)+MDCY(KCA,3)-1.AND.RBR.GT.0.) GOTO 250
217 ENDIF
218
219C...Start readout of decay channel: matrix element, reset counters.
220 MMAT=MDME(IDC,2)
221 260 NTRY=NTRY+1
222 IF(NTRY.GT.1000) THEN
223 CALL LUERRM(14,'(LUDECY:) caught in infinite loop')
224 IF(MSTU(21).GE.1) RETURN
225 ENDIF
226 I=N
227 NP=0
228 NQ=0
229 MBST=0
230 IF(MMAT.GE.11.AND.MMAT.NE.46.AND.P(IP,4).GT.20.*P(IP,5)) MBST=1
231 DO 270 J=1,4
232 PV(1,J)=0.
233 IF(MBST.EQ.0) PV(1,J)=P(IP,J)
234 270 CONTINUE
235 IF(MBST.EQ.1) PV(1,4)=P(IP,5)
236 PV(1,5)=P(IP,5)
237 PS=0.
238 PSQ=0.
239 MREM=0
240 MHADDY=0
241 IF(KFA.GT.80) MHADDY=1
242
243C...Read out decay products. Convert to standard flavour code.
244 JTMAX=5
245 IF(MDME(IDC+1,2).EQ.101) JTMAX=10
246 DO 280 JT=1,JTMAX
247 IF(JT.LE.5) KP=KFDP(IDC,JT)
248 IF(JT.GE.6) KP=KFDP(IDC+1,JT-5)
249 IF(KP.EQ.0) GOTO 280
250 KPA=IABS(KP)
251 KCP=LUCOMP(KPA)
252 IF(KPA.GT.80) MHADDY=1
253 IF(KCHG(KCP,3).EQ.0.AND.KPA.NE.81.AND.KPA.NE.82) THEN
254 KFP=KP
255 ELSEIF(KPA.NE.81.AND.KPA.NE.82) THEN
256 KFP=KFS*KP
257 ELSEIF(KPA.EQ.81.AND.MOD(KFA/1000,10).EQ.0) THEN
258 KFP=-KFS*MOD(KFA/10,10)
259 ELSEIF(KPA.EQ.81.AND.MOD(KFA/100,10).GE.MOD(KFA/10,10)) THEN
260 KFP=KFS*(100*MOD(KFA/10,100)+3)
261 ELSEIF(KPA.EQ.81) THEN
262 KFP=KFS*(1000*MOD(KFA/10,10)+100*MOD(KFA/100,10)+1)
263 ELSEIF(KP.EQ.82) THEN
264 CALL LUKFDI(-KFS*INT(1.+(2.+PARJ(2))*RLU(0)),0,KFP,KDUMP)
265 IF(KFP.EQ.0) GOTO 260
266 MSTJ(93)=1
267 IF(PV(1,5).LT.PARJ(32)+2.*ULMASS(KFP)) GOTO 260
268 ELSEIF(KP.EQ.-82) THEN
269 KFP=-KFP
270 IF(IABS(KFP).GT.10) KFP=KFP+ISIGN(10000,KFP)
271 ENDIF
272 IF(KPA.EQ.81.OR.KPA.EQ.82) KCP=LUCOMP(KFP)
273
274C...Add decay product to event record or to quark flavour list.
275 KFPA=IABS(KFP)
276 KQP=KCHG(KCP,2)
277 IF(MMAT.GE.11.AND.MMAT.LE.30.AND.KQP.NE.0) THEN
278 NQ=NQ+1
279 KFLO(NQ)=KFP
280 MSTJ(93)=2
281 PSQ=PSQ+ULMASS(KFLO(NQ))
282 ELSEIF((MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.48).AND.NP.EQ.3.AND.
283 &MOD(NQ,2).EQ.1) THEN
284 NQ=NQ-1
285 PS=PS-P(I,5)
286 K(I,1)=1
287 KFI=K(I,2)
288 CALL LUKFDI(KFP,KFI,KFLDMP,K(I,2))
289 IF(K(I,2).EQ.0) GOTO 260
290 MSTJ(93)=1
291 P(I,5)=ULMASS(K(I,2))
292 PS=PS+P(I,5)
293 ELSE
294 I=I+1
295 NP=NP+1
296 IF(MMAT.NE.33.AND.KQP.NE.0) NQ=NQ+1
297 IF(MMAT.EQ.33.AND.KQP.NE.0.AND.KQP.NE.2) NQ=NQ+1
298 K(I,1)=1+MOD(NQ,2)
299 IF(MMAT.EQ.4.AND.JT.LE.2.AND.KFP.EQ.21) K(I,1)=2
300 IF(MMAT.EQ.4.AND.JT.EQ.3) K(I,1)=1
301 K(I,2)=KFP
302 K(I,3)=IP
303 K(I,4)=0
304 K(I,5)=0
305 P(I,5)=ULMASS(KFP)
306 IF(MMAT.EQ.45.AND.KFPA.EQ.89) P(I,5)=PARJ(32)
307 PS=PS+P(I,5)
308 ENDIF
309 280 CONTINUE
310
311C...Check masses for resonance decays.
312 IF(MHADDY.EQ.0) THEN
313 IF(PS+PARJ(64).GT.PV(1,5)) GOTO 240
314 ENDIF
315
316C...Choose decay multiplicity in phase space model.
317 290 IF(MMAT.GE.11.AND.MMAT.LE.30) THEN
318 PSP=PS
319 CNDE=PARJ(61)*LOG(MAX((PV(1,5)-PS-PSQ)/PARJ(62),1.1))
320 IF(MMAT.EQ.12) CNDE=CNDE+PARJ(63)
321 300 NTRY=NTRY+1
322 IF(NTRY.GT.1000) THEN
323 CALL LUERRM(14,'(LUDECY:) caught in infinite loop')
324 IF(MSTU(21).GE.1) RETURN
325 ENDIF
326 IF(MMAT.LE.20) THEN
327 GAUSS=SQRT(-2.*CNDE*LOG(MAX(1E-10,RLU(0))))*
328 & SIN(PARU(2)*RLU(0))
329 ND=0.5+0.5*NP+0.25*NQ+CNDE+GAUSS
330 IF(ND.LT.NP+NQ/2.OR.ND.LT.2.OR.ND.GT.10) GOTO 300
331 IF(MMAT.EQ.13.AND.ND.EQ.2) GOTO 300
332 IF(MMAT.EQ.14.AND.ND.LE.3) GOTO 300
333 IF(MMAT.EQ.15.AND.ND.LE.4) GOTO 300
334 ELSE
335 ND=MMAT-20
336 ENDIF
337
338C...Form hadrons from flavour content.
339 DO 310 JT=1,4
340 KFL1(JT)=KFLO(JT)
341 310 CONTINUE
342 IF(ND.EQ.NP+NQ/2) GOTO 330
343 DO 320 I=N+NP+1,N+ND-NQ/2
344 JT=1+INT((NQ-1)*RLU(0))
345 CALL LUKFDI(KFL1(JT),0,KFL2,K(I,2))
346 IF(K(I,2).EQ.0) GOTO 300
347 KFL1(JT)=-KFL2
348 320 CONTINUE
349 330 JT=2
350 JT2=3
351 JT3=4
352 IF(NQ.EQ.4.AND.RLU(0).LT.PARJ(66)) JT=4
353 IF(JT.EQ.4.AND.ISIGN(1,KFL1(1)*(10-IABS(KFL1(1))))*
354 & ISIGN(1,KFL1(JT)*(10-IABS(KFL1(JT)))).GT.0) JT=3
355 IF(JT.EQ.3) JT2=2
356 IF(JT.EQ.4) JT3=2
357 CALL LUKFDI(KFL1(1),KFL1(JT),KFLDMP,K(N+ND-NQ/2+1,2))
358 IF(K(N+ND-NQ/2+1,2).EQ.0) GOTO 300
359 IF(NQ.EQ.4) CALL LUKFDI(KFL1(JT2),KFL1(JT3),KFLDMP,K(N+ND,2))
360 IF(NQ.EQ.4.AND.K(N+ND,2).EQ.0) GOTO 300
361
362C...Check that sum of decay product masses not too large.
363 PS=PSP
364 DO 340 I=N+NP+1,N+ND
365 K(I,1)=1
366 K(I,3)=IP
367 K(I,4)=0
368 K(I,5)=0
369 P(I,5)=ULMASS(K(I,2))
370 PS=PS+P(I,5)
371 340 CONTINUE
372 IF(PS+PARJ(64).GT.PV(1,5)) GOTO 300
373
374C...Rescale energy to subtract off spectator quark mass.
375 ELSEIF((MMAT.EQ.31.OR.MMAT.EQ.33.OR.MMAT.EQ.44.OR.MMAT.EQ.45)
376 &.AND.NP.GE.3) THEN
377 PS=PS-P(N+NP,5)
378 PQT=(P(N+NP,5)+PARJ(65))/PV(1,5)
379 DO 350 J=1,5
380 P(N+NP,J)=PQT*PV(1,J)
381 PV(1,J)=(1.-PQT)*PV(1,J)
382 350 CONTINUE
383 IF(PS+PARJ(64).GT.PV(1,5)) GOTO 260
384 ND=NP-1
385 MREM=1
386
387C...Phase space factors imposed in W decay.
388 ELSEIF(MMAT.EQ.46) THEN
389 MSTJ(93)=1
390 PSMC=ULMASS(K(N+1,2))
391 MSTJ(93)=1
392 PSMC=PSMC+ULMASS(K(N+2,2))
393 IF(MAX(PS,PSMC)+PARJ(32).GT.PV(1,5)) GOTO 240
394 HR1=(P(N+1,5)/PV(1,5))**2
395 HR2=(P(N+2,5)/PV(1,5))**2
396 IF((1.-HR1-HR2)*(2.+HR1+HR2)*SQRT((1.-HR1-HR2)**2-4.*HR1*HR2)
397 & .LT.2.*RLU(0)) GOTO 240
398 ND=NP
399
400C...Fully specified final state: check mass broadening effects.
401 ELSE
402 IF(NP.GE.2.AND.PS+PARJ(64).GT.PV(1,5)) GOTO 260
403 ND=NP
404 ENDIF
405
406C...Select W mass in decay Q -> W + q, without W propagator.
407 IF(MMAT.EQ.45.AND.MSTJ(25).LE.0) THEN
408 HLQ=(PARJ(32)/PV(1,5))**2
409 HUQ=(1.-(P(N+2,5)+PARJ(64))/PV(1,5))**2
410 HRQ=(P(N+2,5)/PV(1,5))**2
411 360 HW=HLQ+RLU(0)*(HUQ-HLQ)
412 IF(HMEPS(HW).LT.RLU(0)) GOTO 360
413 P(N+1,5)=PV(1,5)*SQRT(HW)
414
415C...Ditto, including W propagator. Divide mass range into three regions.
416 ELSEIF(MMAT.EQ.45) THEN
417 HQW=(PV(1,5)/PMAS(24,1))**2
418 HLW=(PARJ(32)/PMAS(24,1))**2
419 HUW=((PV(1,5)-P(N+2,5)-PARJ(64))/PMAS(24,1))**2
420 HRQ=(P(N+2,5)/PV(1,5))**2
421 HG=PMAS(24,2)/PMAS(24,1)
422 HATL=ATAN((HLW-1.)/HG)
423 HM=MIN(1.,HUW-0.001)
424 HMV1=HMEPS(HM/HQW)/((HM-1.)**2+HG**2)
425 370 HM=HM-HG
426 HMV2=HMEPS(HM/HQW)/((HM-1.)**2+HG**2)
427 IF(HMV2.GT.HMV1.AND.HM-HG.GT.HLW) THEN
428 HMV1=HMV2
429 GOTO 370
430 ENDIF
431 HMV=MIN(2.*HMV1,HMEPS(HM/HQW)/HG**2)
432 HM1=1.-SQRT(1./HMV-HG**2)
433 IF(HM1.GT.HLW.AND.HM1.LT.HM) THEN
434 HM=HM1
435 ELSEIF(HMV2.LE.HMV1) THEN
436 HM=MAX(HLW,HM-MIN(0.1,1.-HM))
437 ENDIF
438 HATM=ATAN((HM-1.)/HG)
439 HWT1=(HATM-HATL)/HG
440 HWT2=HMV*(MIN(1.,HUW)-HM)
441 HWT3=0.
442 IF(HUW.GT.1.) THEN
443 HATU=ATAN((HUW-1.)/HG)
444 HMP1=HMEPS(1./HQW)
445 HWT3=HMP1*HATU/HG
446 ENDIF
447
448C...Select mass region and W mass there. Accept according to weight.
449 380 HREG=RLU(0)*(HWT1+HWT2+HWT3)
450 IF(HREG.LE.HWT1) THEN
451 HW=1.+HG*TAN(HATL+RLU(0)*(HATM-HATL))
452 HACC=HMEPS(HW/HQW)
453 ELSEIF(HREG.LE.HWT1+HWT2) THEN
454 HW=HM+RLU(0)*(MIN(1.,HUW)-HM)
455 HACC=HMEPS(HW/HQW)/((HW-1.)**2+HG**2)/HMV
456 ELSE
457 HW=1.+HG*TAN(RLU(0)*HATU)
458 HACC=HMEPS(HW/HQW)/HMP1
459 ENDIF
460 IF(HACC.LT.RLU(0)) GOTO 380
461 P(N+1,5)=PMAS(24,1)*SQRT(HW)
462 ENDIF
463
464C...Determine position of grandmother, number of sisters, Q -> W sign.
465 NM=0
466 KFAS=0
467 MSGN=0
468 IF(MMAT.EQ.3.OR.MMAT.EQ.46) THEN
469 IM=K(IP,3)
470 IF(IM.LT.0.OR.IM.GE.IP) IM=0
471 IF(MMAT.EQ.46.AND.MSTJ(27).EQ.1) THEN
472 IM=0
473 ELSEIF(MMAT.EQ.46.AND.MSTJ(27).GE.2.AND.IM.NE.0) THEN
474 IF(K(IM,2).EQ.94) THEN
475 IM=K(K(IM,3),3)
476 IF(IM.LT.0.OR.IM.GE.IP) IM=0
477 ENDIF
478 ENDIF
479 IF(IM.NE.0) KFAM=IABS(K(IM,2))
480 IF(IM.NE.0.AND.MMAT.EQ.3) THEN
481 DO 390 IL=MAX(IP-2,IM+1),MIN(IP+2,N)
482 IF(K(IL,3).EQ.IM) NM=NM+1
483 IF(K(IL,3).EQ.IM.AND.IL.NE.IP) ISIS=IL
484 390 CONTINUE
485 IF(NM.NE.2.OR.KFAM.LE.100.OR.MOD(KFAM,10).NE.1.OR.
486 & MOD(KFAM/1000,10).NE.0) NM=0
487 IF(NM.EQ.2) THEN
488 KFAS=IABS(K(ISIS,2))
489 IF((KFAS.LE.100.OR.MOD(KFAS,10).NE.1.OR.
490 & MOD(KFAS/1000,10).NE.0).AND.KFAS.NE.22) NM=0
491 ENDIF
492 ELSEIF(IM.NE.0.AND.MMAT.EQ.46) THEN
493 MSGN=ISIGN(1,K(IM,2)*K(IP,2))
494 IF(KFAM.GT.100.AND.MOD(KFAM/1000,10).EQ.0) MSGN=
495 & MSGN*(-1)**MOD(KFAM/100,10)
496 ENDIF
497 ENDIF
498
499C...Kinematics of one-particle decays.
500 IF(ND.EQ.1) THEN
501 DO 400 J=1,4
502 P(N+1,J)=P(IP,J)
503 400 CONTINUE
504 GOTO 660
505 ENDIF
506
507C...Calculate maximum weight ND-particle decay.
508 PV(ND,5)=P(N+ND,5)
509 IF(ND.GE.3) THEN
510 WTMAX=1./WTCOR(ND-2)
511 PMAX=PV(1,5)-PS+P(N+ND,5)
512 PMIN=0.
513 DO 410 IL=ND-1,1,-1
514 PMAX=PMAX+P(N+IL,5)
515 PMIN=PMIN+P(N+IL+1,5)
516 WTMAX=WTMAX*PAWT(PMAX,PMIN,P(N+IL,5))
517 410 CONTINUE
518 ENDIF
519
520C...Find virtual gamma mass in Dalitz decay.
521 420 IF(ND.EQ.2) THEN
522 ELSEIF(MMAT.EQ.2) THEN
523 PMES=4.*PMAS(11,1)**2
524 PMRHO2=PMAS(131,1)**2
525 PGRHO2=PMAS(131,2)**2
526 430 PMST=PMES*(P(IP,5)**2/PMES)**RLU(0)
527 WT=(1+0.5*PMES/PMST)*SQRT(MAX(0.,1.-PMES/PMST))*
528 & (1.-PMST/P(IP,5)**2)**3*(1.+PGRHO2/PMRHO2)/
529 & ((1.-PMST/PMRHO2)**2+PGRHO2/PMRHO2)
530 IF(WT.LT.RLU(0)) GOTO 430
531 PV(2,5)=MAX(2.00001*PMAS(11,1),SQRT(PMST))
532
533C...M-generator gives weight. If rejected, try again.
534 ELSE
535 440 RORD(1)=1.
536 DO 470 IL1=2,ND-1
537 RSAV=RLU(0)
538 DO 450 IL2=IL1-1,1,-1
539 IF(RSAV.LE.RORD(IL2)) GOTO 460
540 RORD(IL2+1)=RORD(IL2)
541 450 CONTINUE
542 460 RORD(IL2+1)=RSAV
543 470 CONTINUE
544 RORD(ND)=0.
545 WT=1.
546 DO 480 IL=ND-1,1,-1
547 PV(IL,5)=PV(IL+1,5)+P(N+IL,5)+(RORD(IL)-RORD(IL+1))*(PV(1,5)-PS)
548 WT=WT*PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5))
549 480 CONTINUE
550 IF(WT.LT.RLU(0)*WTMAX) GOTO 440
551 ENDIF
552
553C...Perform two-particle decays in respective CM frame.
554 490 DO 510 IL=1,ND-1
555 PA=PAWT(PV(IL,5),PV(IL+1,5),P(N+IL,5))
556 UE(3)=2.*RLU(0)-1.
557 PHI=PARU(2)*RLU(0)
558 UE(1)=SQRT(1.-UE(3)**2)*COS(PHI)
559 UE(2)=SQRT(1.-UE(3)**2)*SIN(PHI)
560 DO 500 J=1,3
561 P(N+IL,J)=PA*UE(J)
562 PV(IL+1,J)=-PA*UE(J)
563 500 CONTINUE
564 P(N+IL,4)=SQRT(PA**2+P(N+IL,5)**2)
565 PV(IL+1,4)=SQRT(PA**2+PV(IL+1,5)**2)
566 510 CONTINUE
567
568C...Lorentz transform decay products to lab frame.
569 DO 520 J=1,4
570 P(N+ND,J)=PV(ND,J)
571 520 CONTINUE
572 DO 560 IL=ND-1,1,-1
573 DO 530 J=1,3
574 BE(J)=PV(IL,J)/PV(IL,4)
575 530 CONTINUE
576 GA=PV(IL,4)/PV(IL,5)
577 DO 550 I=N+IL,N+ND
578 BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
579 DO 540 J=1,3
580 P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
581 540 CONTINUE
582 P(I,4)=GA*(P(I,4)+BEP)
583 550 CONTINUE
584 560 CONTINUE
585
586C...Check that no infinite loop in matrix element weight.
587 NTRY=NTRY+1
588 IF(NTRY.GT.800) GOTO 590
589
590C...Matrix elements for omega and phi decays.
591 IF(MMAT.EQ.1) THEN
592 WT=(P(N+1,5)*P(N+2,5)*P(N+3,5))**2-(P(N+1,5)*FOUR(N+2,N+3))**2
593 & -(P(N+2,5)*FOUR(N+1,N+3))**2-(P(N+3,5)*FOUR(N+1,N+2))**2
594 & +2.*FOUR(N+1,N+2)*FOUR(N+1,N+3)*FOUR(N+2,N+3)
595 IF(MAX(WT*WTCOR(9)/P(IP,5)**6,0.001).LT.RLU(0)) GOTO 420
596
597C...Matrix elements for pi0 or eta Dalitz decay to gamma e+ e-.
598 ELSEIF(MMAT.EQ.2) THEN
599 FOUR12=FOUR(N+1,N+2)
600 FOUR13=FOUR(N+1,N+3)
601 WT=(PMST-0.5*PMES)*(FOUR12**2+FOUR13**2)+
602 & PMES*(FOUR12*FOUR13+FOUR12**2+FOUR13**2)
603 IF(WT.LT.RLU(0)*0.25*PMST*(P(IP,5)**2-PMST)**2) GOTO 490
604
605C...Matrix element for S0 -> S1 + V1 -> S1 + S2 + S3 (S scalar,
606C...V vector), of form cos**2(theta02) in V1 rest frame, and for
607C...S0 -> gamma + V1 -> gamma + S2 + S3, of form sin**2(theta02).
608 ELSEIF(MMAT.EQ.3.AND.NM.EQ.2) THEN
609 FOUR10=FOUR(IP,IM)
610 FOUR12=FOUR(IP,N+1)
611 FOUR02=FOUR(IM,N+1)
612 PMS1=P(IP,5)**2
613 PMS0=P(IM,5)**2
614 PMS2=P(N+1,5)**2
615 IF(KFAS.NE.22) HNUM=(FOUR10*FOUR12-PMS1*FOUR02)**2
616 IF(KFAS.EQ.22) HNUM=PMS1*(2.*FOUR10*FOUR12*FOUR02-
617 & PMS1*FOUR02**2-PMS0*FOUR12**2-PMS2*FOUR10**2+PMS1*PMS0*PMS2)
618 HNUM=MAX(1E-6*PMS1**2*PMS0*PMS2,HNUM)
619 HDEN=(FOUR10**2-PMS1*PMS0)*(FOUR12**2-PMS1*PMS2)
620 IF(HNUM.LT.RLU(0)*HDEN) GOTO 490
621
622C...Matrix element for "onium" -> g + g + g or gamma + g + g.
623 ELSEIF(MMAT.EQ.4) THEN
624 HX1=2.*FOUR(IP,N+1)/P(IP,5)**2
625 HX2=2.*FOUR(IP,N+2)/P(IP,5)**2
626 HX3=2.*FOUR(IP,N+3)/P(IP,5)**2
627 WT=((1.-HX1)/(HX2*HX3))**2+((1.-HX2)/(HX1*HX3))**2+
628 & ((1.-HX3)/(HX1*HX2))**2
629 IF(WT.LT.2.*RLU(0)) GOTO 420
630 IF(K(IP+1,2).EQ.22.AND.(1.-HX1)*P(IP,5)**2.LT.4.*PARJ(32)**2)
631 & GOTO 420
632
633C...Effective matrix element for nu spectrum in tau -> nu + hadrons.
634 ELSEIF(MMAT.EQ.41) THEN
635 HX1=2.*FOUR(IP,N+1)/P(IP,5)**2
636 HXM=MIN(0.75,2.*(1.-PS/P(IP,5)))
637 IF(HX1*(3.-2.*HX1).LT.RLU(0)*HXM*(3.-2.*HXM)) GOTO 420
638
639C...Matrix elements for weak decays (only semileptonic for c and b)
640 ELSEIF((MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.44.OR.MMAT.EQ.48)
641 &.AND.ND.EQ.3) THEN
642 IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+3)
643 IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+3)
644 IF(WT.LT.RLU(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 420
645 ELSEIF(MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.44.OR.MMAT.EQ.48) THEN
646 DO 580 J=1,4
647 P(N+NP+1,J)=0.
648 DO 570 IS=N+3,N+NP
649 P(N+NP+1,J)=P(N+NP+1,J)+P(IS,J)
650 570 CONTINUE
651 580 CONTINUE
652 IF(MBST.EQ.0) WT=FOUR(IP,N+1)*FOUR(N+2,N+NP+1)
653 IF(MBST.EQ.1) WT=P(IP,5)*P(N+1,4)*FOUR(N+2,N+NP+1)
654 IF(WT.LT.RLU(0)*P(IP,5)*PV(1,5)**3/WTCOR(10)) GOTO 420
655
656C...Angular distribution in W decay.
657 ELSEIF(MMAT.EQ.46.AND.MSGN.NE.0) THEN
658 IF(MSGN.GT.0) WT=FOUR(IM,N+1)*FOUR(N+2,IP+1)
659 IF(MSGN.LT.0) WT=FOUR(IM,N+2)*FOUR(N+1,IP+1)
660 IF(WT.LT.RLU(0)*P(IM,5)**4/WTCOR(10)) GOTO 490
661 ENDIF
662
663C...Scale back energy and reattach spectator.
664 590 IF(MREM.EQ.1) THEN
665 DO 600 J=1,5
666 PV(1,J)=PV(1,J)/(1.-PQT)
667 600 CONTINUE
668 ND=ND+1
669 MREM=0
670 ENDIF
671
672C...Low invariant mass for system with spectator quark gives particle,
673C...not two jets. Readjust momenta accordingly.
674 IF((MMAT.EQ.31.OR.MMAT.EQ.45).AND.ND.EQ.3) THEN
675 MSTJ(93)=1
676 PM2=ULMASS(K(N+2,2))
677 MSTJ(93)=1
678 PM3=ULMASS(K(N+3,2))
679 IF(P(N+2,5)**2+P(N+3,5)**2+2.*FOUR(N+2,N+3).GE.
680 & (PARJ(32)+PM2+PM3)**2) GOTO 660
681 K(N+2,1)=1
682 KFTEMP=K(N+2,2)
683 CALL LUKFDI(KFTEMP,K(N+3,2),KFLDMP,K(N+2,2))
684 IF(K(N+2,2).EQ.0) GOTO 260
685 P(N+2,5)=ULMASS(K(N+2,2))
686 PS=P(N+1,5)+P(N+2,5)
687 PV(2,5)=P(N+2,5)
688 MMAT=0
689 ND=2
690 GOTO 490
691 ELSEIF(MMAT.EQ.44) THEN
692 MSTJ(93)=1
693 PM3=ULMASS(K(N+3,2))
694 MSTJ(93)=1
695 PM4=ULMASS(K(N+4,2))
696 IF(P(N+3,5)**2+P(N+4,5)**2+2.*FOUR(N+3,N+4).GE.
697 & (PARJ(32)+PM3+PM4)**2) GOTO 630
698 K(N+3,1)=1
699 KFTEMP=K(N+3,2)
700 CALL LUKFDI(KFTEMP,K(N+4,2),KFLDMP,K(N+3,2))
701 IF(K(N+3,2).EQ.0) GOTO 260
702 P(N+3,5)=ULMASS(K(N+3,2))
703 DO 610 J=1,3
704 P(N+3,J)=P(N+3,J)+P(N+4,J)
705 610 CONTINUE
706 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)
707 HA=P(N+1,4)**2-P(N+2,4)**2
708 HB=HA-(P(N+1,5)**2-P(N+2,5)**2)
709 HC=(P(N+1,1)-P(N+2,1))**2+(P(N+1,2)-P(N+2,2))**2+
710 & (P(N+1,3)-P(N+2,3))**2
711 HD=(PV(1,4)-P(N+3,4))**2
712 HE=HA**2-2.*HD*(P(N+1,4)**2+P(N+2,4)**2)+HD**2
713 HF=HD*HC-HB**2
714 HG=HD*HC-HA*HB
715 HH=(SQRT(HG**2+HE*HF)-HG)/(2.*HF)
716 DO 620 J=1,3
717 PCOR=HH*(P(N+1,J)-P(N+2,J))
718 P(N+1,J)=P(N+1,J)+PCOR
719 P(N+2,J)=P(N+2,J)-PCOR
720 620 CONTINUE
721 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)
722 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)
723 ND=ND-1
724 ENDIF
725
726C...Check invariant mass of W jets. May give one particle or start over.
727 630 IF((MMAT.EQ.42.OR.MMAT.EQ.43.OR.MMAT.EQ.44.OR.MMAT.EQ.48)
728 &.AND.IABS(K(N+1,2)).LT.10) THEN
729 PMR=SQRT(MAX(0.,P(N+1,5)**2+P(N+2,5)**2+2.*FOUR(N+1,N+2)))
730 MSTJ(93)=1
731 PM1=ULMASS(K(N+1,2))
732 MSTJ(93)=1
733 PM2=ULMASS(K(N+2,2))
734 IF(PMR.GT.PARJ(32)+PM1+PM2) GOTO 640
735 KFLDUM=INT(1.5+RLU(0))
736 CALL LUKFDI(K(N+1,2),-ISIGN(KFLDUM,K(N+1,2)),KFLDMP,KF1)
737 CALL LUKFDI(K(N+2,2),-ISIGN(KFLDUM,K(N+2,2)),KFLDMP,KF2)
738 IF(KF1.EQ.0.OR.KF2.EQ.0) GOTO 260
739 PSM=ULMASS(KF1)+ULMASS(KF2)
740 IF((MMAT.EQ.42.OR.MMAT.EQ.48).AND.PMR.GT.PARJ(64)+PSM) GOTO 640
741 IF(MMAT.GE.43.AND.PMR.GT.0.2*PARJ(32)+PSM) GOTO 640
742 IF(MMAT.EQ.48) GOTO 420
743 IF(ND.EQ.4.OR.KFA.EQ.15) GOTO 260
744 K(N+1,1)=1
745 KFTEMP=K(N+1,2)
746 CALL LUKFDI(KFTEMP,K(N+2,2),KFLDMP,K(N+1,2))
747 IF(K(N+1,2).EQ.0) GOTO 260
748 P(N+1,5)=ULMASS(K(N+1,2))
749 K(N+2,2)=K(N+3,2)
750 P(N+2,5)=P(N+3,5)
751 PS=P(N+1,5)+P(N+2,5)
752 IF(PS+PARJ(64).GT.PV(1,5)) GOTO 260
753 PV(2,5)=P(N+3,5)
754 MMAT=0
755 ND=2
756 GOTO 490
757 ENDIF
758
759C...Phase space decay of partons from W decay.
760 640 IF((MMAT.EQ.42.OR.MMAT.EQ.48).AND.IABS(K(N+1,2)).LT.10) THEN
761 KFLO(1)=K(N+1,2)
762 KFLO(2)=K(N+2,2)
763 K(N+1,1)=K(N+3,1)
764 K(N+1,2)=K(N+3,2)
765 DO 650 J=1,5
766 PV(1,J)=P(N+1,J)+P(N+2,J)
767 P(N+1,J)=P(N+3,J)
768 650 CONTINUE
769 PV(1,5)=PMR
770 N=N+1
771 NP=0
772 NQ=2
773 PS=0.
774 MSTJ(93)=2
775 PSQ=ULMASS(KFLO(1))
776 MSTJ(93)=2
777 PSQ=PSQ+ULMASS(KFLO(2))
778 MMAT=11
779 GOTO 290
780 ENDIF
781
782C...Boost back for rapidly moving particle.
783 660 N=N+ND
784 IF(MBST.EQ.1) THEN
785 DO 670 J=1,3
786 BE(J)=P(IP,J)/P(IP,4)
787 670 CONTINUE
788 GA=P(IP,4)/P(IP,5)
789 DO 690 I=NSAV+1,N
790 BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
791 DO 680 J=1,3
792 P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
793 680 CONTINUE
794 P(I,4)=GA*(P(I,4)+BEP)
795 690 CONTINUE
796 ENDIF
797
798C...Fill in position of decay vertex.
799 DO 710 I=NSAV+1,N
800 DO 700 J=1,4
801 V(I,J)=VDCY(J)
802 700 CONTINUE
803 V(I,5)=0.
804 710 CONTINUE
805
806C...Set up for parton shower evolution from jets.
807 IF(MSTJ(23).GE.1.AND.MMAT.EQ.4.AND.K(NSAV+1,2).EQ.21) THEN
808 K(NSAV+1,1)=3
809 K(NSAV+2,1)=3
810 K(NSAV+3,1)=3
811 K(NSAV+1,4)=MSTU(5)*(NSAV+2)
812 K(NSAV+1,5)=MSTU(5)*(NSAV+3)
813 K(NSAV+2,4)=MSTU(5)*(NSAV+3)
814 K(NSAV+2,5)=MSTU(5)*(NSAV+1)
815 K(NSAV+3,4)=MSTU(5)*(NSAV+1)
816 K(NSAV+3,5)=MSTU(5)*(NSAV+2)
817 MSTJ(92)=-(NSAV+1)
818 ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.4) THEN
819 K(NSAV+2,1)=3
820 K(NSAV+3,1)=3
821 K(NSAV+2,4)=MSTU(5)*(NSAV+3)
822 K(NSAV+2,5)=MSTU(5)*(NSAV+3)
823 K(NSAV+3,4)=MSTU(5)*(NSAV+2)
824 K(NSAV+3,5)=MSTU(5)*(NSAV+2)
825 MSTJ(92)=NSAV+2
826 ELSEIF(MSTJ(23).GE.1.AND.(MMAT.EQ.32.OR.MMAT.EQ.44.OR.MMAT.EQ.46)
827 &.AND.IABS(K(NSAV+1,2)).LE.10.AND.IABS(K(NSAV+2,2)).LE.10) THEN
828 K(NSAV+1,1)=3
829 K(NSAV+2,1)=3
830 K(NSAV+1,4)=MSTU(5)*(NSAV+2)
831 K(NSAV+1,5)=MSTU(5)*(NSAV+2)
832 K(NSAV+2,4)=MSTU(5)*(NSAV+1)
833 K(NSAV+2,5)=MSTU(5)*(NSAV+1)
834 MSTJ(92)=NSAV+1
835 ELSEIF(MSTJ(23).GE.1.AND.(MMAT.EQ.32.OR.MMAT.EQ.44.OR.MMAT.EQ.46)
836 &.AND.IABS(K(NSAV+1,2)).LE.20.AND.IABS(K(NSAV+2,2)).LE.20) THEN
837 MSTJ(92)=NSAV+1
838 ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33.AND.IABS(K(NSAV+2,2)).EQ.21)
839 &THEN
840 K(NSAV+1,1)=3
841 K(NSAV+2,1)=3
842 K(NSAV+3,1)=3
843 KCP=LUCOMP(K(NSAV+1,2))
844 KQP=KCHG(KCP,2)*ISIGN(1,K(NSAV+1,2))
845 JCON=4
846 IF(KQP.LT.0) JCON=5
847 K(NSAV+1,JCON)=MSTU(5)*(NSAV+2)
848 K(NSAV+2,9-JCON)=MSTU(5)*(NSAV+1)
849 K(NSAV+2,JCON)=MSTU(5)*(NSAV+3)
850 K(NSAV+3,9-JCON)=MSTU(5)*(NSAV+2)
851 MSTJ(92)=NSAV+1
852 ELSEIF(MSTJ(23).GE.1.AND.MMAT.EQ.33) THEN
853 K(NSAV+1,1)=3
854 K(NSAV+3,1)=3
855 K(NSAV+1,4)=MSTU(5)*(NSAV+3)
856 K(NSAV+1,5)=MSTU(5)*(NSAV+3)
857 K(NSAV+3,4)=MSTU(5)*(NSAV+1)
858 K(NSAV+3,5)=MSTU(5)*(NSAV+1)
859 MSTJ(92)=NSAV+1
860
861C...Set up for parton shower evolution in t -> W + b.
862 ELSEIF(MSTJ(27).GE.1.AND.MMAT.EQ.45.AND.ND.EQ.3) THEN
863 K(NSAV+2,1)=3
864 K(NSAV+3,1)=3
865 K(NSAV+2,4)=MSTU(5)*(NSAV+3)
866 K(NSAV+2,5)=MSTU(5)*(NSAV+3)
867 K(NSAV+3,4)=MSTU(5)*(NSAV+2)
868 K(NSAV+3,5)=MSTU(5)*(NSAV+2)
869 MSTJ(92)=NSAV+1
870 ENDIF
871
872C...Mark decayed particle; special option for B-B~ mixing.
873 IF(K(IP,1).EQ.5) K(IP,1)=15
874 IF(K(IP,1).LE.10) K(IP,1)=11
875 IF(MMIX.EQ.1.AND.MSTJ(26).EQ.2.AND.K(IP,1).EQ.11) K(IP,1)=12
876 K(IP,4)=NSAV+1
877 K(IP,5)=N
878
879 RETURN
880 END