7 common/lujets/
n,
k(9000,5),
p(9000,5),
v(9000,5)
9 common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
11 common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
13 common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)
15 dimension vdcy(4),kflo(4),kfl1(4),pv(10,5),rord(10),ue(3),be(3),
17 DATA wtcor/2.,5.,15.,60.,250.,1500.,1.2e4,1.2e5,150.,16./
21 pawt(a,b,
c)=sqrt((a**2-(b+
c)**2)*(a**2-(b-
c)**2))/(2.*a)
23 hmeps(ha)=((1.-hrq-ha)**2+3.*ha*(1.+hrq-ha))*
24 &sqrt((1.-hrq-ha)**2-4.*hrq*ha)
37 ELSEIF(
k(ip,1).NE.4)
THEN
41 100 vdcy(
j)=
v(ip,
j)+
v(ip,5)*
p(ip,
j)/
p(ip,5)
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
53 IF(mout.EQ.1.AND.
k(ip,1).NE.5)
THEN
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
64 IF(mdcy(kca,2).LE.0.OR.mdcy(kca,3).LE.0)
THEN
65 CALL
luerrm(9,
'(LUDECY:) no decay channel defined')
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
72 IF(
rlu(0).GT.0.5) kfs=-kfs
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
92 CALL
luerrm(2,
'(LUDECY:) all decay channels closed by user')
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
108 IF(idl.LT.mdcy(kca,2)+mdcy(kca,3)-1.AND.rbr.GT.0.) goto 140
114 IF(
ntry.GT.1000)
THEN
115 CALL
luerrm(14,
'(LUDECY:) caught in infinite loop')
116 IF(mstu(21).GE.1)
RETURN
122 IF(mmat.GE.11.AND.mmat.NE.46.AND.
p(ip,4).GT.20.*
p(ip,5)) mbst=1
125 160
IF(mbst.EQ.0) pv(1,
j)=
p(ip,
j)
126 IF(mbst.EQ.1) pv(1,4)=
p(ip,5)
134 IF(mdme(idc+1,2).EQ.101) jtmax=10
136 IF(jt.LE.5) kp=kfdp(idc,jt)
137 IF(jt.GE.6) kp=kfdp(idc+1,jt-5)
141 IF(kchg(kcp,3).EQ.0.AND.kpa.NE.81.AND.kpa.NE.82)
THEN
143 ELSEIF(kpa.NE.81.AND.kpa.NE.82)
THEN
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
155 IF(pv(1,5).LT.parj(32)+2.*
ulmass(kfp)) goto 150
156 ELSEIF(kp.EQ.-82)
THEN
158 IF(iabs(kfp).GT.10) kfp=kfp+isign(10000,kfp)
160 IF(kpa.EQ.81.OR.kpa.EQ.82) kcp=
lucomp(kfp)
165 IF(mmat.GE.11.AND.mmat.LE.30.AND.kqp.NE.0)
THEN
170 ELSEIF(mmat.GE.42.AND.mmat.LE.43.AND.
np.EQ.3.AND.mod(nq,2).EQ.1)
177 IF(
k(
i,2).EQ.0) goto 150
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
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
194 IF(mmat.EQ.45.AND.kfpa.EQ.89)
p(
i,5)=parj(32)
200 180
IF(mmat.GE.11.AND.mmat.LE.30)
THEN
202 cnde=parj(61)*
log(
max((pv(1,5)-
ps-psq)/parj(62),1.1))
203 IF(mmat.EQ.12) cnde=cnde+parj(63)
205 IF(
ntry.GT.1000)
THEN
206 CALL
luerrm(14,
'(LUDECY:) caught in infinite loop')
207 IF(mstu(21).GE.1)
RETURN
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
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))
228 IF(
k(
i,2).EQ.0) goto 190
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
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
252 IF(
ps+parj(64).GT.pv(1,5)) goto 190
255 ELSEIF((mmat.EQ.31.OR.mmat.EQ.33.OR.mmat.EQ.44.OR.mmat.EQ.45).
258 pqt=(
p(
n+
np,5)+parj(65))/pv(1,5)
261 240 pv(1,
j)=(1.-pqt)*pv(1,
j)
262 IF(
ps+parj(64).GT.pv(1,5)) goto 150
267 ELSEIF(mmat.EQ.46)
THEN
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
281 IF(
np.GE.2.AND.
ps+parj(64).GT.pv(1,5)) goto 150
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)
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)
303 hmv1=hmeps(hm/hqw)/((hm-1.)**2+hg**2)
305 hmv2=hmeps(hm/hqw)/((hm-1.)**2+hg**2)
307 hsav2=1./((hm-1.)**2+hg**2)
308 IF(hmv2.GT.hmv1.AND.hm-hg.GT.hlw)
THEN
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
316 ELSEIF(hmv2.LE.hmv1)
THEN
317 hm=
max(hlw,hm-min(0.1,1.-hm))
319 hatm=atan((hm-1.)/hg)
321 hwt2=hmv*(min(1.,huw)-hm)
324 hatu=atan((huw-1.)/hg)
330 270 hreg=
rlu(0)*(hwt1+hwt2+hwt3)
331 IF(hreg.LE.hwt1)
THEN
332 hw=1.+hg*tan(hatl+
rlu(0)*(hatm-hatl))
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
338 hw=1.+hg*tan(
rlu(0)*hatu)
339 hacc=hmeps(hw/hqw)/hmp1
341 IF(hacc.LT.
rlu(0)) goto 270
342 p(
n+1,5)=pmas(24,1)*sqrt(hw)
348 IF(mmat.EQ.3.OR.mmat.EQ.46)
THEN
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)
375 pmax=pv(1,5)-
ps+
p(
n+nd,5)
379 pmin=pmin+
p(
n+il+1,5)
380 300 wtmax=wtmax*pawt(pmax,pmin,
p(
n+il,5))
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))
401 DO 340 il2=il1-1,1,-1
402 IF(rsav.LE.rord(il2)) goto 350
403 340 rord(il2+1)=rord(il2)
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
415 pa=pawt(pv(il,5),pv(il+1,5),
p(
n+il,5))
418 ue(1)=sqrt(1.-ue(3)**2)*cos(
phi)
419 ue(2)=sqrt(1.-ue(3)**2)*sin(
phi)
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)
428 400
p(
n+nd,
j)=pv(nd,
j)
431 410 be(
j)=pv(il,
j)/pv(il,4)
434 bep=be(1)*
p(
i,1)+be(2)*
p(
i,2)+be(3)*
p(
i,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)
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
447 ELSEIF(mmat.EQ.2)
THEN
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
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
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)
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
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
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
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
502 450 pv(1,
j)=pv(1,
j)/(1.-pqt)
509 IF((mmat.EQ.31.OR.mmat.EQ.45).AND.nd.EQ.3)
THEN
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
519 IF(
k(
n+2,2).EQ.0) goto 150
526 ELSEIF(mmat.EQ.44)
THEN
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
536 IF(
k(
n+3,2).EQ.0) goto 150
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
549 hh=(sqrt(hg**2+he*hf)-hg)/(2.*hf)
551 pcor=hh*(
p(
n+1,
j)-
p(
n+2,
j))
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)
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)))
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
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
578 IF(
k(
n+1,2).EQ.0) goto 150
590 490
IF(mmat.EQ.42.AND.iabs(
k(
n+1,2)).LT.10)
THEN
615 520 be(
j)=
p(ip,
j)/
p(ip,4)
618 bep=be(1)*
p(
i,1)+be(2)*
p(
i,2)+be(3)*
p(
i,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)
631 IF(mstj(23).GE.1.AND.mmat.EQ.4.AND.
k(nsav+1,2).EQ.21)
THEN
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)
642 ELSEIF(mstj(23).GE.1.AND.mmat.EQ.4)
THEN
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)
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
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)
659 ELSEIF(mstj(23).GE.1.AND.mmat.EQ.33.AND.iabs(
k(nsav+2,2)).EQ.21)
665 kqp=kchg(kcp,2)*isign(1,
k(nsav+1,2))
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)
673 ELSEIF(mstj(23).GE.1.AND.mmat.EQ.33)
THEN
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)
684 IF(
k(ip,1).EQ.5)
k(ip,1)=15
685 IF(
k(ip,1).LE.10)
k(ip,1)=11