Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pydecy.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pydecy.f
1 
2 C*********************************************************************
3 
4 C...PYDECY
5 C...Handles the decay of unstable particles.
6 
7  SUBROUTINE pydecy(IP)
8 
9 C...Double precision and integer declarations.
10  IMPLICIT DOUBLE PRECISION(a-h, o-z)
11  IMPLICIT INTEGER(i-n)
12  INTEGER pyk,pychge,pycomp
13 C...Commonblocks.
14  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
15  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
16  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
17  common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
18  SAVE /pyjets/,/pydat1/,/pydat2/,/pydat3/
19 C...Local arrays.
20  dimension vdcy(4),kflo(4),kfl1(4),pv(10,5),rord(10),ue(3),be(3),
21  &wtcor(10),ptau(4),pcmtau(4),dbetau(3)
22  CHARACTER cidc*4
23  DATA wtcor/2d0,5d0,15d0,60d0,250d0,1500d0,1.2d4,1.2d5,150d0,16d0/
24 
25 C...Functions: momentum in two-particle decays and four-product.
26  pawt(a,b,c)=sqrt((a**2-(b+c)**2)*(a**2-(b-c)**2))/(2d0*a)
27  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)
28 
29 C...Initial values.
30  ntry=0
31  nsav=n
32  kfa=iabs(k(ip,2))
33  kfs=isign(1,k(ip,2))
34  kc=pycomp(kfa)
35  mstj(92)=0
36 
37 C...Choose lifetime and determine decay vertex.
38  IF(k(ip,1).EQ.5) THEN
39  v(ip,5)=0d0
40  ELSEIF(k(ip,1).NE.4) THEN
41  v(ip,5)=-pmas(kc,4)*log(pyr(0))
42  ENDIF
43  DO 100 j=1,4
44  vdcy(j)=v(ip,j)+v(ip,5)*p(ip,j)/p(ip,5)
45  100 CONTINUE
46 
47 C...Determine whether decay allowed or not.
48  mout=0
49  IF(mstj(22).EQ.2) THEN
50  IF(pmas(kc,4).GT.parj(71)) mout=1
51  ELSEIF(mstj(22).EQ.3) THEN
52  IF(vdcy(1)**2+vdcy(2)**2+vdcy(3)**2.GT.parj(72)**2) mout=1
53  ELSEIF(mstj(22).EQ.4) THEN
54  IF(vdcy(1)**2+vdcy(2)**2.GT.parj(73)**2) mout=1
55  IF(abs(vdcy(3)).GT.parj(74)) mout=1
56  ENDIF
57  IF(mout.EQ.1.AND.k(ip,1).NE.5) THEN
58  k(ip,1)=4
59  RETURN
60  ENDIF
61 
62 C...Interface to external tau decay library (for tau polarization).
63  IF(kfa.EQ.15.AND.mstj(28).GE.1) THEN
64 
65 C...Starting values for pointers and momenta.
66  itau=ip
67  DO 110 j=1,4
68  ptau(j)=p(itau,j)
69  pcmtau(j)=p(itau,j)
70  110 CONTINUE
71 
72 C...Iterate to find position and code of mother of tau.
73  imtau=itau
74  120 imtau=k(imtau,3)
75 
76  IF(imtau.EQ.0) THEN
77 C...If no known origin then impossible to do anything further.
78  kforig=0
79  iorig=0
80 
81  ELSEIF(k(imtau,2).EQ.k(itau,2)) THEN
82 C...If tau -> tau + gamma then add gamma energy and loop.
83  IF(k(k(imtau,4),2).EQ.22) THEN
84  DO 130 j=1,4
85  pcmtau(j)=pcmtau(j)+p(k(imtau,4),j)
86  130 CONTINUE
87  ELSEIF(k(k(imtau,5),2).EQ.22) THEN
88  DO 140 j=1,4
89  pcmtau(j)=pcmtau(j)+p(k(imtau,5),j)
90  140 CONTINUE
91  ENDIF
92  goto 120
93 
94  ELSEIF(iabs(k(imtau,2)).GT.100) THEN
95 C...If coming from weak decay of hadron then W is not stored in record,
96 C...but can be reconstructed by adding neutrino momentum.
97  kforig=-isign(24,k(itau,2))
98  iorig=0
99  DO 160 ii=k(imtau,4),k(imtau,5)
100  IF(k(ii,2)*isign(1,k(itau,2)).EQ.-16) THEN
101  DO 150 j=1,4
102  pcmtau(j)=pcmtau(j)+p(ii,j)
103  150 CONTINUE
104  ENDIF
105  160 CONTINUE
106 
107  ELSE
108 C...If coming from resonance decay then find latest copy of this
109 C...resonance (may not completely agree).
110  kforig=k(imtau,2)
111  iorig=imtau
112  DO 170 ii=imtau+1,ip-1
113  IF(k(ii,2).EQ.kforig.AND.k(ii,3).EQ.iorig.AND.
114  & abs(p(ii,5)-p(iorig,5)).LT.1d-5*p(iorig,5)) iorig=ii
115  170 CONTINUE
116  DO 180 j=1,4
117  pcmtau(j)=p(iorig,j)
118  180 CONTINUE
119  ENDIF
120 
121 C...Boost tau to rest frame of production process (where known)
122 C...and rotate it to sit along +z axis.
123  DO 190 j=1,3
124  dbetau(j)=pcmtau(j)/pcmtau(4)
125  190 CONTINUE
126  IF(kforig.NE.0) CALL pyrobo(itau,itau,0d0,0d0,-dbetau(1),
127  & -dbetau(2),-dbetau(3))
128  phitau=pyangl(p(itau,1),p(itau,2))
129  CALL pyrobo(itau,itau,0d0,-phitau,0d0,0d0,0d0)
130  thetau=pyangl(p(itau,3),p(itau,1))
131  CALL pyrobo(itau,itau,-thetau,0d0,0d0,0d0,0d0)
132 
133 C...Call tau decay routine (if meaningful) and fill extra info.
134  IF(kforig.NE.0.OR.mstj(28).EQ.2) THEN
135  CALL pytaud(itau,iorig,kforig,ndecay)
136  DO 200 ii=nsav+1,nsav+ndecay
137  k(ii,1)=1
138  k(ii,3)=ip
139  k(ii,4)=0
140  k(ii,5)=0
141  200 CONTINUE
142  n=nsav+ndecay
143  ENDIF
144 
145 C...Boost back decay tau and decay products.
146  DO 210 j=1,4
147  p(itau,j)=ptau(j)
148  210 CONTINUE
149  IF(kforig.NE.0.OR.mstj(28).EQ.2) THEN
150  CALL pyrobo(nsav+1,n,thetau,phitau,0d0,0d0,0d0)
151  IF(kforig.NE.0) CALL pyrobo(nsav+1,n,0d0,0d0,dbetau(1),
152  & dbetau(2),dbetau(3))
153 
154 C...Skip past ordinary tau decay treatment.
155  mmat=0
156  mbst=0
157  nd=0
158  goto 630
159  ENDIF
160  ENDIF
161 
162 C...B-Bbar mixing: flip sign of meson appropriately.
163  mmix=0
164  IF((kfa.EQ.511.OR.kfa.EQ.531).AND.mstj(26).GE.1) THEN
165  xbbmix=parj(76)
166  IF(kfa.EQ.531) xbbmix=parj(77)
167  IF(sin(0.5d0*xbbmix*v(ip,5)/pmas(kc,4))**2.GT.pyr(0)) mmix=1
168  IF(mmix.EQ.1) kfs=-kfs
169  ENDIF
170 
171 C...Check existence of decay channels. Particle/antiparticle rules.
172  kca=kc
173  IF(mdcy(kc,2).GT.0) THEN
174  mdmdcy=mdme(mdcy(kc,2),2)
175  IF(mdmdcy.GT.80.AND.mdmdcy.LE.90) kca=mdmdcy
176  ENDIF
177  IF(mdcy(kca,2).LE.0.OR.mdcy(kca,3).LE.0) THEN
178  CALL pyerrm(9,'(PYDECY:) no decay channel defined')
179  RETURN
180  ENDIF
181  IF(mod(kfa/1000,10).EQ.0.AND.kca.EQ.85) kfs=-kfs
182  IF(kchg(kc,3).EQ.0) THEN
183  kfsp=1
184  kfsn=0
185  IF(pyr(0).GT.0.5d0) kfs=-kfs
186  ELSEIF(kfs.GT.0) THEN
187  kfsp=1
188  kfsn=0
189  ELSE
190  kfsp=0
191  kfsn=1
192  ENDIF
193 
194 C...Sum branching ratios of allowed decay channels.
195  220 nope=0
196  brsu=0d0
197  DO 230 idl=mdcy(kca,2),mdcy(kca,2)+mdcy(kca,3)-1
198  IF(mdme(idl,1).NE.1.AND.kfsp*mdme(idl,1).NE.2.AND.
199  & kfsn*mdme(idl,1).NE.3) goto 230
200  IF(mdme(idl,2).GT.100) goto 230
201  nope=nope+1
202  brsu=brsu+brat(idl)
203  230 CONTINUE
204  IF(nope.EQ.0) THEN
205  CALL pyerrm(2,'(PYDECY:) all decay channels closed by user')
206  RETURN
207  ENDIF
208 
209 C...Select decay channel among allowed ones.
210  240 rbr=brsu*pyr(0)
211  idl=mdcy(kca,2)-1
212  250 idl=idl+1
213  IF(mdme(idl,1).NE.1.AND.kfsp*mdme(idl,1).NE.2.AND.
214  &kfsn*mdme(idl,1).NE.3) THEN
215  IF(idl.LT.mdcy(kca,2)+mdcy(kca,3)-1) goto 250
216  ELSEIF(mdme(idl,2).GT.100) THEN
217  IF(idl.LT.mdcy(kca,2)+mdcy(kca,3)-1) goto 250
218  ELSE
219  idc=idl
220  rbr=rbr-brat(idl)
221  IF(idl.LT.mdcy(kca,2)+mdcy(kca,3)-1.AND.rbr.GT.0d0) goto 250
222  ENDIF
223 
224 C...Start readout of decay channel: matrix element, reset counters.
225  mmat=mdme(idc,2)
226  260 ntry=ntry+1
227  IF(mod(ntry,200).EQ.0) THEN
228  WRITE(cidc,'(I4)') idc
229 C...Do not print warning for some well-known special cases.
230  IF(kfa.NE.113.AND.kfa.NE.115.AND.kfa.NE.215)
231  & CALL pyerrm(4,'(PYDECY:) caught in loop for decay channel'//
232  & cidc)
233  goto 240
234  ENDIF
235  IF(ntry.GT.1000) THEN
236  CALL pyerrm(14,'(PYDECY:) caught in infinite loop')
237  IF(mstu(21).GE.1) RETURN
238  ENDIF
239  i=n
240  np=0
241  nq=0
242  mbst=0
243  IF(mmat.GE.11.AND.p(ip,4).GT.20d0*p(ip,5)) mbst=1
244  DO 270 j=1,4
245  pv(1,j)=0d0
246  IF(mbst.EQ.0) pv(1,j)=p(ip,j)
247  270 CONTINUE
248  IF(mbst.EQ.1) pv(1,4)=p(ip,5)
249  pv(1,5)=p(ip,5)
250  ps=0d0
251  psq=0d0
252  mrem=0
253  mhaddy=0
254  IF(kfa.GT.80) mhaddy=1
255 C.. Random flavour and popcorn system memory.
256  irndmo=0
257  jtmo=0
258  mstu(121)=0
259  mstu(125)=10
260 
261 C...Read out decay products. Convert to standard flavour code.
262  jtmax=5
263  IF(mdme(idc+1,2).EQ.101) jtmax=10
264  DO 280 jt=1,jtmax
265  IF(jt.LE.5) kp=kfdp(idc,jt)
266  IF(jt.GE.6) kp=kfdp(idc+1,jt-5)
267  IF(kp.EQ.0) goto 280
268  kpa=iabs(kp)
269  kcp=pycomp(kpa)
270  IF(kpa.GT.80) mhaddy=1
271  IF(kchg(kcp,3).EQ.0.AND.kpa.NE.81.AND.kpa.NE.82) THEN
272  kfp=kp
273  ELSEIF(kpa.NE.81.AND.kpa.NE.82) THEN
274  kfp=kfs*kp
275  ELSEIF(kpa.EQ.81.AND.mod(kfa/1000,10).EQ.0) THEN
276  kfp=-kfs*mod(kfa/10,10)
277  ELSEIF(kpa.EQ.81.AND.mod(kfa/100,10).GE.mod(kfa/10,10)) THEN
278  kfp=kfs*(100*mod(kfa/10,100)+3)
279  ELSEIF(kpa.EQ.81) THEN
280  kfp=kfs*(1000*mod(kfa/10,10)+100*mod(kfa/100,10)+1)
281  ELSEIF(kp.EQ.82) THEN
282  CALL pydcyk(-kfs*int(1d0+(2d0+parj(2))*pyr(0)),0,kfp,kdump)
283  IF(kfp.EQ.0) goto 260
284  kfp=-kfp
285  irndmo=1
286  mstj(93)=1
287  IF(pv(1,5).LT.parj(32)+2d0*pymass(kfp)) goto 260
288  ELSEIF(kp.EQ.-82) THEN
289  kfp=mstu(124)
290  ENDIF
291  IF(kpa.EQ.81.OR.kpa.EQ.82) kcp=pycomp(kfp)
292 
293 C...Add decay product to event record or to quark flavour list.
294  kfpa=iabs(kfp)
295  kqp=kchg(kcp,2)
296  IF(mmat.GE.11.AND.mmat.LE.30.AND.kqp.NE.0) THEN
297  nq=nq+1
298  kflo(nq)=kfp
299 C...set rndmflav popcorn system pointer
300  IF(kp.EQ.82.AND.mstu(121).GT.0) jtmo=nq
301  mstj(93)=2
302  psq=psq+pymass(kflo(nq))
303  ELSEIF((mmat.EQ.42.OR.mmat.EQ.43.OR.mmat.EQ.48).AND.np.EQ.3.AND.
304  & mod(nq,2).EQ.1) THEN
305  nq=nq-1
306  ps=ps-p(i,5)
307  k(i,1)=1
308  kfi=k(i,2)
309  CALL pykfdi(kfp,kfi,kfldmp,k(i,2))
310  IF(k(i,2).EQ.0) goto 260
311  mstj(93)=1
312  p(i,5)=pymass(k(i,2))
313  ps=ps+p(i,5)
314  ELSE
315  i=i+1
316  np=np+1
317  IF(mmat.NE.33.AND.kqp.NE.0) nq=nq+1
318  IF(mmat.EQ.33.AND.kqp.NE.0.AND.kqp.NE.2) nq=nq+1
319  k(i,1)=1+mod(nq,2)
320  IF(mmat.EQ.4.AND.jt.LE.2.AND.kfp.EQ.21) k(i,1)=2
321  IF(mmat.EQ.4.AND.jt.EQ.3) k(i,1)=1
322  k(i,2)=kfp
323  k(i,3)=ip
324  k(i,4)=0
325  k(i,5)=0
326  p(i,5)=pymass(kfp)
327  ps=ps+p(i,5)
328  ENDIF
329  280 CONTINUE
330 
331 C...Check masses for resonance decays.
332  IF(mhaddy.EQ.0) THEN
333  IF(ps+parj(64).GT.pv(1,5)) goto 240
334  ENDIF
335 
336 C...Choose decay multiplicity in phase space model.
337  290 IF(mmat.GE.11.AND.mmat.LE.30) THEN
338  psp=ps
339  cnde=parj(61)*log(max((pv(1,5)-ps-psq)/parj(62),1.1d0))
340  IF(mmat.EQ.12) cnde=cnde+parj(63)
341  300 ntry=ntry+1
342 C...Reset popcorn flags if new attempt. Re-select rndmflav if failed.
343  IF(irndmo.EQ.0) THEN
344  mstu(121)=0
345  jtmo=0
346  ELSEIF(irndmo.EQ.1) THEN
347  irndmo=2
348  ELSE
349  goto 260
350  ENDIF
351  IF(ntry.GT.1000) THEN
352  CALL pyerrm(14,'(PYDECY:) caught in infinite loop')
353  IF(mstu(21).GE.1) RETURN
354  ENDIF
355  IF(mmat.LE.20) THEN
356  gauss=sqrt(-2d0*cnde*log(max(1d-10,pyr(0))))*
357  & sin(paru(2)*pyr(0))
358  nd=0.5d0+0.5d0*np+0.25d0*nq+cnde+gauss
359  IF(nd.LT.np+nq/2.OR.nd.LT.2.OR.nd.GT.10) goto 300
360  IF(mmat.EQ.13.AND.nd.EQ.2) goto 300
361  IF(mmat.EQ.14.AND.nd.LE.3) goto 300
362  IF(mmat.EQ.15.AND.nd.LE.4) goto 300
363  ELSE
364  nd=mmat-20
365  ENDIF
366 C.. Set maximum popcorn meson number. Test rndmflav popcorn size.
367  mstu(125)=nd-nq/2
368  IF(mstu(121).GT.mstu(125)) goto 300
369 
370 C...Form hadrons from flavour content.
371  DO 310 jt=1,nq
372  kfl1(jt)=kflo(jt)
373  310 CONTINUE
374  IF(nd.EQ.np+nq/2) goto 330
375  DO 320 i=n+np+1,n+nd-nq/2
376 C.. Stick to started popcorn system, else pick side at random
377  jt=jtmo
378  IF(jt.EQ.0) jt=1+int((nq-1)*pyr(0))
379  CALL pydcyk(kfl1(jt),0,kfl2,k(i,2))
380  IF(k(i,2).EQ.0) goto 300
381  mstu(125)=mstu(125)-1
382  jtmo=0
383  IF(mstu(121).GT.0) jtmo=jt
384  kfl1(jt)=-kfl2
385  320 CONTINUE
386  330 jt=2
387  jt2=3
388  jt3=4
389  IF(nq.EQ.4.AND.pyr(0).LT.parj(66)) jt=4
390  IF(jt.EQ.4.AND.isign(1,kfl1(1)*(10-iabs(kfl1(1))))*
391  & isign(1,kfl1(jt)*(10-iabs(kfl1(jt)))).GT.0) jt=3
392  IF(jt.EQ.3) jt2=2
393  IF(jt.EQ.4) jt3=2
394  CALL pydcyk(kfl1(1),kfl1(jt),kfldmp,k(n+nd-nq/2+1,2))
395  IF(k(n+nd-nq/2+1,2).EQ.0) goto 300
396  IF(nq.EQ.4) CALL pydcyk(kfl1(jt2),kfl1(jt3),kfldmp,k(n+nd,2))
397  IF(nq.EQ.4.AND.k(n+nd,2).EQ.0) goto 300
398 
399 C...Check that sum of decay product masses not too large.
400  ps=psp
401  DO 340 i=n+np+1,n+nd
402  k(i,1)=1
403  k(i,3)=ip
404  k(i,4)=0
405  k(i,5)=0
406  p(i,5)=pymass(k(i,2))
407  ps=ps+p(i,5)
408  340 CONTINUE
409  IF(ps+parj(64).GT.pv(1,5)) goto 300
410 
411 C...Rescale energy to subtract off spectator quark mass.
412  ELSEIF((mmat.EQ.31.OR.mmat.EQ.33.OR.mmat.EQ.44)
413  & .AND.np.GE.3) THEN
414  ps=ps-p(n+np,5)
415  pqt=(p(n+np,5)+parj(65))/pv(1,5)
416  DO 350 j=1,5
417  p(n+np,j)=pqt*pv(1,j)
418  pv(1,j)=(1d0-pqt)*pv(1,j)
419  350 CONTINUE
420  IF(ps+parj(64).GT.pv(1,5)) goto 260
421  nd=np-1
422  mrem=1
423 
424 C...Fully specified final state: check mass broadening effects.
425  ELSE
426  IF(np.GE.2.AND.ps+parj(64).GT.pv(1,5)) goto 260
427  nd=np
428  ENDIF
429 
430 C...Determine position of grandmother, number of sisters.
431  nm=0
432  kfas=0
433  msgn=0
434  IF(mmat.EQ.3) THEN
435  im=k(ip,3)
436  IF(im.LT.0.OR.im.GE.ip) im=0
437  IF(im.NE.0) kfam=iabs(k(im,2))
438  IF(im.NE.0) THEN
439  DO 360 il=max(ip-2,im+1),min(ip+2,n)
440  IF(k(il,3).EQ.im) nm=nm+1
441  IF(k(il,3).EQ.im.AND.il.NE.ip) isis=il
442  360 CONTINUE
443  IF(nm.NE.2.OR.kfam.LE.100.OR.mod(kfam,10).NE.1.OR.
444  & mod(kfam/1000,10).NE.0) nm=0
445  IF(nm.EQ.2) THEN
446  kfas=iabs(k(isis,2))
447  IF((kfas.LE.100.OR.mod(kfas,10).NE.1.OR.
448  & mod(kfas/1000,10).NE.0).AND.kfas.NE.22) nm=0
449  ENDIF
450  ENDIF
451  ENDIF
452 
453 C...Kinematics of one-particle decays.
454  IF(nd.EQ.1) THEN
455  DO 370 j=1,4
456  p(n+1,j)=p(ip,j)
457  370 CONTINUE
458  goto 630
459  ENDIF
460 
461 C...Calculate maximum weight ND-particle decay.
462  pv(nd,5)=p(n+nd,5)
463  IF(nd.GE.3) THEN
464  wtmax=1d0/wtcor(nd-2)
465  pmax=pv(1,5)-ps+p(n+nd,5)
466  pmin=0d0
467  DO 380 il=nd-1,1,-1
468  pmax=pmax+p(n+il,5)
469  pmin=pmin+p(n+il+1,5)
470  wtmax=wtmax*pawt(pmax,pmin,p(n+il,5))
471  380 CONTINUE
472  ENDIF
473 
474 C...Find virtual gamma mass in Dalitz decay.
475  390 IF(nd.EQ.2) THEN
476  ELSEIF(mmat.EQ.2) THEN
477  pmes=4d0*pmas(11,1)**2
478  pmrho2=pmas(131,1)**2
479  pgrho2=pmas(131,2)**2
480  400 pmst=pmes*(p(ip,5)**2/pmes)**pyr(0)
481  wt=(1+0.5d0*pmes/pmst)*sqrt(max(0d0,1d0-pmes/pmst))*
482  & (1d0-pmst/p(ip,5)**2)**3*(1d0+pgrho2/pmrho2)/
483  & ((1d0-pmst/pmrho2)**2+pgrho2/pmrho2)
484  IF(wt.LT.pyr(0)) goto 400
485  pv(2,5)=max(2.00001d0*pmas(11,1),sqrt(pmst))
486 
487 C...M-generator gives weight. If rejected, try again.
488  ELSE
489  410 rord(1)=1d0
490  DO 440 il1=2,nd-1
491  rsav=pyr(0)
492  DO 420 il2=il1-1,1,-1
493  IF(rsav.LE.rord(il2)) goto 430
494  rord(il2+1)=rord(il2)
495  420 CONTINUE
496  430 rord(il2+1)=rsav
497  440 CONTINUE
498  rord(nd)=0d0
499  wt=1d0
500  DO 450 il=nd-1,1,-1
501  pv(il,5)=pv(il+1,5)+p(n+il,5)+(rord(il)-rord(il+1))*
502  & (pv(1,5)-ps)
503  wt=wt*pawt(pv(il,5),pv(il+1,5),p(n+il,5))
504  450 CONTINUE
505  IF(wt.LT.pyr(0)*wtmax) goto 410
506  ENDIF
507 
508 C...Perform two-particle decays in respective CM frame.
509  460 DO 480 il=1,nd-1
510  pa=pawt(pv(il,5),pv(il+1,5),p(n+il,5))
511  ue(3)=2d0*pyr(0)-1d0
512  phi=paru(2)*pyr(0)
513  ue(1)=sqrt(1d0-ue(3)**2)*cos(phi)
514  ue(2)=sqrt(1d0-ue(3)**2)*sin(phi)
515  DO 470 j=1,3
516  p(n+il,j)=pa*ue(j)
517  pv(il+1,j)=-pa*ue(j)
518  470 CONTINUE
519  p(n+il,4)=sqrt(pa**2+p(n+il,5)**2)
520  pv(il+1,4)=sqrt(pa**2+pv(il+1,5)**2)
521  480 CONTINUE
522 
523 C...Lorentz transform decay products to lab frame.
524  DO 490 j=1,4
525  p(n+nd,j)=pv(nd,j)
526  490 CONTINUE
527  DO 530 il=nd-1,1,-1
528  DO 500 j=1,3
529  be(j)=pv(il,j)/pv(il,4)
530  500 CONTINUE
531  ga=pv(il,4)/pv(il,5)
532  DO 520 i=n+il,n+nd
533  bep=be(1)*p(i,1)+be(2)*p(i,2)+be(3)*p(i,3)
534  DO 510 j=1,3
535  p(i,j)=p(i,j)+ga*(ga*bep/(1d0+ga)+p(i,4))*be(j)
536  510 CONTINUE
537  p(i,4)=ga*(p(i,4)+bep)
538  520 CONTINUE
539  530 CONTINUE
540 
541 C...Check that no infinite loop in matrix element weight.
542  ntry=ntry+1
543  IF(ntry.GT.800) goto 560
544 
545 C...Matrix elements for omega and phi decays.
546  IF(mmat.EQ.1) THEN
547  wt=(p(n+1,5)*p(n+2,5)*p(n+3,5))**2-(p(n+1,5)*four(n+2,n+3))**2
548  & -(p(n+2,5)*four(n+1,n+3))**2-(p(n+3,5)*four(n+1,n+2))**2
549  & +2d0*four(n+1,n+2)*four(n+1,n+3)*four(n+2,n+3)
550  IF(max(wt*wtcor(9)/p(ip,5)**6,0.001d0).LT.pyr(0)) goto 390
551 
552 C...Matrix elements for pi0 or eta Dalitz decay to gamma e+ e-.
553  ELSEIF(mmat.EQ.2) THEN
554  four12=four(n+1,n+2)
555  four13=four(n+1,n+3)
556  wt=(pmst-0.5d0*pmes)*(four12**2+four13**2)+
557  & pmes*(four12*four13+four12**2+four13**2)
558  IF(wt.LT.pyr(0)*0.25d0*pmst*(p(ip,5)**2-pmst)**2) goto 460
559 
560 C...Matrix element for S0 -> S1 + V1 -> S1 + S2 + S3 (S scalar,
561 C...V vector), of form cos**2(theta02) in V1 rest frame, and for
562 C...S0 -> gamma + V1 -> gamma + S2 + S3, of form sin**2(theta02).
563  ELSEIF(mmat.EQ.3.AND.nm.EQ.2) THEN
564  four10=four(ip,im)
565  four12=four(ip,n+1)
566  four02=four(im,n+1)
567  pms1=p(ip,5)**2
568  pms0=p(im,5)**2
569  pms2=p(n+1,5)**2
570  IF(kfas.NE.22) hnum=(four10*four12-pms1*four02)**2
571  IF(kfas.EQ.22) hnum=pms1*(2d0*four10*four12*four02-
572  & pms1*four02**2-pms0*four12**2-pms2*four10**2+pms1*pms0*pms2)
573  hnum=max(1d-6*pms1**2*pms0*pms2,hnum)
574  hden=(four10**2-pms1*pms0)*(four12**2-pms1*pms2)
575  IF(hnum.LT.pyr(0)*hden) goto 460
576 
577 C...Matrix element for "onium" -> g + g + g or gamma + g + g.
578  ELSEIF(mmat.EQ.4) THEN
579  hx1=2d0*four(ip,n+1)/p(ip,5)**2
580  hx2=2d0*four(ip,n+2)/p(ip,5)**2
581  hx3=2d0*four(ip,n+3)/p(ip,5)**2
582  wt=((1d0-hx1)/(hx2*hx3))**2+((1d0-hx2)/(hx1*hx3))**2+
583  & ((1d0-hx3)/(hx1*hx2))**2
584  IF(wt.LT.2d0*pyr(0)) goto 390
585  IF(k(ip+1,2).EQ.22.AND.(1d0-hx1)*p(ip,5)**2.LT.4d0*parj(32)**2)
586  & goto 390
587 
588 C...Effective matrix element for nu spectrum in tau -> nu + hadrons.
589  ELSEIF(mmat.EQ.41) THEN
590  IF(mbst.EQ.0) hx1=2d0*four(ip,n+1)/p(ip,5)**2
591  IF(mbst.EQ.1) hx1=2d0*p(n+1,4)/p(ip,5)
592  hxm=min(0.75d0,2d0*(1d0-ps/p(ip,5)))
593  IF(hx1*(3d0-2d0*hx1).LT.pyr(0)*hxm*(3d0-2d0*hxm)) goto 390
594 
595 C...Matrix elements for weak decays (only semileptonic for c and b)
596  ELSEIF((mmat.EQ.42.OR.mmat.EQ.43.OR.mmat.EQ.44.OR.mmat.EQ.48)
597  & .AND.nd.EQ.3) THEN
598  IF(mbst.EQ.0) wt=four(ip,n+1)*four(n+2,n+3)
599  IF(mbst.EQ.1) wt=p(ip,5)*p(n+1,4)*four(n+2,n+3)
600  IF(wt.LT.pyr(0)*p(ip,5)*pv(1,5)**3/wtcor(10)) goto 390
601  ELSEIF(mmat.EQ.42.OR.mmat.EQ.43.OR.mmat.EQ.44.OR.mmat.EQ.48) THEN
602  DO 550 j=1,4
603  p(n+np+1,j)=0d0
604  DO 540 is=n+3,n+np
605  p(n+np+1,j)=p(n+np+1,j)+p(is,j)
606  540 CONTINUE
607  550 CONTINUE
608  IF(mbst.EQ.0) wt=four(ip,n+1)*four(n+2,n+np+1)
609  IF(mbst.EQ.1) wt=p(ip,5)*p(n+1,4)*four(n+2,n+np+1)
610  IF(wt.LT.pyr(0)*p(ip,5)*pv(1,5)**3/wtcor(10)) goto 390
611  ENDIF
612 
613 C...Scale back energy and reattach spectator.
614  560 IF(mrem.EQ.1) THEN
615  DO 570 j=1,5
616  pv(1,j)=pv(1,j)/(1d0-pqt)
617  570 CONTINUE
618  nd=nd+1
619  mrem=0
620  ENDIF
621 
622 C...Low invariant mass for system with spectator quark gives particle,
623 C...not two jets. Readjust momenta accordingly.
624  IF(mmat.EQ.31.AND.nd.EQ.3) THEN
625  mstj(93)=1
626  pm2=pymass(k(n+2,2))
627  mstj(93)=1
628  pm3=pymass(k(n+3,2))
629  IF(p(n+2,5)**2+p(n+3,5)**2+2d0*four(n+2,n+3).GE.
630  & (parj(32)+pm2+pm3)**2) goto 630
631  k(n+2,1)=1
632  kftemp=k(n+2,2)
633  CALL pykfdi(kftemp,k(n+3,2),kfldmp,k(n+2,2))
634  IF(k(n+2,2).EQ.0) goto 260
635  p(n+2,5)=pymass(k(n+2,2))
636  ps=p(n+1,5)+p(n+2,5)
637  pv(2,5)=p(n+2,5)
638  mmat=0
639  nd=2
640  goto 460
641  ELSEIF(mmat.EQ.44) THEN
642  mstj(93)=1
643  pm3=pymass(k(n+3,2))
644  mstj(93)=1
645  pm4=pymass(k(n+4,2))
646  IF(p(n+3,5)**2+p(n+4,5)**2+2d0*four(n+3,n+4).GE.
647  & (parj(32)+pm3+pm4)**2) goto 600
648  k(n+3,1)=1
649  kftemp=k(n+3,2)
650  CALL pykfdi(kftemp,k(n+4,2),kfldmp,k(n+3,2))
651  IF(k(n+3,2).EQ.0) goto 260
652  p(n+3,5)=pymass(k(n+3,2))
653  DO 580 j=1,3
654  p(n+3,j)=p(n+3,j)+p(n+4,j)
655  580 CONTINUE
656  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)
657  ha=p(n+1,4)**2-p(n+2,4)**2
658  hb=ha-(p(n+1,5)**2-p(n+2,5)**2)
659  hc=(p(n+1,1)-p(n+2,1))**2+(p(n+1,2)-p(n+2,2))**2+
660  & (p(n+1,3)-p(n+2,3))**2
661  hd=(pv(1,4)-p(n+3,4))**2
662  he=ha**2-2d0*hd*(p(n+1,4)**2+p(n+2,4)**2)+hd**2
663  hf=hd*hc-hb**2
664  hg=hd*hc-ha*hb
665  hh=(sqrt(hg**2+he*hf)-hg)/(2d0*hf)
666  DO 590 j=1,3
667  pcor=hh*(p(n+1,j)-p(n+2,j))
668  p(n+1,j)=p(n+1,j)+pcor
669  p(n+2,j)=p(n+2,j)-pcor
670  590 CONTINUE
671  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)
672  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)
673  nd=nd-1
674  ENDIF
675 
676 C...Check invariant mass of W jets. May give one particle or start over.
677  600 IF((mmat.EQ.42.OR.mmat.EQ.43.OR.mmat.EQ.44.OR.mmat.EQ.48)
678  &.AND.iabs(k(n+1,2)).LT.10) THEN
679  pmr=sqrt(max(0d0,p(n+1,5)**2+p(n+2,5)**2+2d0*four(n+1,n+2)))
680  mstj(93)=1
681  pm1=pymass(k(n+1,2))
682  mstj(93)=1
683  pm2=pymass(k(n+2,2))
684  IF(pmr.GT.parj(32)+pm1+pm2) goto 610
685  kfldum=int(1.5d0+pyr(0))
686  CALL pykfdi(k(n+1,2),-isign(kfldum,k(n+1,2)),kfldmp,kf1)
687  CALL pykfdi(k(n+2,2),-isign(kfldum,k(n+2,2)),kfldmp,kf2)
688  IF(kf1.EQ.0.OR.kf2.EQ.0) goto 260
689  psm=pymass(kf1)+pymass(kf2)
690  IF((mmat.EQ.42.OR.mmat.EQ.48).AND.pmr.GT.parj(64)+psm) goto 610
691  IF(mmat.GE.43.AND.pmr.GT.0.2d0*parj(32)+psm) goto 610
692  IF(mmat.EQ.48) goto 390
693  IF(nd.EQ.4.OR.kfa.EQ.15) goto 260
694  k(n+1,1)=1
695  kftemp=k(n+1,2)
696  CALL pykfdi(kftemp,k(n+2,2),kfldmp,k(n+1,2))
697  IF(k(n+1,2).EQ.0) goto 260
698  p(n+1,5)=pymass(k(n+1,2))
699  k(n+2,2)=k(n+3,2)
700  p(n+2,5)=p(n+3,5)
701  ps=p(n+1,5)+p(n+2,5)
702  IF(ps+parj(64).GT.pv(1,5)) goto 260
703  pv(2,5)=p(n+3,5)
704  mmat=0
705  nd=2
706  goto 460
707  ENDIF
708 
709 C...Phase space decay of partons from W decay.
710  610 IF((mmat.EQ.42.OR.mmat.EQ.48).AND.iabs(k(n+1,2)).LT.10) THEN
711  kflo(1)=k(n+1,2)
712  kflo(2)=k(n+2,2)
713  k(n+1,1)=k(n+3,1)
714  k(n+1,2)=k(n+3,2)
715  DO 620 j=1,5
716  pv(1,j)=p(n+1,j)+p(n+2,j)
717  p(n+1,j)=p(n+3,j)
718  620 CONTINUE
719  pv(1,5)=pmr
720  n=n+1
721  np=0
722  nq=2
723  ps=0d0
724  mstj(93)=2
725  psq=pymass(kflo(1))
726  mstj(93)=2
727  psq=psq+pymass(kflo(2))
728  mmat=11
729  goto 290
730  ENDIF
731 
732 C...Boost back for rapidly moving particle.
733  630 n=n+nd
734  IF(mbst.EQ.1) THEN
735  DO 640 j=1,3
736  be(j)=p(ip,j)/p(ip,4)
737  640 CONTINUE
738  ga=p(ip,4)/p(ip,5)
739  DO 660 i=nsav+1,n
740  bep=be(1)*p(i,1)+be(2)*p(i,2)+be(3)*p(i,3)
741  DO 650 j=1,3
742  p(i,j)=p(i,j)+ga*(ga*bep/(1d0+ga)+p(i,4))*be(j)
743  650 CONTINUE
744  p(i,4)=ga*(p(i,4)+bep)
745  660 CONTINUE
746  ENDIF
747 
748 C...Fill in position of decay vertex.
749  DO 680 i=nsav+1,n
750  DO 670 j=1,4
751  v(i,j)=vdcy(j)
752  670 CONTINUE
753  v(i,5)=0d0
754  680 CONTINUE
755 
756 C...Set up for parton shower evolution from jets.
757  IF(mstj(23).GE.1.AND.mmat.EQ.4.AND.k(nsav+1,2).EQ.21) THEN
758  k(nsav+1,1)=3
759  k(nsav+2,1)=3
760  k(nsav+3,1)=3
761  k(nsav+1,4)=mstu(5)*(nsav+2)
762  k(nsav+1,5)=mstu(5)*(nsav+3)
763  k(nsav+2,4)=mstu(5)*(nsav+3)
764  k(nsav+2,5)=mstu(5)*(nsav+1)
765  k(nsav+3,4)=mstu(5)*(nsav+1)
766  k(nsav+3,5)=mstu(5)*(nsav+2)
767  mstj(92)=-(nsav+1)
768  ELSEIF(mstj(23).GE.1.AND.mmat.EQ.4) THEN
769  k(nsav+2,1)=3
770  k(nsav+3,1)=3
771  k(nsav+2,4)=mstu(5)*(nsav+3)
772  k(nsav+2,5)=mstu(5)*(nsav+3)
773  k(nsav+3,4)=mstu(5)*(nsav+2)
774  k(nsav+3,5)=mstu(5)*(nsav+2)
775  mstj(92)=nsav+2
776  ELSEIF(mstj(23).GE.1.AND.(mmat.EQ.32.OR.mmat.EQ.44).AND.
777  & iabs(k(nsav+1,2)).LE.10.AND.iabs(k(nsav+2,2)).LE.10) THEN
778  k(nsav+1,1)=3
779  k(nsav+2,1)=3
780  k(nsav+1,4)=mstu(5)*(nsav+2)
781  k(nsav+1,5)=mstu(5)*(nsav+2)
782  k(nsav+2,4)=mstu(5)*(nsav+1)
783  k(nsav+2,5)=mstu(5)*(nsav+1)
784  mstj(92)=nsav+1
785  ELSEIF(mstj(23).GE.1.AND.(mmat.EQ.32.OR.mmat.EQ.44).AND.
786  & iabs(k(nsav+1,2)).LE.20.AND.iabs(k(nsav+2,2)).LE.20) THEN
787  mstj(92)=nsav+1
788  ELSEIF(mstj(23).GE.1.AND.mmat.EQ.33.AND.iabs(k(nsav+2,2)).EQ.21)
789  & THEN
790  k(nsav+1,1)=3
791  k(nsav+2,1)=3
792  k(nsav+3,1)=3
793  kcp=pycomp(k(nsav+1,2))
794  kqp=kchg(kcp,2)*isign(1,k(nsav+1,2))
795  jcon=4
796  IF(kqp.LT.0) jcon=5
797  k(nsav+1,jcon)=mstu(5)*(nsav+2)
798  k(nsav+2,9-jcon)=mstu(5)*(nsav+1)
799  k(nsav+2,jcon)=mstu(5)*(nsav+3)
800  k(nsav+3,9-jcon)=mstu(5)*(nsav+2)
801  mstj(92)=nsav+1
802  ELSEIF(mstj(23).GE.1.AND.mmat.EQ.33) THEN
803  k(nsav+1,1)=3
804  k(nsav+3,1)=3
805  k(nsav+1,4)=mstu(5)*(nsav+3)
806  k(nsav+1,5)=mstu(5)*(nsav+3)
807  k(nsav+3,4)=mstu(5)*(nsav+1)
808  k(nsav+3,5)=mstu(5)*(nsav+1)
809  mstj(92)=nsav+1
810  ENDIF
811 
812 C...Mark decayed particle; special option for B-Bbar mixing.
813  IF(k(ip,1).EQ.5) k(ip,1)=15
814  IF(k(ip,1).LE.10) k(ip,1)=11
815  IF(mmix.EQ.1.AND.mstj(26).EQ.2.AND.k(ip,1).EQ.11) k(ip,1)=12
816  k(ip,4)=nsav+1
817  k(ip,5)=n
818 
819  RETURN
820  END