Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ludecy.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ludecy.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE ludecy(IP)
5 
6 C...Purpose: to handle the decay of unstable particles.
7  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
8  SAVE /lujets/
9  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
10  SAVE /ludat1/
11  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
12  SAVE /ludat2/
13  common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)
14  SAVE /ludat3/
15  dimension vdcy(4),kflo(4),kfl1(4),pv(10,5),rord(10),ue(3),be(3),
16  &wtcor(10)
17  DATA wtcor/2.,5.,15.,60.,250.,1500.,1.2e4,1.2e5,150.,16./
18 
19 C...Functions: momentum in two-particle decays, four-product and
20 C...matrix element times phase space in weak decays.
21  pawt(a,b,c)=sqrt((a**2-(b+c)**2)*(a**2-(b-c)**2))/(2.*a)
22  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)
23  hmeps(ha)=((1.-hrq-ha)**2+3.*ha*(1.+hrq-ha))*
24  &sqrt((1.-hrq-ha)**2-4.*hrq*ha)
25 
26 C...Initial values.
27  ntry=0
28  nsav=n
29  kfa=iabs(k(ip,2))
30  kfs=isign(1,k(ip,2))
31  kc=lucomp(kfa)
32  mstj(92)=0
33 
34 C...Choose lifetime and determine decay vertex.
35  IF(k(ip,1).EQ.5) THEN
36  v(ip,5)=0.
37  ELSEIF(k(ip,1).NE.4) THEN
38  v(ip,5)=-pmas(kc,4)*log(rlu(0))
39  ENDIF
40  DO 100 j=1,4
41  100 vdcy(j)=v(ip,j)+v(ip,5)*p(ip,j)/p(ip,5)
42 
43 C...Determine whether decay allowed or not.
44  mout=0
45  IF(mstj(22).EQ.2) THEN
46  IF(pmas(kc,4).GT.parj(71)) mout=1
47  ELSEIF(mstj(22).EQ.3) THEN
48  IF(vdcy(1)**2+vdcy(2)**2+vdcy(3)**2.GT.parj(72)**2) mout=1
49  ELSEIF(mstj(22).EQ.4) THEN
50  IF(vdcy(1)**2+vdcy(2)**2.GT.parj(73)**2) mout=1
51  IF(abs(vdcy(3)).GT.parj(74)) mout=1
52  ENDIF
53  IF(mout.EQ.1.AND.k(ip,1).NE.5) THEN
54  k(ip,1)=4
55  RETURN
56  ENDIF
57 
58 C...Check existence of decay channels. Particle/antiparticle rules.
59  kca=kc
60  IF(mdcy(kc,2).GT.0) THEN
61  mdmdcy=mdme(mdcy(kc,2),2)
62  IF(mdmdcy.GT.80.AND.mdmdcy.LE.90) kca=mdmdcy
63  ENDIF
64  IF(mdcy(kca,2).LE.0.OR.mdcy(kca,3).LE.0) THEN
65  CALL luerrm(9,'(LUDECY:) no decay channel defined')
66  RETURN
67  ENDIF
68  IF(mod(kfa/1000,10).EQ.0.AND.(kca.EQ.85.OR.kca.EQ.87)) kfs=-kfs
69  IF(kchg(kc,3).EQ.0) THEN
70  kfsp=1
71  kfsn=0
72  IF(rlu(0).GT.0.5) kfs=-kfs
73  ELSEIF(kfs.GT.0) THEN
74  kfsp=1
75  kfsn=0
76  ELSE
77  kfsp=0
78  kfsn=1
79  ENDIF
80 
81 C...Sum branching ratios of allowed decay channels.
82  110 nope=0
83  brsu=0.
84  DO 120 idl=mdcy(kca,2),mdcy(kca,2)+mdcy(kca,3)-1
85  IF(mdme(idl,1).NE.1.AND.kfsp*mdme(idl,1).NE.2.AND.
86  &kfsn*mdme(idl,1).NE.3) goto 120
87  IF(mdme(idl,2).GT.100) goto 120
88  nope=nope+1
89  brsu=brsu+brat(idl)
90  120 CONTINUE
91  IF(nope.EQ.0) THEN
92  CALL luerrm(2,'(LUDECY:) all decay channels closed by user')
93  RETURN
94  ENDIF
95 
96 C...Select decay channel among allowed ones.
97  130 rbr=brsu*rlu(0)
98  idl=mdcy(kca,2)-1
99  140 idl=idl+1
100  IF(mdme(idl,1).NE.1.AND.kfsp*mdme(idl,1).NE.2.AND.
101  &kfsn*mdme(idl,1).NE.3) THEN
102  IF(idl.LT.mdcy(kca,2)+mdcy(kca,3)-1) goto 140
103  ELSEIF(mdme(idl,2).GT.100) THEN
104  IF(idl.LT.mdcy(kca,2)+mdcy(kca,3)-1) goto 140
105  ELSE
106  idc=idl
107  rbr=rbr-brat(idl)
108  IF(idl.LT.mdcy(kca,2)+mdcy(kca,3)-1.AND.rbr.GT.0.) goto 140
109  ENDIF
110 
111 C...Start readout of decay channel: matrix element, reset counters.
112  mmat=mdme(idc,2)
113  150 ntry=ntry+1
114  IF(ntry.GT.1000) THEN
115  CALL luerrm(14,'(LUDECY:) 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(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(-kfs*int(1.+(2.+parj(2))*rlu(0)),0,kfp,kdump)
153  IF(kfp.EQ.0) goto 150
154  mstj(93)=1
155  IF(pv(1,5).LT.parj(32)+2.*ulmass(kfp)) goto 150
156  ELSEIF(kp.EQ.-82) THEN
157  kfp=-kfp
158  IF(iabs(kfp).GT.10) kfp=kfp+isign(10000,kfp)
159  ENDIF
160  IF(kpa.EQ.81.OR.kpa.EQ.82) kcp=lucomp(kfp)
161 
162 C...Add decay product to event record or to quark flavour list.
163  kfpa=iabs(kfp)
164  kqp=kchg(kcp,2)
165  IF(mmat.GE.11.AND.mmat.LE.30.AND.kqp.NE.0) THEN
166  nq=nq+1
167  kflo(nq)=kfp
168  mstj(93)=2
169  psq=psq+ulmass(kflo(nq))
170  ELSEIF(mmat.GE.42.AND.mmat.LE.43.AND.np.EQ.3.AND.mod(nq,2).EQ.1)
171  &THEN
172  nq=nq-1
173  ps=ps-p(i,5)
174  k(i,1)=1
175  kfi=k(i,2)
176  CALL lukfdi(kfp,kfi,kfldmp,k(i,2))
177  IF(k(i,2).EQ.0) goto 150
178  mstj(93)=1
179  p(i,5)=ulmass(k(i,2))
180  ps=ps+p(i,5)
181  ELSE
182  i=i+1
183  np=np+1
184  IF(mmat.NE.33.AND.kqp.NE.0) nq=nq+1
185  IF(mmat.EQ.33.AND.kqp.NE.0.AND.kqp.NE.2) nq=nq+1
186  k(i,1)=1+mod(nq,2)
187  IF(mmat.EQ.4.AND.jt.LE.2.AND.kfp.EQ.21) k(i,1)=2
188  IF(mmat.EQ.4.AND.jt.EQ.3) k(i,1)=1
189  k(i,2)=kfp
190  k(i,3)=ip
191  k(i,4)=0
192  k(i,5)=0
193  p(i,5)=ulmass(kfp)
194  IF(mmat.EQ.45.AND.kfpa.EQ.89) p(i,5)=parj(32)
195  ps=ps+p(i,5)
196  ENDIF
197  170 CONTINUE
198 
199 C...Choose decay multiplicity in phase space model.
200  180 IF(mmat.GE.11.AND.mmat.LE.30) THEN
201  psp=ps
202  cnde=parj(61)*log(max((pv(1,5)-ps-psq)/parj(62),1.1))
203  IF(mmat.EQ.12) cnde=cnde+parj(63)
204  190 ntry=ntry+1
205  IF(ntry.GT.1000) THEN
206  CALL luerrm(14,'(LUDECY:) caught in infinite loop')
207  IF(mstu(21).GE.1) RETURN
208  ENDIF
209  IF(mmat.LE.20) THEN
210  gauss=sqrt(-2.*cnde*log(max(1e-10,rlu(0))))*
211  & sin(paru(2)*rlu(0))
212  nd=0.5+0.5*np+0.25*nq+cnde+gauss
213  IF(nd.LT.np+nq/2.OR.nd.LT.2.OR.nd.GT.10) goto 190
214  IF(mmat.EQ.13.AND.nd.EQ.2) goto 190
215  IF(mmat.EQ.14.AND.nd.LE.3) goto 190
216  IF(mmat.EQ.15.AND.nd.LE.4) goto 190
217  ELSE
218  nd=mmat-20
219  ENDIF
220 
221 C...Form hadrons from flavour content.
222  DO 200 jt=1,4
223  200 kfl1(jt)=kflo(jt)
224  IF(nd.EQ.np+nq/2) goto 220
225  DO 210 i=n+np+1,n+nd-nq/2
226  jt=1+int((nq-1)*rlu(0))
227  CALL lukfdi(kfl1(jt),0,kfl2,k(i,2))
228  IF(k(i,2).EQ.0) goto 190
229  210 kfl1(jt)=-kfl2
230  220 jt=2
231  jt2=3
232  jt3=4
233  IF(nq.EQ.4.AND.rlu(0).LT.parj(66)) jt=4
234  IF(jt.EQ.4.AND.isign(1,kfl1(1)*(10-iabs(kfl1(1))))*
235  & isign(1,kfl1(jt)*(10-iabs(kfl1(jt)))).GT.0) jt=3
236  IF(jt.EQ.3) jt2=2
237  IF(jt.EQ.4) jt3=2
238  CALL lukfdi(kfl1(1),kfl1(jt),kfldmp,k(n+nd-nq/2+1,2))
239  IF(k(n+nd-nq/2+1,2).EQ.0) goto 190
240  IF(nq.EQ.4) CALL lukfdi(kfl1(jt2),kfl1(jt3),kfldmp,k(n+nd,2))
241  IF(nq.EQ.4.AND.k(n+nd,2).EQ.0) goto 190
242 
243 C...Check that sum of decay product masses not too large.
244  ps=psp
245  DO 230 i=n+np+1,n+nd
246  k(i,1)=1
247  k(i,3)=ip
248  k(i,4)=0
249  k(i,5)=0
250  p(i,5)=ulmass(k(i,2))
251  230 ps=ps+p(i,5)
252  IF(ps+parj(64).GT.pv(1,5)) goto 190
253 
254 C...Rescale energy to subtract off spectator quark mass.
255  ELSEIF((mmat.EQ.31.OR.mmat.EQ.33.OR.mmat.EQ.44.OR.mmat.EQ.45).
256  &and.np.GE.3) THEN
257  ps=ps-p(n+np,5)
258  pqt=(p(n+np,5)+parj(65))/pv(1,5)
259  DO 240 j=1,5
260  p(n+np,j)=pqt*pv(1,j)
261  240 pv(1,j)=(1.-pqt)*pv(1,j)
262  IF(ps+parj(64).GT.pv(1,5)) goto 150
263  nd=np-1
264  mrem=1
265 
266 C...Phase space factors imposed in W decay.
267  ELSEIF(mmat.EQ.46) THEN
268  mstj(93)=1
269  psmc=ulmass(k(n+1,2))
270  mstj(93)=1
271  psmc=psmc+ulmass(k(n+2,2))
272  IF(max(ps,psmc)+parj(32).GT.pv(1,5)) goto 130
273  hr1=(p(n+1,5)/pv(1,5))**2
274  hr2=(p(n+2,5)/pv(1,5))**2
275  IF((1.-hr1-hr2)*(2.+hr1+hr2)*sqrt((1.-hr1-hr2)**2-4.*hr1*hr2).
276  & lt.2.*rlu(0)) goto 130
277  nd=np
278 
279 C...Fully specified final state: check mass broadening effects.
280  ELSE
281  IF(np.GE.2.AND.ps+parj(64).GT.pv(1,5)) goto 150
282  nd=np
283  ENDIF
284 
285 C...Select W mass in decay Q -> W + q, without W propagator.
286  IF(mmat.EQ.45.AND.mstj(25).LE.0) THEN
287  hlq=(parj(32)/pv(1,5))**2
288  huq=(1.-(p(n+2,5)+parj(64))/pv(1,5))**2
289  hrq=(p(n+2,5)/pv(1,5))**2
290  250 hw=hlq+rlu(0)*(huq-hlq)
291  IF(hmeps(hw).LT.rlu(0)) goto 250
292  p(n+1,5)=pv(1,5)*sqrt(hw)
293 
294 C...Ditto, including W propagator. Divide mass range into three regions.
295  ELSEIF(mmat.EQ.45) THEN
296  hqw=(pv(1,5)/pmas(24,1))**2
297  hlw=(parj(32)/pmas(24,1))**2
298  huw=((pv(1,5)-p(n+2,5)-parj(64))/pmas(24,1))**2
299  hrq=(p(n+2,5)/pv(1,5))**2
300  hg=pmas(24,2)/pmas(24,1)
301  hatl=atan((hlw-1.)/hg)
302  hm=min(1.,huw-0.001)
303  hmv1=hmeps(hm/hqw)/((hm-1.)**2+hg**2)
304  260 hm=hm-hg
305  hmv2=hmeps(hm/hqw)/((hm-1.)**2+hg**2)
306  hsav1=hmeps(hm/hqw)
307  hsav2=1./((hm-1.)**2+hg**2)
308  IF(hmv2.GT.hmv1.AND.hm-hg.GT.hlw) THEN
309  hmv1=hmv2
310  goto 260
311  ENDIF
312  hmv=min(2.*hmv1,hmeps(hm/hqw)/hg**2)
313  hm1=1.-sqrt(1./hmv-hg**2)
314  IF(hm1.GT.hlw.AND.hm1.LT.hm) THEN
315  hm=hm1
316  ELSEIF(hmv2.LE.hmv1) THEN
317  hm=max(hlw,hm-min(0.1,1.-hm))
318  ENDIF
319  hatm=atan((hm-1.)/hg)
320  hwt1=(hatm-hatl)/hg
321  hwt2=hmv*(min(1.,huw)-hm)
322  hwt3=0.
323  IF(huw.GT.1.) THEN
324  hatu=atan((huw-1.)/hg)
325  hmp1=hmeps(1./hqw)
326  hwt3=hmp1*hatu/hg
327  ENDIF
328 
329 C...Select mass region and W mass there. Accept according to weight.
330  270 hreg=rlu(0)*(hwt1+hwt2+hwt3)
331  IF(hreg.LE.hwt1) THEN
332  hw=1.+hg*tan(hatl+rlu(0)*(hatm-hatl))
333  hacc=hmeps(hw/hqw)
334  ELSEIF(hreg.LE.hwt1+hwt2) THEN
335  hw=hm+rlu(0)*(min(1.,huw)-hm)
336  hacc=hmeps(hw/hqw)/((hw-1.)**2+hg**2)/hmv
337  ELSE
338  hw=1.+hg*tan(rlu(0)*hatu)
339  hacc=hmeps(hw/hqw)/hmp1
340  ENDIF
341  IF(hacc.LT.rlu(0)) goto 270
342  p(n+1,5)=pmas(24,1)*sqrt(hw)
343  ENDIF
344 
345 C...Determine position of grandmother, number of sisters, Q -> W sign.
346  nm=0
347  msgn=0
348  IF(mmat.EQ.3.OR.mmat.EQ.46) THEN
349  im=k(ip,3)
350  IF(im.LT.0.OR.im.GE.ip) im=0
351  IF(im.NE.0) kfam=iabs(k(im,2))
352  IF(im.NE.0.AND.mmat.EQ.3) THEN
353  DO 280 il=max(ip-2,im+1),min(ip+2,n)
354  280 IF(k(il,3).EQ.im) nm=nm+1
355  IF(nm.NE.2.OR.kfam.LE.100.OR.mod(kfam,10).NE.1.OR.
356  & mod(kfam/1000,10).NE.0) nm=0
357  ELSEIF(im.NE.0.AND.mmat.EQ.46) THEN
358  msgn=isign(1,k(im,2)*k(ip,2))
359  IF(kfam.GT.100.AND.mod(kfam/1000,10).EQ.0) msgn=
360  & msgn*(-1)**mod(kfam/100,10)
361  ENDIF
362  ENDIF
363 
364 C...Kinematics of one-particle decays.
365  IF(nd.EQ.1) THEN
366  DO 290 j=1,4
367  290 p(n+1,j)=p(ip,j)
368  goto 510
369  ENDIF
370 
371 C...Calculate maximum weight ND-particle decay.
372  pv(nd,5)=p(n+nd,5)
373  IF(nd.GE.3) THEN
374  wtmax=1./wtcor(nd-2)
375  pmax=pv(1,5)-ps+p(n+nd,5)
376  pmin=0.
377  DO 300 il=nd-1,1,-1
378  pmax=pmax+p(n+il,5)
379  pmin=pmin+p(n+il+1,5)
380  300 wtmax=wtmax*pawt(pmax,pmin,p(n+il,5))
381  ENDIF
382 
383 C...Find virtual gamma mass in Dalitz decay.
384  310 IF(nd.EQ.2) THEN
385  ELSEIF(mmat.EQ.2) THEN
386  pmes=4.*pmas(11,1)**2
387  pmrho2=pmas(131,1)**2
388  pgrho2=pmas(131,2)**2
389  320 pmst=pmes*(p(ip,5)**2/pmes)**rlu(0)
390  wt=(1+0.5*pmes/pmst)*sqrt(max(0.,1.-pmes/pmst))*
391  & (1.-pmst/p(ip,5)**2)**3*(1.+pgrho2/pmrho2)/
392  & ((1.-pmst/pmrho2)**2+pgrho2/pmrho2)
393  IF(wt.LT.rlu(0)) goto 320
394  pv(2,5)=max(2.00001*pmas(11,1),sqrt(pmst))
395 
396 C...M-generator gives weight. If rejected, try again.
397  ELSE
398  330 rord(1)=1.
399  DO 350 il1=2,nd-1
400  rsav=rlu(0)
401  DO 340 il2=il1-1,1,-1
402  IF(rsav.LE.rord(il2)) goto 350
403  340 rord(il2+1)=rord(il2)
404  350 rord(il2+1)=rsav
405  rord(nd)=0.
406  wt=1.
407  DO 360 il=nd-1,1,-1
408  pv(il,5)=pv(il+1,5)+p(n+il,5)+(rord(il)-rord(il+1))*(pv(1,5)-ps)
409  360 wt=wt*pawt(pv(il,5),pv(il+1,5),p(n+il,5))
410  IF(wt.LT.rlu(0)*wtmax) goto 330
411  ENDIF
412 
413 C...Perform two-particle decays in respective CM frame.
414  370 DO 390 il=1,nd-1
415  pa=pawt(pv(il,5),pv(il+1,5),p(n+il,5))
416  ue(3)=2.*rlu(0)-1.
417  phi=paru(2)*rlu(0)
418  ue(1)=sqrt(1.-ue(3)**2)*cos(phi)
419  ue(2)=sqrt(1.-ue(3)**2)*sin(phi)
420  DO 380 j=1,3
421  p(n+il,j)=pa*ue(j)
422  380 pv(il+1,j)=-pa*ue(j)
423  p(n+il,4)=sqrt(pa**2+p(n+il,5)**2)
424  390 pv(il+1,4)=sqrt(pa**2+pv(il+1,5)**2)
425 
426 C...Lorentz transform decay products to lab frame.
427  DO 400 j=1,4
428  400 p(n+nd,j)=pv(nd,j)
429  DO 430 il=nd-1,1,-1
430  DO 410 j=1,3
431  410 be(j)=pv(il,j)/pv(il,4)
432  ga=pv(il,4)/pv(il,5)
433  DO 430 i=n+il,n+nd
434  bep=be(1)*p(i,1)+be(2)*p(i,2)+be(3)*p(i,3)
435  DO 420 j=1,3
436  420 p(i,j)=p(i,j)+ga*(ga*bep/(1.+ga)+p(i,4))*be(j)
437  430 p(i,4)=ga*(p(i,4)+bep)
438 
439 C...Matrix elements for omega and phi decays.
440  IF(mmat.EQ.1) THEN
441  wt=(p(n+1,5)*p(n+2,5)*p(n+3,5))**2-(p(n+1,5)*four(n+2,n+3))**2
442  & -(p(n+2,5)*four(n+1,n+3))**2-(p(n+3,5)*four(n+1,n+2))**2
443  & +2.*four(n+1,n+2)*four(n+1,n+3)*four(n+2,n+3)
444  IF(max(wt*wtcor(9)/p(ip,5)**6,0.001).LT.rlu(0)) goto 310
445 
446 C...Matrix elements for pi0 or eta Dalitz decay to gamma e+ e-.
447  ELSEIF(mmat.EQ.2) THEN
448  four12=four(n+1,n+2)
449  four13=four(n+1,n+3)
450  four23=0.5*pmst-0.25*pmes
451  wt=(pmst-0.5*pmes)*(four12**2+four13**2)+
452  & pmes*(four12*four13+four12**2+four13**2)
453  IF(wt.LT.rlu(0)*0.25*pmst*(p(ip,5)**2-pmst)**2) goto 370
454 
455 C...Matrix element for S0 -> S1 + V1 -> S1 + S2 + S3 (S scalar,
456 C...V vector), of form cos**2(theta02) in V1 rest frame.
457  ELSEIF(mmat.EQ.3.AND.nm.EQ.2) THEN
458  IF((p(ip,5)**2*four(im,n+1)-four(ip,im)*four(ip,n+1))**2.LE.
459  & rlu(0)*(four(ip,im)**2-(p(ip,5)*p(im,5))**2)*(four(ip,n+1)**2-
460  & (p(ip,5)*p(n+1,5))**2)) goto 370
461 
462 C...Matrix element for "onium" -> g + g + g or gamma + g + g.
463  ELSEIF(mmat.EQ.4) THEN
464  hx1=2.*four(ip,n+1)/p(ip,5)**2
465  hx2=2.*four(ip,n+2)/p(ip,5)**2
466  hx3=2.*four(ip,n+3)/p(ip,5)**2
467  wt=((1.-hx1)/(hx2*hx3))**2+((1.-hx2)/(hx1*hx3))**2+
468  & ((1.-hx3)/(hx1*hx2))**2
469  IF(wt.LT.2.*rlu(0)) goto 310
470  IF(k(ip+1,2).EQ.22.AND.(1.-hx1)*p(ip,5)**2.LT.4.*parj(32)**2)
471  & goto 310
472 
473 C...Effective matrix element for nu spectrum in tau -> nu + hadrons.
474  ELSEIF(mmat.EQ.41) THEN
475  hx1=2.*four(ip,n+1)/p(ip,5)**2
476  IF(8.*hx1*(3.-2.*hx1)/9..LT.rlu(0)) goto 310
477 
478 C...Matrix elements for weak decays (only semileptonic for c and b)
479  ELSEIF(mmat.GE.42.AND.mmat.LE.44.AND.nd.EQ.3) THEN
480  IF(mbst.EQ.0) wt=four(ip,n+1)*four(n+2,n+3)
481  IF(mbst.EQ.1) wt=p(ip,5)*p(n+1,4)*four(n+2,n+3)
482  IF(wt.LT.rlu(0)*p(ip,5)*pv(1,5)**3/wtcor(10)) goto 310
483  ELSEIF(mmat.GE.42.AND.mmat.LE.44) THEN
484  DO 440 j=1,4
485  p(n+np+1,j)=0.
486  DO 440 is=n+3,n+np
487  440 p(n+np+1,j)=p(n+np+1,j)+p(is,j)
488  IF(mbst.EQ.0) wt=four(ip,n+1)*four(n+2,n+np+1)
489  IF(mbst.EQ.1) wt=p(ip,5)*p(n+1,4)*four(n+2,n+np+1)
490  IF(wt.LT.rlu(0)*p(ip,5)*pv(1,5)**3/wtcor(10)) goto 310
491 
492 C...Angular distribution in W decay.
493  ELSEIF(mmat.EQ.46.AND.msgn.NE.0) THEN
494  IF(msgn.GT.0) wt=four(im,n+1)*four(n+2,ip+1)
495  IF(msgn.LT.0) wt=four(im,n+2)*four(n+1,ip+1)
496  IF(wt.LT.rlu(0)*p(im,5)**4/wtcor(10)) goto 370
497  ENDIF
498 
499 C...Scale back energy and reattach spectator.
500  IF(mrem.EQ.1) THEN
501  DO 450 j=1,5
502  450 pv(1,j)=pv(1,j)/(1.-pqt)
503  nd=nd+1
504  mrem=0
505  ENDIF
506 
507 C...Low invariant mass for system with spectator quark gives particle,
508 C...not two jets. Readjust momenta accordingly.
509  IF((mmat.EQ.31.OR.mmat.EQ.45).AND.nd.EQ.3) THEN
510  mstj(93)=1
511  pm2=ulmass(k(n+2,2))
512  mstj(93)=1
513  pm3=ulmass(k(n+3,2))
514  IF(p(n+2,5)**2+p(n+3,5)**2+2.*four(n+2,n+3).GE.
515  & (parj(32)+pm2+pm3)**2) goto 510
516  k(n+2,1)=1
517  kftemp=k(n+2,2)
518  CALL lukfdi(kftemp,k(n+3,2),kfldmp,k(n+2,2))
519  IF(k(n+2,2).EQ.0) goto 150
520  p(n+2,5)=ulmass(k(n+2,2))
521  ps=p(n+1,5)+p(n+2,5)
522  pv(2,5)=p(n+2,5)
523  mmat=0
524  nd=2
525  goto 370
526  ELSEIF(mmat.EQ.44) THEN
527  mstj(93)=1
528  pm3=ulmass(k(n+3,2))
529  mstj(93)=1
530  pm4=ulmass(k(n+4,2))
531  IF(p(n+3,5)**2+p(n+4,5)**2+2.*four(n+3,n+4).GE.
532  & (parj(32)+pm3+pm4)**2) goto 480
533  k(n+3,1)=1
534  kftemp=k(n+3,2)
535  CALL lukfdi(kftemp,k(n+4,2),kfldmp,k(n+3,2))
536  IF(k(n+3,2).EQ.0) goto 150
537  p(n+3,5)=ulmass(k(n+3,2))
538  DO 460 j=1,3
539  460 p(n+3,j)=p(n+3,j)+p(n+4,j)
540  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)
541  ha=p(n+1,4)**2-p(n+2,4)**2
542  hb=ha-(p(n+1,5)**2-p(n+2,5)**2)
543  hc=(p(n+1,1)-p(n+2,1))**2+(p(n+1,2)-p(n+2,2))**2+
544  & (p(n+1,3)-p(n+2,3))**2
545  hd=(pv(1,4)-p(n+3,4))**2
546  he=ha**2-2.*hd*(p(n+1,4)**2+p(n+2,4)**2)+hd**2
547  hf=hd*hc-hb**2
548  hg=hd*hc-ha*hb
549  hh=(sqrt(hg**2+he*hf)-hg)/(2.*hf)
550  DO 470 j=1,3
551  pcor=hh*(p(n+1,j)-p(n+2,j))
552  p(n+1,j)=p(n+1,j)+pcor
553  470 p(n+2,j)=p(n+2,j)-pcor
554  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)
555  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)
556  nd=nd-1
557  ENDIF
558 
559 C...Check invariant mass of W jets. May give one particle or start over.
560  480 IF(mmat.GE.42.AND.mmat.LE.44.AND.iabs(k(n+1,2)).LT.10) THEN
561  pmr=sqrt(max(0.,p(n+1,5)**2+p(n+2,5)**2+2.*four(n+1,n+2)))
562  mstj(93)=1
563  pm1=ulmass(k(n+1,2))
564  mstj(93)=1
565  pm2=ulmass(k(n+2,2))
566  IF(pmr.GT.parj(32)+pm1+pm2) goto 490
567  kfldum=int(1.5+rlu(0))
568  CALL lukfdi(k(n+1,2),-isign(kfldum,k(n+1,2)),kfldmp,kf1)
569  CALL lukfdi(k(n+2,2),-isign(kfldum,k(n+2,2)),kfldmp,kf2)
570  IF(kf1.EQ.0.OR.kf2.EQ.0) goto 150
571  psm=ulmass(kf1)+ulmass(kf2)
572  IF(mmat.EQ.42.AND.pmr.GT.parj(64)+psm) goto 490
573  IF(mmat.GE.43.AND.pmr.GT.0.2*parj(32)+psm) goto 490
574  IF(nd.EQ.4.OR.kfa.EQ.15) goto 150
575  k(n+1,1)=1
576  kftemp=k(n+1,2)
577  CALL lukfdi(kftemp,k(n+2,2),kfldmp,k(n+1,2))
578  IF(k(n+1,2).EQ.0) goto 150
579  p(n+1,5)=ulmass(k(n+1,2))
580  k(n+2,2)=k(n+3,2)
581  p(n+2,5)=p(n+3,5)
582  ps=p(n+1,5)+p(n+2,5)
583  pv(2,5)=p(n+3,5)
584  mmat=0
585  nd=2
586  goto 370
587  ENDIF
588 
589 C...Phase space decay of partons from W decay.
590  490 IF(mmat.EQ.42.AND.iabs(k(n+1,2)).LT.10) THEN
591  kflo(1)=k(n+1,2)
592  kflo(2)=k(n+2,2)
593  k(n+1,1)=k(n+3,1)
594  k(n+1,2)=k(n+3,2)
595  DO 500 j=1,5
596  pv(1,j)=p(n+1,j)+p(n+2,j)
597  500 p(n+1,j)=p(n+3,j)
598  pv(1,5)=pmr
599  n=n+1
600  np=0
601  nq=2
602  ps=0.
603  mstj(93)=2
604  psq=ulmass(kflo(1))
605  mstj(93)=2
606  psq=psq+ulmass(kflo(2))
607  mmat=11
608  goto 180
609  ENDIF
610 
611 C...Boost back for rapidly moving particle.
612  510 n=n+nd
613  IF(mbst.EQ.1) THEN
614  DO 520 j=1,3
615  520 be(j)=p(ip,j)/p(ip,4)
616  ga=p(ip,4)/p(ip,5)
617  DO 540 i=nsav+1,n
618  bep=be(1)*p(i,1)+be(2)*p(i,2)+be(3)*p(i,3)
619  DO 530 j=1,3
620  530 p(i,j)=p(i,j)+ga*(ga*bep/(1.+ga)+p(i,4))*be(j)
621  540 p(i,4)=ga*(p(i,4)+bep)
622  ENDIF
623 
624 C...Fill in position of decay vertex.
625  DO 560 i=nsav+1,n
626  DO 550 j=1,4
627  550 v(i,j)=vdcy(j)
628  560 v(i,5)=0.
629 
630 C...Set up for parton shower evolution from jets.
631  IF(mstj(23).GE.1.AND.mmat.EQ.4.AND.k(nsav+1,2).EQ.21) THEN
632  k(nsav+1,1)=3
633  k(nsav+2,1)=3
634  k(nsav+3,1)=3
635  k(nsav+1,4)=mstu(5)*(nsav+2)
636  k(nsav+1,5)=mstu(5)*(nsav+3)
637  k(nsav+2,4)=mstu(5)*(nsav+3)
638  k(nsav+2,5)=mstu(5)*(nsav+1)
639  k(nsav+3,4)=mstu(5)*(nsav+1)
640  k(nsav+3,5)=mstu(5)*(nsav+2)
641  mstj(92)=-(nsav+1)
642  ELSEIF(mstj(23).GE.1.AND.mmat.EQ.4) THEN
643  k(nsav+2,1)=3
644  k(nsav+3,1)=3
645  k(nsav+2,4)=mstu(5)*(nsav+3)
646  k(nsav+2,5)=mstu(5)*(nsav+3)
647  k(nsav+3,4)=mstu(5)*(nsav+2)
648  k(nsav+3,5)=mstu(5)*(nsav+2)
649  mstj(92)=nsav+2
650  ELSEIF(mstj(23).GE.1.AND.(mmat.EQ.32.OR.mmat.EQ.44.OR.mmat.EQ.46).
651  &and.iabs(k(nsav+1,2)).LE.10.AND.iabs(k(nsav+2,2)).LE.10) THEN
652  k(nsav+1,1)=3
653  k(nsav+2,1)=3
654  k(nsav+1,4)=mstu(5)*(nsav+2)
655  k(nsav+1,5)=mstu(5)*(nsav+2)
656  k(nsav+2,4)=mstu(5)*(nsav+1)
657  k(nsav+2,5)=mstu(5)*(nsav+1)
658  mstj(92)=nsav+1
659  ELSEIF(mstj(23).GE.1.AND.mmat.EQ.33.AND.iabs(k(nsav+2,2)).EQ.21)
660  &THEN
661  k(nsav+1,1)=3
662  k(nsav+2,1)=3
663  k(nsav+3,1)=3
664  kcp=lucomp(k(nsav+1,2))
665  kqp=kchg(kcp,2)*isign(1,k(nsav+1,2))
666  jcon=4
667  IF(kqp.LT.0) jcon=5
668  k(nsav+1,jcon)=mstu(5)*(nsav+2)
669  k(nsav+2,9-jcon)=mstu(5)*(nsav+1)
670  k(nsav+2,jcon)=mstu(5)*(nsav+3)
671  k(nsav+3,9-jcon)=mstu(5)*(nsav+2)
672  mstj(92)=nsav+1
673  ELSEIF(mstj(23).GE.1.AND.mmat.EQ.33) THEN
674  k(nsav+1,1)=3
675  k(nsav+3,1)=3
676  k(nsav+1,4)=mstu(5)*(nsav+3)
677  k(nsav+1,5)=mstu(5)*(nsav+3)
678  k(nsav+3,4)=mstu(5)*(nsav+1)
679  k(nsav+3,5)=mstu(5)*(nsav+1)
680  mstj(92)=nsav+1
681  ENDIF
682 
683 C...Mark decayed particle.
684  IF(k(ip,1).EQ.5) k(ip,1)=15
685  IF(k(ip,1).LE.10) k(ip,1)=11
686  k(ip,4)=nsav+1
687  k(ip,5)=n
688 
689  RETURN
690  END