10 IMPLICIT DOUBLE PRECISION(a-
h, o-
z)
14 parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
15 &kexcit=4000000,kdimen=5000000)
16 parameter(maxnur=1000)
18 common/pypart/
npart,npartd,ipart(maxnur),ptpart(maxnur)
20 common/
pydat1/mstu(200),paru(200),mstj(200),parj(200)
21 common/
pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
26 dimension pmth(5,140),
ps(5),pma(100),pmsd(100),iep(100),ipa(100),
27 &kfla(100),kfld(100),kfl(100),itry(100),isi(100),isl(100),dp(100),
28 &dpt(5,4),ksh(0:140),kcii(2),niis(2),iiis(2,2),theiis(2,2),
29 &phiiis(2,2),isii(2),isset(2),iscol(0:140),ischg(0:140),
33 IF(mstj(41).LE.0)
THEN
35 ELSEIF(mstj(41).EQ.1.OR.mstj(41).EQ.11)
THEN
36 IF(qmax.LE.parj(82).AND.ip2.GE.-80)
RETURN
38 IF(qmax.LE.min(parj(82),parj(83),parj(90)).AND.ip2.GE.-80)
44 IF(ip1.GT.0.AND.ip1.LE.min(
n,mstu(4)-mstu(32)).AND.ip2.EQ.0)
THEN
47 ELSEIF(min(ip1,ip2).GT.0.AND.
max(ip1,ip2).LE.min(
n,mstu(4)-
52 ELSEIF(ip1.GT.0.AND.ip1.LE.min(
n,mstu(4)-mstu(32)).AND.ip2.LT.0
53 & .AND.ip2.GE.-80)
THEN
58 ELSEIF(ip1.GT.0.AND.ip1.LE.min(
n,mstu(4)-mstu(32)).AND.
66 &
'(PYSHOW:) failed to reconstruct showering system')
67 IF(mstu(21).GE.1)
RETURN
72 IF((mstj(41).EQ.11.OR.mstj(41).EQ.12).AND.npa.GE.2.AND.
79 CALL
pyptfs(2,0.5d0*qmax,0d0,ptgen)
92 pmth(2,21)=sqrt(pmth(1,21)**2+0.25d0*parj(82)**2)
93 pmth(3,21)=2d0*pmth(2,21)
97 pmth(2,22)=sqrt(pmth(1,22)**2+0.25d0*parj(83)**2)
98 pmth(3,22)=2d0*pmth(2,22)
100 pmth(5,22)=pmth(3,22)
102 IF(mstj(41).GE.2) pmqth1=min(parj(82),parj(83))
103 pmqt1e=min(pmqth1,parj(90))
105 IF(mstj(41).GE.2) pmqth2=min(pmth(2,21),pmth(2,22))
106 pmqt2e=min(pmqth2,0.5d0*parj(90))
109 IF(mstj(41).GE.2) ischg(ifl)=1
112 pmth(2,ifl)=sqrt(pmth(1,ifl)**2+0.25d0*pmqth1**2)
113 pmth(3,ifl)=pmth(2,ifl)+pmqth2
114 pmth(4,ifl)=sqrt(pmth(1,ifl)**2+0.25d0*parj(82)**2)+pmth(2,21)
115 pmth(5,ifl)=sqrt(pmth(1,ifl)**2+0.25d0*parj(83)**2)+pmth(2,22)
118 IF(mstj(41).EQ.2.OR.mstj(41).GE.4) ischg(ifl)=1
119 IF(mstj(41).EQ.2.OR.mstj(41).GE.4) ksh(ifl)=1
121 pmth(2,ifl)=sqrt(pmth(1,ifl)**2+0.25d0*parj(90)**2)
122 pmth(3,ifl)=pmth(2,ifl)+0.5d0*parj(90)
123 pmth(4,ifl)=pmth(3,ifl)
124 pmth(5,ifl)=pmth(3,ifl)
126 pt2min=
max(0.5d0*parj(82),1.1d0*parj(81))**2
128 alfm=
log(pt2min/alams)
138 kfla(
i)=iabs(
k(ipa(
i),2))
143 IF(npa.LE.1) iref(
i)=ir
144 IF(npa.GE.2) iref(
i+1)=ir
148 IF(kfla(
i).LE.8)
THEN
150 IF(mstj(41).GE.2) ischg(ir)=1
151 ELSEIF(kfla(
i).EQ.11.OR.kfla(
i).EQ.13.OR.kfla(
i).EQ.15.OR.
152 & kfla(
i).EQ.17)
THEN
153 IF(mstj(41).EQ.2.OR.mstj(41).GE.4) ischg(ir)=1
154 ELSEIF(kfla(
i).EQ.21)
THEN
156 ELSEIF((kfla(
i).GE.ksusy1+1.AND.kfla(
i).LE.ksusy1+8).OR.
157 & (kfla(
i).GE.ksusy2+1.AND.kfla(
i).LE.ksusy2+8))
THEN
159 ELSEIF(kfla(
i).EQ.ksusy1+21)
THEN
163 ELSEIF(
mstp(148).GE.1.AND.(kfla(
i).EQ.9900443.OR.
164 & kfla(
i).EQ.9900553))
THEN
168 IF(iscol(ir).EQ.1.OR.ischg(ir).EQ.1) ksh(ir)=1
170 IF(iscol(ir).EQ.1.AND.ischg(ir).EQ.1)
THEN
171 pmth(2,ir)=sqrt(pmth(1,ir)**2+0.25d0*pmqth1**2)
172 pmth(3,ir)=pmth(2,ir)+pmqth2
173 pmth(4,ir)=sqrt(pmth(1,ir)**2+0.25d0*parj(82)**2)+pmth(2,21)
174 pmth(5,ir)=sqrt(pmth(1,ir)**2+0.25d0*parj(83)**2)+pmth(2,22)
175 ELSEIF(iscol(ir).EQ.1)
THEN
176 pmth(2,ir)=sqrt(pmth(1,ir)**2+0.25d0*parj(82)**2)
177 pmth(3,ir)=pmth(2,ir)+0.5d0*parj(82)
178 pmth(4,ir)=pmth(3,ir)
179 pmth(5,ir)=pmth(3,ir)
180 ELSEIF(ischg(ir).EQ.1)
THEN
181 pmth(2,ir)=sqrt(pmth(1,ir)**2+0.25d0*parj(90)**2)
182 pmth(3,ir)=pmth(2,ir)+0.5d0*parj(90)
183 pmth(4,ir)=pmth(3,ir)
184 pmth(5,ir)=pmth(3,ir)
186 IF(ksh(ir).EQ.1) pma(
i)=pmth(3,ir)
188 IF(ksh(ir).EQ.0.OR.pma(
i).GT.10d0*qmax) irej=irej+1
193 IF(irej.EQ.npa.AND.ip2.GE.-7)
RETURN
195 IF(npa.EQ.1)
ps(5)=
ps(4)
196 IF(
ps(5).LE.pm+pmqt1e)
RETURN
201 ELSEIF(
k(ip1,3).EQ.
k(ip2,3).AND.
k(ip1,3).GT.0)
THEN
202 kfsrce=iabs(
k(
k(ip1,3),2))
204 ipar1=
max(1,
k(ip1,3))
205 ipar2=
max(1,
k(ip2,3))
206 IF(
k(ipar1,3).EQ.
k(ipar2,3).AND.
k(ipar1,3).GT.0)
207 & kfsrce=iabs(
k(
k(ipar1,3),2))
210 IF(kfsrce.GE.1.AND.kfsrce.LE.8) itypes=1
211 IF(kfsrce.GE.ksusy1+1.AND.kfsrce.LE.ksusy1+8) itypes=2
212 IF(kfsrce.GE.ksusy2+1.AND.kfsrce.LE.ksusy2+8) itypes=2
213 IF(kfsrce.GE.21.AND.kfsrce.LE.24) itypes=3
214 IF(kfsrce.GE.32.AND.kfsrce.LE.34) itypes=3
215 IF(kfsrce.EQ.25.OR.(kfsrce.GE.35.AND.kfsrce.LE.37)) itypes=4
216 IF(kfsrce.GE.ksusy1+22.AND.kfsrce.LE.ksusy1+37) itypes=5
217 IF(kfsrce.EQ.ksusy1+21) itypes=6
221 IF(kfla(1).GE.1.AND.kfla(1).LE.8) itype1=1
222 IF(kfla(1).GE.ksusy1+1.AND.kfla(1).LE.ksusy1+8) itype1=2
223 IF(kfla(1).GE.ksusy2+1.AND.kfla(1).LE.ksusy2+8) itype1=2
224 IF(kfla(1).GE.21.AND.kfla(1).LE.24) itype1=3
225 IF(kfla(1).GE.32.AND.kfla(1).LE.34) itype1=3
226 IF(kfla(1).EQ.25.OR.(kfla(1).GE.35.AND.kfla(1).LE.37)) itype1=4
227 IF(kfla(1).GE.ksusy1+22.AND.kfla(1).LE.ksusy1+37) itype1=5
228 IF(kfla(1).EQ.ksusy1+21) itype1=6
230 IF(kfla(2).GE.1.AND.kfla(2).LE.8) itype2=1
231 IF(kfla(2).GE.ksusy1+1.AND.kfla(2).LE.ksusy1+8) itype2=2
232 IF(kfla(2).GE.ksusy2+1.AND.kfla(2).LE.ksusy2+8) itype2=2
233 IF(kfla(2).GE.21.AND.kfla(2).LE.24) itype2=3
234 IF(kfla(2).GE.32.AND.kfla(2).LE.34) itype2=3
235 IF(kfla(2).EQ.25.OR.(kfla(2).GE.35.AND.kfla(2).LE.37)) itype2=4
236 IF(kfla(2).GE.ksusy1+22.AND.kfla(2).LE.ksusy1+37) itype2=5
237 IF(kfla(2).EQ.ksusy1+21) itype2=6
240 itypmn=min(itype1,itype2)
241 itypmx=
max(itype1,itype2)
243 IF(itype1.GT.itype2) iord=2
245 IF(itype1.EQ.6.OR.itype2.EQ.6) iglui=1
250 IF(npa.EQ.2.AND.mstj(47).GE.1.AND.mpspd.EQ.0)
THEN
251 IF(mstj(38).NE.0)
THEN
255 ELSEIF(mstj(47).GE.6)
THEN
262 IF(itypmn.EQ.1.AND.itypmx.EQ.1.AND.(itypes.EQ.0.OR.
265 IF(kfsrce.EQ.21.OR.kfsrce.EQ.22)
THEN
267 ELSEIF(kfsrce.EQ.23.OR.(kfsrce.EQ.0.AND.
268 &
k(ipa(1),2)+
k(ipa(2),2).EQ.0))
THEN
271 IF(kfsrce.EQ.23)
THEN
274 kannfl=iabs(
k(iannfl,2))
275 IF(kannfl.GE.1.AND.kannfl.LE.18)
ei=kchg(kannfl,1)/3d0
279 vi=ai-4d0*
ei*paru(102)
280 ef=kchg(kfla(1),1)/3d0
281 af=
sign(1d0,ef+0.1d0)
282 vf=af-4d0*ef*paru(102)
283 xwc=1d0/(16d0*paru(102)*(1d0-paru(102)))
286 sqwz=
ps(5)*pmas(23,2)
287 sbwz=1d0/((sh-sqmz)**2+sqwz**2)
288 vect=
ei**2*ef**2+2d0*
ei*vi*ef*vf*xwc*sh*(sh-sqmz)*sbwz+
289 & (vi**2+ai**2)*vf**2*xwc**2*sh**2*sbwz
290 axiv=(vi**2+ai**2)*af**2*xwc**2*sh**2*sbwz
292 alpha=vect/(vect+axiv)
293 ELSEIF(kfsrce.EQ.24.OR.kfsrce.EQ.0)
THEN
297 ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.1.AND.itypes.EQ.5)
THEN
299 ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.3.AND.(itypes.EQ.0.OR.
304 ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.1.AND.itypes.EQ.4)
THEN
306 IF(kfsrce.EQ.25.OR.kfsrce.EQ.35.OR.kfsrce.EQ.37)
THEN
308 ELSEIF(kfsrce.EQ.36)
THEN
311 ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.4.AND.(itypes.EQ.0.OR.
316 ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.2.AND.(itypes.EQ.0.OR.
319 ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.3.AND.(itypes.EQ.0.OR.
322 ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.2.AND.itypes.EQ.4)
THEN
324 ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.4.AND.(itypes.EQ.0.OR.
329 ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.2.AND.(itypes.EQ.0.OR.
332 ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.5.AND.(itypes.EQ.0.OR.
335 ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.5.AND.(itypes.EQ.0.OR.
340 ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.2.AND.itypes.EQ.6)
THEN
342 ELSEIF(itypmn.EQ.1.AND.itypmx.EQ.6.AND.(itypes.EQ.0.OR.
345 ELSEIF(itypmn.EQ.2.AND.itypmx.EQ.6.AND.(itypes.EQ.0.OR.
350 ELSEIF(itypmn.EQ.6.AND.itypmx.EQ.6.AND.itypes.EQ.0)
THEN
359 IF(mstj(50).GE.1.AND.mstj(50).LE.3.AND.npa.EQ.2.AND.kfsrce.EQ.0
360 &.AND.mpspd.EQ.0) miis=mstj(50)
361 IF(mstj(50).GE.4.AND.mstj(50).LE.6.AND.npa.EQ.2.AND.mpspd.EQ.0)
367 IF(kca.NE.0) kcii(
i)=kchg(kca,2)*isign(1,
k(ipa(
i),2))
369 IF(kcii(
i).NE.0)
THEN
371 icsi=mod(
k(ipa(
i),3+
j)/mstu(5),mstu(5))
372 IF(icsi.GT.0.AND.icsi.NE.ipa(1).AND.icsi.NE.ipa(2).AND.
373 & (kcii(
i).EQ.(-1)**(
j+1).OR.kcii(
i).EQ.2))
THEN
380 IF(niis(1)+niis(2).EQ.0) miis=0
400 DO 250
i=3+niis(1),2+niis(1)+niis(2)
402 k(
n+
i,
j)=
k(iiis(2,
i-2-niis(1)),
j)
403 p(
n+
i,
j)=
p(iiis(2,
i-2-niis(1)),
j)
407 CALL
pyrobo(
n+1,
n+2+niis(1)+niis(2),0d0,0d0,-
ps(1)/
ps(4),
410 CALL
pyrobo(
n+1,
n+2+niis(1)+niis(2),0d0,-
phi,0d0,0d0,0d0)
412 CALL
pyrobo(
n+1,
n+2+niis(1)+niis(2),-the,0d0,0d0,0d0,0d0)
417 DO 270
i=3+niis(1),2+niis(1)+niis(2)
418 theiis(2,
i-2-niis(1))=paru(1)-
pyangl(
p(
n+
i,3),
419 & sqrt(
p(
n+
i,1)**2+
p(
n+
i,2)**2))
425 IF(npa.GE.3) CALL
pyrobo(ipa(1),ipa(npa),0d0,0d0,-
ps(1)/
ps(4),
430 IF(
n.GT.mstu(4)-mstu(32)-10)
THEN
431 CALL
pyerrm(11,
'(PYSHOW:) no more memory left in PYJETS')
432 IF(mstu(21).GE.1)
RETURN
460 IF(ksh(ir).EQ.0) goto 290
461 IF(
p(im,5).LT.pmth(2,ir)) goto 290
466 IF(
n+nep.GT.mstu(4)-mstu(32)-10)
THEN
467 CALL
pyerrm(11,
'(PYSHOW:) no more memory left in PYJETS')
468 IF(mstu(21).GE.1)
RETURN
475 IF(
k(im-1,3).EQ.igm) iau=im-1
476 IF(
n.GE.im+1.AND.
k(im+1,3).EQ.igm) iau=im+1
490 ELSEIF(kflm.NE.21)
THEN
493 iref(
n+1-ns)=iref(im-ns)
494 iref(
n+2-ns)=iabs(
k(
n+2,2))
495 ELSEIF(
k(im,5).EQ.21)
THEN
503 iref(
n+1-ns)=iabs(
k(
n+1,2))
504 iref(
n+2-ns)=iabs(
k(
n+2,2))
512 kfld(ip)=iabs(
k(
n+ip,2))
513 IF(kchg(
pycomp(kfld(ip)),2).EQ.0)
k(
n+ip,1)=1
517 IF(ksh(iref(
n+ip-ns)).EQ.1) isi(ip)=1
524 IF(npa.GE.3)
p(
n+
i,4)=
p(ipa(
i),4)
525 p(
n+
i,5)=min(qmax,
ps(5))
527 IF(ip2.LE.-8)
p(
n+
i,5)=
max(
p(
n+
i,5),2d0*pmth(3,ir))
528 IF(isi(
i).EQ.0)
p(
n+
i,5)=
p(ipa(
i),5)
531 IF(mstj(43).LE.2) pem=
v(im,2)
532 IF(mstj(43).GE.3) pem=
p(im,4)
533 p(
n+1,5)=min(
p(im,5),
v(im,1)*pem)
534 p(
n+2,5)=min(
p(im,5),(1d0-
v(im,1))*pem)
535 IF(
k(
n+2,2).EQ.22)
p(
n+2,5)=pmth(1,22)
541 IF(
p(
n+
i,5).LE.pmth(3,ir))
p(
n+
i,5)=pmth(1,ir)
550 IF(inum.EQ.0.AND.isl(
i).EQ.1) inum=
i
553 IF(inum.EQ.0.AND.itry(
i).EQ.0.AND.isi(
i).EQ.1)
THEN
555 IF(
p(
n+
i,5).GE.pmth(2,ir)) inum=
i
561 IF(isi(
i).EQ.1.AND.pmsd(
i).GE.pmqt2e)
THEN
564 IF(rpm.GT.rmax.AND.
p(
n+
i,5).GE.pmth(2,ir))
THEN
575 IF(mpspd.EQ.1.AND.igm.EQ.0.AND.itry(inumt).GE.1)
THEN
576 IF(
k(ip1-1+inum,4).GT.0) inum=3-inum
577 ELSEIF(mpspd.EQ.1.AND.im.EQ.ns+2.AND.itry(inumt).GE.1)
THEN
578 IF(kfld(inumt).NE.21.AND.
k(ip1+2,4).GT.0) inum=3-inum
579 IF(kfld(inumt).EQ.21.AND.
k(ip1+3,4).GT.0) inum=3-inum
586 IF(iep(
i).GT.
n+nep) iep(
i)=
n+1
589 kfl(
i)=iabs(
k(iep(
i),2))
591 itry(inum)=itry(inum)+1
592 IF(itry(inum).GT.200)
THEN
593 CALL
pyerrm(14,
'(PYSHOW:) caught in infinite loop')
594 IF(mstu(21).GE.1)
RETURN
598 IF(ksh(ir).EQ.0) goto 450
599 IF(
p(iep(1),5).LT.pmth(2,ir)) goto 450
603 IF(mpspd.EQ.1.AND.igm.EQ.0)
THEN
604 IF(
k(ip1-1+inum,4).GT.0) ipspd=ip1-1+inum
605 ELSEIF(mpspd.EQ.1.AND.im.EQ.ns+2)
THEN
606 IF(kfl(1).NE.21.AND.
k(ip1+2,4).GT.0) ipspd=ip1+2
607 IF(kfl(1).EQ.21.AND.
k(ip1+3,4).GT.0) ipspd=ip1+3
609 IF(inum.EQ.1.OR.inum.EQ.2)
THEN
611 IF(ipspd.NE.0) isset(inum)=1
615 IF(miis.GE.1.AND.iep(1).LE.ns+3)
THEN
618 IF(iabs(kcii(iii)).EQ.1.AND.niis(iii).EQ.1)
THEN
620 ELSEIF(kcii(iii).EQ.2.AND.niis(iii).EQ.1)
THEN
621 IF(
pyr(0).GT.0.5d0) isii(iii)=1
622 ELSEIF(kcii(iii).EQ.2.AND.niis(iii).EQ.2)
THEN
624 IF(
pyr(0).GT.0.5d0) isii(iii)=2
631 ELSEIF(igm.EQ.0.OR.mstj(43).LE.2)
THEN
634 IF(inum.EQ.1) pmed=
v(im,1)*pem
635 IF(inum.EQ.2) pmed=(1d0-
v(im,1))*pem
637 IF(mod(mstj(43),2).EQ.1)
THEN
640 IF(iscol(ir).EQ.0) zce=0.5d0*parj(90)/pmed
642 zc=0.5d0*(1d0-sqrt(
max(0d0,1d0-(2d0*pmth(2,21)/pmed)**2)))
643 IF(zc.LT.1d-6) zc=(pmth(2,21)/pmed)**2
645 IF(iscol(ir).EQ.0) pmtmpe=0.5d0*parj(90)
646 zce=0.5d0*(1d0-sqrt(
max(0d0,1d0-(2d0*pmtmpe/pmed)**2)))
647 IF(zce.LT.1d-6) zce=(pmtmpe/pmed)**2
650 zce=min(zce,0.49991d0)
651 IF(((mstj(41).EQ.1.AND.zc.GT.0.49d0).OR.(mstj(41).GE.2.AND.
652 &min(zc,zce).GT.0.4999d0)).AND.ipspd.EQ.0)
THEN
653 p(iep(1),5)=pmth(1,ir)
654 v(iep(1),5)=
p(iep(1),5)**2
660 IF(mstj(49).EQ.0.AND.kfl(1).EQ.21)
THEN
661 fbr=6d0*
log((1d0-zc)/zc)+mstj(45)*0.5d0
664 ELSEIF(mstj(49).EQ.0.AND.
mstp(149).GE.0.AND.
665 & (kfl(1).EQ.9900443.OR.kfl(1).EQ.9900553))
THEN
666 fbr=6d0*
log((1d0-zc)/zc)
668 ELSEIF(mstj(49).EQ.0)
THEN
669 fbr=(8d0/3d0)*
log((1d0-zc)/zc)
670 IF(iglui.EQ.1.AND.ir.GE.31) fbr=fbr*(9d0/4d0)
673 ELSEIF(mstj(49).EQ.1.AND.kfl(1).EQ.21)
THEN
674 fbr=(parj(87)+mstj(45)*parj(88))*(1d0-2d0*zc)
675 ELSEIF(mstj(49).EQ.1)
THEN
677 IF(igm.EQ.0.AND.m3jc.GE.1) fbr=4d0*fbr
680 ELSEIF(kfl(1).EQ.21)
THEN
681 fbr=6d0*mstj(45)*(0.5d0-zc)
683 fbr=2d0*
log((1d0-zc)/zc)
687 IF(iscol(ir).EQ.0) fbr=0d0
691 IF(mstj(41).GE.2.AND.ischg(ir).EQ.1)
THEN
692 IF(kfl(1).LE.18)
THEN
693 fbre=(kchg(kfl(1),1)/3d0)**2*2d0*
log((1d0-zce)/zce)
695 IF(mstj(41).EQ.10) fbre=parj(84)*fbre
705 IF(ksh(iri).EQ.1) pm=pmth(2,iri)
708 pms=min(pms,(
p(im,5)-pm2)**2)
713 DO 430 iff=4,mstj(45)
714 IF(pms.GT.4d0*pmth(2,iff)**2) b0=(33d0-2d0*iff)/6d0
717 pmsc=
max(0.5d0*parj(82),pms-pmth(1,ir)**2)
721 ELSEIF(fbr.LT.1d-3)
THEN
723 ELSEIF(mstj(44).LE.0)
THEN
724 pmsqcd=pmsc*exp(
max(-50d0,
log(
pyr(0))*paru(2)/(paru(111)*fbr)))
725 ELSEIF(mstj(44).EQ.1)
THEN
726 pmsqcd=4d0*alams*(0.25d0*pmsc/alams)**(
pyr(0)**(b0/fbr))
728 pmsqcd=pmsc*exp(
max(-50d0,alfm*b0*
log(
pyr(0))/fbr))
731 IF(ipspd.EQ.0) pmsqcd=pmsqcd+pmth(1,ir)**2
732 IF(zc.GT.0.49d0.OR.pmsqcd.LE.pmth(4,ir)**2) pmsqcd=pmth(2,ir)**2
737 IF(mstj(41).GE.2.AND.ischg(ir).EQ.1.AND.ipspd.EQ.0)
THEN
739 pmse=
max(0.5d0*parj(83),pms-pmth(1,ir)**2)
740 IF(fbre.LT.1d-3)
THEN
743 pmsqed=pmse*exp(
max(-50d0,
log(
pyr(0))*paru(2)/
747 pmsqed=pmsqed+pmth(1,ir)**2
748 IF(zce.GT.0.4999d0.OR.pmsqed.LE.pmth(5,ir)**2) pmsqed=
750 IF(pmsqed.GT.pmsqcd)
THEN
757 p(iep(1),5)=sqrt(
v(iep(1),5))
758 IF(
p(iep(1),5).LE.pmth(3,ir))
THEN
759 p(iep(1),5)=pmth(1,ir)
760 v(iep(1),5)=
p(iep(1),5)**2
768 pmsgd1=
p(ipsgd1,5)**2
769 pmsgd2=
p(ipsgd2,5)**2
770 alamps=sqrt(
max(1d-10,(pmsqcd-pmsgd1-pmsgd2)**2-
771 & 4d0*pmsgd1*pmsgd2))
772 z=0.5d0*(pmsqcd*(2d0*
p(ipsgd1,4)/
p(ipspd,4)-1d0)+alamps-
773 & pmsgd1+pmsgd2)/alamps
774 z=
max(0.00001d0,min(0.99999d0,
z))
775 IF(kfl(1).NE.21)
THEN
778 k(iep(1),5)=iabs(
k(ipsgd1,2))
782 ELSEIF(mce.EQ.2)
THEN
783 z=1d0-(1d0-zce)*(zce/(1d0-zce))**
pyr(0)
784 IF(1d0+
z**2.LT.2d0*
pyr(0)) goto 410
789 ELSEIF(mstj(49).EQ.0.AND.
790 & (kfl(1).EQ.9900443.OR.kfl(1).EQ.9900553))
THEN
791 z=(1d0-zc)*(zc/(1d0-zc))**
pyr(0)
793 IF(
mstp(149).LE.0.OR.
pyr(0).GT.0.5d0)
z=1d0-
z
794 IF((1d0-
z*(1d0-
z))**2.LT.
pyr(0)) goto 410
799 ELSEIF(mstj(49).NE.1.AND.kfl(1).NE.21)
THEN
800 z=1d0-(1d0-zc)*(zc/(1d0-zc))**
pyr(0)
802 IF(m3jc.EQ.0.AND.1d0+
z**2.LT.2d0*
pyr(0)) goto 410
804 ELSEIF(mstj(49).EQ.0.AND.mstj(45)*0.5d0.LT.
pyr(0)*fbr)
THEN
805 z=(1d0-zc)*(zc/(1d0-zc))**
pyr(0)
806 IF(
pyr(0).GT.0.5d0)
z=1d0-
z
807 IF((1d0-
z*(1d0-
z))**2.LT.
pyr(0)) goto 410
809 ELSEIF(mstj(49).NE.1)
THEN
811 IF(
z**2+(1d0-
z)**2.LT.
pyr(0)) goto 410
812 kflb=1+int(mstj(45)*
pyr(0))
813 pmq=4d0*pmth(2,kflb)**2/
v(iep(1),5)
814 IF(pmq.GE.1d0) goto 410
815 IF(mstj(44).LE.2.OR.mstj(44).EQ.4)
THEN
816 IF(
z.LT.zc.OR.
z.GT.1d0-zc) goto 410
817 pmq0=4d0*pmth(2,21)**2/
v(iep(1),5)
818 IF(mod(mstj(43),2).EQ.0.AND.(1d0+0.5d0*pmq)*sqrt(1d0-pmq)
819 & .LT.
pyr(0)*(1d0+0.5d0*pmq0)*sqrt(1d0-pmq0)) goto 410
821 IF((1d0+0.5d0*pmq)*sqrt(1d0-pmq).LT.
pyr(0)) goto 410
826 ELSEIF(kfl(1).NE.21)
THEN
827 z=1d0-sqrt(zc**2+
pyr(0)*(1d0-2d0*zc))
829 ELSEIF(
pyr(0)*(parj(87)+mstj(45)*parj(88)).LE.parj(87))
THEN
830 z=zc+(1d0-2d0*zc)*
pyr(0)
833 z=zc+(1d0-2d0*zc)*
pyr(0)
834 kflb=1+int(mstj(45)*
pyr(0))
835 pmq=4d0*pmth(2,kflb)**2/
v(iep(1),5)
836 IF(pmq.GE.1d0) goto 410
841 IF(mce.EQ.1.AND.mstj(44).GE.2.AND.ipspd.EQ.0)
THEN
842 IF(kfl(1).EQ.21.AND.
k(iep(1),5).LT.10.AND.
843 & (mstj(44).EQ.3.OR.mstj(44).EQ.5))
THEN
844 IF(alfm/
log(
v(iep(1),5)*0.25d0/alams).LT.
pyr(0)) goto 410
846 pt2app=
z*(1d0-
z)*
v(iep(1),5)
847 IF(mstj(44).GE.4) pt2app=pt2app*
848 & (1d0-pmth(1,ir)**2/
v(iep(1),5))**2
849 IF(pt2app.LT.pt2min) goto 410
850 IF(alfm/
log(pt2app/alams).LT.
pyr(0)) goto 410
855 IF(kfl(1).EQ.21)
THEN
856 irgd1=iabs(
k(iep(1),5))
860 irgd2=iabs(
k(iep(1),5))
864 ELSEIF(nep.GE.3)
THEN
866 ELSEIF(igm.EQ.0.OR.mstj(43).LE.2)
THEN
867 ped=0.5d0*(
v(im,5)+
v(iep(1),5)-pm2**2)/
p(im,5)
869 IF(iep(1).EQ.
n+1) ped=
v(im,1)*pem
870 IF(iep(1).EQ.
n+2) ped=(1d0-
v(im,1))*pem
872 IF(mod(mstj(43),2).EQ.1)
THEN
873 pmqth3=0.5d0*parj(82)
874 IF(irgd2.EQ.22) pmqth3=0.5d0*parj(83)
875 IF(irgd2.EQ.22.AND.iscol(ir).EQ.0) pmqth3=0.5d0*parj(90)
876 pmq1=(pmth(1,irgd1)**2+pmqth3**2)/
v(iep(1),5)
877 pmq2=(pmth(1,irgd2)**2+pmqth3**2)/
v(iep(1),5)
878 zd=sqrt(
max(0d0,(1d0-
v(iep(1),5)/ped**2)*((1d0-pmq1-pmq2)**2-
882 zd=sqrt(
max(0d0,1d0-
v(iep(1),5)/ped**2))
885 IF(kfl(1).EQ.21.AND.
k(iep(1),5).LT.10.AND.
886 &(mstj(44).EQ.3.OR.mstj(44).EQ.5))
THEN
887 ELSEIF(ipspd.NE.0)
THEN
891 IF(
z.LT.zl.OR.
z.GT.zu) goto 410
893 IF(kfl(1).EQ.21)
v(iep(1),3)=
log(zu*(1d0-zl)/
max(1d-20,zl*
895 IF(kfl(1).NE.21)
v(iep(1),3)=
log((1d0-zl)/
max(1d-10,1d0-zu))
898 IF(mstj(40).NE.0.AND.kfl(1).NE.21.AND.ipspd.EQ.0)
THEN
900 eglu=0.5d0*
ps(5)*(1d0-
z)*(1d0+
v(iep(1),5)/
v(ns+1,5))
904 chi=parj(89)**2/(parj(89)**2+eglu**2)
905 IF(mstj(40).EQ.1)
THEN
906 IF(chi.LT.
pyr(0)) goto 410
907 ELSEIF(mstj(40).EQ.2)
THEN
908 IF(1d0-chi.LT.
pyr(0)) goto 410
918 IF(mce.EQ.2.AND.igm.EQ.0)
THEN
919 x1=
z*(1d0+
v(iep(1),5)/
v(ns+1,5))
920 x2=1d0-
v(iep(1),5)/
v(ns+1,5)
924 qf1=kchg(
pycomp(ki1),1)*isign(1,ki1)/3d0
925 qf2=kchg(
pycomp(ki2),1)*isign(1,ki2)/3d0
926 wshow=qf1**2*(1d0-x1)/x3*(1d0+(x1/(2d0-x2))**2)+
927 & qf2**2*(1d0-x2)/x3*(1d0+(x2/(2d0-x1))**2)
928 wme=(qf1*(1d0-x1)/x3-qf2*(1d0-x2)/x3)**2*(x1**2+x2**2)
929 ELSEIF(mce.EQ.2)
THEN
932 ELSEIF(mstj(49).NE.1.AND.
k(iep(1),2).NE.21)
THEN
936 IF(ir.GE.31.AND.igm.EQ.0)
THEN
940 ELSEIF(ir.GE.31)
THEN
943 pedme=pem*(
v(im,1)+(1d0-
v(im,1))*ps1me/
v(im,5))
944 ecmme=pedme+sqrt(
max(0d0,pedme**2-ps1me+pm2me**2))
945 ELSEIF(
k(im,2).EQ.21)
THEN
949 IF(iep(1).GT.iep(2)) zmme=1d0-zmme
950 pmlme=sqrt(
max(0d0,(
v(im,5)-ps1me-pm2me**2)**2-
951 & 4d0*ps1me*pm2me**2))
952 pedme=pem*(0.5d0*(
v(im,5)-pmlme+ps1me-pm2me**2)+pmlme*zmme)/
954 ecmme=pedme+sqrt(
max(0d0,pedme**2-ps1me+pm2me**2))
959 pedme=pem*(
v(im,1)+(1d0-
v(im,1))*ps1me/
v(im,5))
960 ecmme=pedme+sqrt(
max(0d0,pedme**2-ps1me+pm2me**2))
966 x1=(1d0+ps1me/ecmme**2-r2me**2)*(
z+(1d0-
z)*pm1me**2/ps1me)
967 x2=1d0+r2me**2-ps1me/ecmme**2
969 IF(ir.EQ.iord+30)
THEN
976 IF((m3jcc.GE.16.AND.m3jcc.LE.19).OR.
977 & (m3jcc.GE.26.AND.m3jcc.LE.29).OR.
978 & (m3jcc.GE.36.AND.m3jcc.LE.39).OR.
979 & (m3jcc.GE.46.AND.m3jcc.LE.49).OR.
980 & (m3jcc.GE.56.AND.m3jcc.LE.64)) isprad=0
981 IF(isprad.EQ.1) wme=wme*
max(1d-10,1d0+r1me**2-r2me**2-x1)/
982 &
max(1d-10,2d0-x1-x2)
984 wshow=2d0/(
max(1d-10,2d0-x1-x2)*
985 &
max(1d-10,1d0+r2me**2-r1me**2-x2))
986 IF(iglui.EQ.1.AND.ir.GE.31) wshow=(9d0/4d0)*wshow
987 ELSEIF(mstj(49).NE.1)
THEN
991 x1=
z*(1d0+
v(iep(1),5)/
v(ns+1,5))
992 x2=1d0-
v(iep(1),5)/
v(ns+1,5)
994 wshow=4d0*x3*((1d0-x1)/(2d0-x2)**2+(1d0-x2)/(2d0-x1)**2)
996 IF(mstj(102).GE.2) wme=x3**2-2d0*(1d0+x3)*(1d0-x1)*(1d0-x2)*
1000 IF(wme.LT.
pyr(0)*wshow) goto 410
1004 IF(mce.EQ.1.AND.igm.GT.0.AND.mstj(42).GE.2.AND.ipspd.EQ.0)
THEN
1005 pemao=
v(im,1)*
p(im,4)
1006 IF(iep(1).EQ.
n+2) pemao=(1d0-
v(im,1))*
p(im,4)
1007 IF(ir.GE.31.AND.mstj(42).GE.5)
THEN
1009 ELSEIF(kfl(1).EQ.21.AND.
k(iep(1),5).LE.10.AND.(mstj(42).EQ.4
1010 & .OR.mstj(42).EQ.7))
THEN
1012 ELSEIF(kfl(1).EQ.21.AND.
k(iep(1),5).LE.10.AND.(mstj(42).EQ.3
1013 & .OR.mstj(42).EQ.6))
THEN
1015 pmdao=pmth(2,
k(iep(1),5))
1016 the2id=
z*(1d0-
z)*pemao**2/(
v(iep(1),5)-4d0*pmdao**2)
1019 the2id=
z*(1d0-
z)*pemao**2/
v(iep(1),5)
1020 IF(mstj(42).GE.3.AND.mstj(42).NE.5) the2id=the2id*
1021 & (1d0+pmth(1,ir)**2*(1d0-
z)/(
v(iep(1),5)*
z))**2
1025 440
IF(
k(iaom,5).EQ.22)
THEN
1027 IF(
k(iaom,3).LE.ns) maom=0
1028 IF(maom.EQ.1) goto 440
1030 IF(maom.EQ.1.AND.maod.EQ.1)
THEN
1031 the2im=
v(iaom,1)*(1d0-
v(iaom,1))*
p(iaom,4)**2/
v(iaom,5)
1032 IF(the2id.LT.the2im) goto 410
1037 IF(mstj(48).EQ.1.AND.ipspd.EQ.0)
THEN
1038 IF(nep.EQ.1.AND.im.EQ.ns)
THEN
1039 the2id=
z*(1d0-
z)*
ps(4)**2/
v(iep(1),5)
1040 IF(parj(85)**2*the2id.LT.1d0) goto 410
1041 ELSEIF(nep.EQ.2.AND.iep(1).EQ.ns+2)
THEN
1042 the2id=
z*(1d0-
z)*(0.5d0*
p(im,4))**2/
v(iep(1),5)
1043 IF(parj(85)**2*the2id.LT.1d0) goto 410
1044 ELSEIF(nep.EQ.2.AND.iep(1).EQ.ns+3)
THEN
1045 the2id=
z*(1d0-
z)*(0.5d0*
p(im,4))**2/
v(iep(1),5)
1046 IF(parj(86)**2*the2id.LT.1d0) goto 410
1052 IF(miis.GE.2.AND.iep(1).LE.ns+3)
THEN
1053 the2d=
max((1d0-
z)/
z,
z/(1d0-
z))*
v(iep(1),5)/(0.5d0*
p(im,4))**2
1054 IF(iep(1).EQ.ns+2.AND.isii(1).GE.1)
THEN
1055 IF(the2d.GT.theiis(1,isii(1))**2) goto 410
1056 ELSEIF(iep(1).EQ.ns+3.AND.isii(2).GE.1)
THEN
1057 IF(the2d.GT.theiis(2,isii(2))**2) goto 410
1065 IF(nep.EQ.1) goto 490
1066 IF(nep.EQ.2.AND.
p(iep(1),5)+
p(iep(2),5).GE.
p(im,5)) goto 350
1069 IF(itry(
i).EQ.0.AND.ksh(ir).EQ.1)
THEN
1070 IF(
p(
n+
i,5).GE.pmth(2,ir)) goto 350
1078 pmsum=pmsum+
p(
n+
i,5)
1080 IF(pmsum.GE.
ps(5)) goto 350
1081 ELSEIF(igm.EQ.0.OR.mstj(43).LE.2.OR.mod(mstj(43),2).EQ.0)
THEN
1084 IF(ksh(irda).EQ.0) goto 480
1085 IF(
p(
i1,5).LT.pmth(2,irda)) goto 480
1094 IF(igm.EQ.0.OR.mstj(43).LE.2)
THEN
1095 ped=0.5d0*(
v(im,5)+
v(
i1,5)-
v(
i2,5))/
p(im,5)
1097 IF(
i1.EQ.
n+1) zm=
v(im,1)
1098 IF(
i1.EQ.
n+2) zm=1d0-
v(im,1)
1099 pml=sqrt((
v(im,5)-
v(
n+1,5)-
v(
n+2,5))**2-
1100 & 4d0*
v(
n+1,5)*
v(
n+2,5))
1101 ped=pem*(0.5d0*(
v(im,5)-pml+
v(
i1,5)-
v(
i2,5))+pml*zm)/
1104 IF(mod(mstj(43),2).EQ.1)
THEN
1105 pmqth3=0.5d0*parj(82)
1106 IF(irgd2.EQ.22) pmqth3=0.5d0*parj(83)
1107 IF(irgd2.EQ.22.AND.iscol(irda).EQ.0) pmqth3=0.5d0*parj(90)
1108 pmq1=(pmth(1,irgd1)**2+pmqth3**2)/
v(
i1,5)
1109 pmq2=(pmth(1,irgd2)**2+pmqth3**2)/
v(
i1,5)
1110 zd=sqrt(
max(0d0,(1d0-
v(
i1,5)/ped**2)*((1d0-pmq1-pmq2)**2-
1114 zd=sqrt(
max(0d0,1d0-
v(
i1,5)/ped**2))
1117 IF(irda.EQ.21.AND.irgd1.LT.10.AND.
1118 & (mstj(44).EQ.3.OR.mstj(44).EQ.5))
THEN
1122 IF(
i1.EQ.
n+1.AND.(
v(
i1,1).LT.zl.OR.
v(
i1,1).GT.zu).AND.
1123 & isset(1).EQ.0)
THEN
1125 ELSEIF(
i1.EQ.
n+2.AND.(
v(
i1,1).LT.zl.OR.
v(
i1,1).GT.zu).AND.
1126 & isset(2).EQ.0)
THEN
1130 IF(irda.EQ.21)
v(
i1,4)=
log(zu*(1d0-zl)/
max(1d-20,
1132 IF(irda.NE.21)
v(
i1,4)=
log((1d0-zl)/
max(1d-10,1d0-zu))
1134 IF(isl(1).EQ.1.AND.isl(2).EQ.1.AND.islm.NE.0)
THEN
1137 ELSEIF(isl(1).EQ.1.AND.isl(2).EQ.1)
THEN
1138 zdr1=
max(0d0,
v(
n+1,3)/
max(1d-6,
v(
n+1,4))-1d0)
1139 zdr2=
max(0d0,
v(
n+2,3)/
max(1d-6,
v(
n+2,4))-1d0)
1140 IF(zdr2.GT.
pyr(0)*(zdr1+zdr2)) isl(1)=0
1141 IF(isl(1).EQ.1) isl(2)=0
1142 IF(isl(1).EQ.0) islm=1
1143 IF(isl(2).EQ.0) islm=2
1145 IF(isl(1).EQ.1.OR.isl(2).EQ.1) goto 350
1150 IF(mod(mstj(43),2).EQ.1.AND.(
p(
n+1,5).GE.
1151 & pmth(2,ird1).OR.
p(
n+2,5).GE.pmth(2,ird2)))
THEN
1152 pmq1=
v(
n+1,5)/
v(im,5)
1153 pmq2=
v(
n+2,5)/
v(im,5)
1154 zd=sqrt(
max(0d0,(1d0-
v(im,5)/pem**2)*((1d0-pmq1-pmq2)**2-
1159 IF(
v(im,1).LT.zl.OR.
v(im,1).GT.zu) goto 350
1169 p(
n+1,3)=sqrt(
max(0d0,(
p(ipa(1),4)+
p(
n+1,5))*(
p(ipa(1),4)-
1171 p(
n+1,4)=
p(ipa(1),4)
1173 ELSEIF(igm.EQ.0.AND.nep.EQ.2)
THEN
1174 ped1=0.5d0*(
v(im,5)+
v(
n+1,5)-
v(
n+2,5))/
p(im,5)
1177 p(
n+1,3)=sqrt(
max(0d0,(ped1+
p(
n+1,5))*(ped1-
p(
n+1,5))))
1182 p(
n+2,4)=
p(im,5)-ped1
1185 ELSEIF(nep.GE.3)
THEN
1195 pqs=pqs+
p(
n+
i,5)**2/
p(
n+
i,4)
1198 fac=(
ps(5)-pqs)/(pes-pqs)
1208 pqs=pqs+
p(
n+
i,5)**2/
p(
n+
i,4)
1210 IF(loop.LT.10.AND.abs(pes-
ps(5)).GT.1d-12*
ps(5)) goto 520
1217 pzm=sqrt(
max(0d0,(pem+
p(im,5))*(pem-
p(im,5))))
1218 pmls=(
v(im,5)-
v(
n+1,5)-
v(
n+2,5))**2-4d0*
v(
n+1,5)*
v(
n+2,5)
1221 ELSEIF(
k(im,2).EQ.21.AND.iabs(
k(
n+1,2)).LE.10.AND.
1222 & (mstj(44).EQ.3.OR.mstj(44).EQ.5))
THEN
1223 pts=pmls*zm*(1d0-zm)/
v(im,5)
1224 ELSEIF(mod(mstj(43),2).EQ.1)
THEN
1225 pts=(pem**2*(zm*(1d0-zm)*
v(im,5)-(1d0-zm)*
v(
n+1,5)-
1226 & zm*
v(
n+2,5))-0.25d0*pmls)/pzm**2
1228 pts=pmls*(zm*(1d0-zm)*pem**2/
v(im,5)-0.25d0)/pzm**2
1230 IF(pts.LT.0d0.AND.looppt.LT.10)
THEN
1233 ELSEIF(pts.LT.0d0)
THEN
1236 pt=sqrt(
max(0d0,pts))
1245 IF(mstj(49).NE.1.AND.mod(mstj(46),2).EQ.1.AND.
k(im,2).EQ.21
1246 & .AND.iau.NE.0)
THEN
1247 IF(
k(igm,3).NE.0) mazip=1
1249 IF(iau.EQ.im+1) zau=1d0-
v(igm,1)
1250 IF(mazip.EQ.0) zau=0d0
1251 IF(
k(igm,2).NE.21)
THEN
1252 hazip=2d0*zau/(1d0+zau**2)
1254 hazip=(zau/(1d0-zau*(1d0-zau)))**2
1256 IF(
k(
n+1,2).NE.21)
THEN
1257 hazip=hazip*(-2d0*zm*(1d0-zm))/(1d0-2d0*zm*(1d0-zm))
1259 hazip=hazip*(zm*(1d0-zm)/(1d0-zm*(1d0-zm)))**2
1266 IF(mstj(49).NE.2.AND.mstj(46).GE.2.AND.(
k(
n+1,2).EQ.21.OR.
1267 &
k(
n+2,2).EQ.21).AND.iau.NE.0)
THEN
1268 IF(
k(igm,3).NE.0) mazic=
n+1
1269 IF(
k(igm,3).NE.0.AND.
k(
n+1,2).NE.21) mazic=
n+2
1270 IF(
k(igm,3).NE.0.AND.
k(
n+1,2).EQ.21.AND.
k(
n+2,2).EQ.21.AND.
1271 & zm.GT.0.5d0) mazic=
n+2
1272 IF(
k(iau,2).EQ.22) mazic=0
1274 IF(mazic.EQ.
n+2) zs=1d0-zm
1276 IF(iau.EQ.im-1) zgm=1d0-
v(igm,1)
1277 IF(mazic.EQ.0) zgm=1d0
1278 IF(mazic.NE.0) hazic=(
p(im,5)/
p(igm,5))*
1279 & sqrt((1d0-zs)*(1d0-zgm)/(zs*zgm))
1280 hazic=min(0.95d0,hazic)
1285 560
IF(nep.EQ.2.AND.igm.GT.0)
THEN
1286 IF(
k(im,2).EQ.21.AND.iabs(
k(
n+1,2)).LE.10.AND.
1287 & (mstj(44).EQ.3.OR.mstj(44).EQ.5))
THEN
1288 p(
n+1,4)=0.5d0*(pem*(
v(im,5)+
v(
n+1,5)-
v(
n+2,5))+
1289 & pzm*sqrt(
max(0d0,pmls))*(2d0*zm-1d0))/
v(im,5)
1290 ELSEIF(mod(mstj(43),2).EQ.1)
THEN
1291 p(
n+1,4)=pem*
v(im,1)
1293 p(
n+1,4)=pem*(0.5d0*(
v(im,5)-sqrt(pmls)+
v(
n+1,5)-
v(
n+2,5))+
1294 & sqrt(pmls)*zm)/
v(im,5)
1299 IF(mpspd.EQ.1.AND.igm.EQ.ns+1)
THEN
1301 IF(
k(ipspd,4).GT.0)
THEN
1309 ELSEIF(mpspd.EQ.1.AND.igm.EQ.ns+2)
THEN
1311 IF(
k(ipspd,4).GT.0)
THEN
1313 phipsm=
pyangl(
p(ipspd,1),
p(ipspd,2))
1314 thepsm=
pyangl(
p(ipspd,3),sqrt(
p(ipspd,1)**2+
p(ipspd,2)**2))
1315 CALL
pyrobo(ipsgd1,ipsgd1,0d0,-phipsm,0d0,0d0,0d0)
1316 CALL
pyrobo(ipsgd1,ipsgd1,-thepsm,0d0,0d0,0d0,0d0)
1318 CALL
pyrobo(ipsgd1,ipsgd1,thepsm,phipsm,0d0,0d0,0d0)
1325 IF(
k(im,2).EQ.21.AND.iabs(
k(
n+1,2)).LE.10.AND.
1326 & (mstj(44).EQ.3.OR.mstj(44).EQ.5))
THEN
1327 p(
n+1,3)=0.5d0*(pzm*(
v(im,5)+
v(
n+1,5)-
v(
n+2,5))+
1328 & pem*sqrt(
max(0d0,pmls))*(2d0*zm-1d0))/
v(im,5)
1329 ELSEIF(pzm.GT.0d0)
THEN
1330 p(
n+1,3)=0.5d0*(
v(
n+2,5)-
v(
n+1,5)-
v(im,5)+
1331 & 2d0*pem*
p(
n+1,4))/pzm
1337 p(
n+2,3)=pzm-
p(
n+1,3)
1338 p(
n+2,4)=pem-
p(
n+1,4)
1339 IF(mstj(43).LE.2)
THEN
1340 v(
n+1,2)=(pem*
p(
n+1,4)-pzm*
p(
n+1,3))/
p(im,5)
1341 v(
n+2,2)=(pem*
p(
n+2,4)-pzm*
p(
n+2,3))/
p(im,5)
1347 IF(mstj(43).LE.2)
THEN
1348 bex=
p(igm,1)/
p(igm,4)
1349 bey=
p(igm,2)/
p(igm,4)
1350 bez=
p(igm,3)/
p(igm,4)
1351 ga=
p(igm,4)/
p(igm,5)
1352 gabep=ga*(ga*(bex*
p(im,1)+bey*
p(im,2)+bez*
p(im,3))/(1d0+ga)-
1361 ptimb=sqrt((
p(im,1)+gabep*bex)**2+(
p(im,2)+gabep*bey)**2)
1362 the=
pyangl(
p(im,3)+gabep*bez,ptimb)
1363 IF(ptimb.GT.1d-4)
THEN
1369 dp(1)=cos(the)*cos(
phi)*
p(
i,1)-sin(
phi)*
p(
i,2)+
1370 & sin(the)*cos(
phi)*
p(
i,3)
1371 dp(2)=cos(the)*sin(
phi)*
p(
i,1)+cos(
phi)*
p(
i,2)+
1372 & sin(the)*sin(
phi)*
p(
i,3)
1373 dp(3)=-sin(the)*
p(
i,1)+cos(the)*
p(
i,3)
1375 dbp=bex*dp(1)+bey*dp(2)+bez*dp(3)
1376 dgabp=ga*(ga*dbp/(1d0+ga)+dp(4))
1377 p(
i,1)=dp(1)+dgabp*bex
1378 p(
i,2)=dp(2)+dgabp*bey
1379 p(
i,3)=dp(3)+dgabp*bez
1380 p(
i,4)=ga*(dp(4)+dbp)
1385 IF(mazip.NE.0.OR.mazic.NE.0)
THEN
1391 dpma=dpt(1,1)*dpt(2,1)+dpt(1,2)*dpt(2,2)+dpt(1,3)*dpt(2,3)
1392 dpmd=dpt(1,1)*dpt(3,1)+dpt(1,2)*dpt(3,2)+dpt(1,3)*dpt(3,3)
1393 dpmm=dpt(1,1)**2+dpt(1,2)**2+dpt(1,3)**2
1395 dpt(4,
j)=dpt(2,
j)-dpma*dpt(1,
j)/
max(1d-10,dpmm)
1396 dpt(5,
j)=dpt(3,
j)-dpmd*dpt(1,
j)/
max(1d-10,dpmm)
1398 dpt(4,4)=sqrt(dpt(4,1)**2+dpt(4,2)**2+dpt(4,3)**2)
1399 dpt(5,4)=sqrt(dpt(5,1)**2+dpt(5,2)**2+dpt(5,3)**2)
1400 IF(min(dpt(4,4),dpt(5,4)).GT.0.1d0*parj(82))
THEN
1401 cad=(dpt(4,1)*dpt(5,1)+dpt(4,2)*dpt(5,2)+
1402 & dpt(4,3)*dpt(5,3))/(dpt(4,4)*dpt(5,4))
1404 IF(1d0+hazip*(2d0*cad**2-1d0).LT.
pyr(0)*(1d0+abs(hazip)))
1408 IF(mazic.EQ.
n+2) cad=-cad
1409 IF((1d0-hazic)*(1d0-hazic*cad)/(1d0+hazic**2-2d0*hazic*cad)
1410 & .LT.
pyr(0)) goto 560
1416 IF(mod(miis,2).EQ.1.AND.igm.EQ.ns+1.AND.(
k(
n+1,2).EQ.21.OR.
1417 &
k(
n+2,2).EQ.21))
THEN
1419 IF(isii(iii).GE.1)
THEN
1421 IF(
k(
n+1,2).NE.21) iaziid=
n+2
1422 IF(
k(
n+1,2).EQ.21.AND.
k(
n+2,2).EQ.21.AND.
1423 &
p(
n+1,4).GT.
p(
n+2,4)) iaziid=
n+2
1424 theiid=
pyangl(
p(iaziid,3),sqrt(
p(iaziid,1)**2+
p(iaziid,2)**2))
1425 IF(iii.EQ.2) theiid=paru(1)-theiid
1426 phiiid=
pyangl(
p(iaziid,1),
p(iaziid,2))
1427 hazii=min(0.95d0,theiid/theiis(iii,isii(iii)))
1428 cad=cos(phiiid-phiiis(iii,isii(iii)))
1429 phirel=abs(phiiid-phiiis(iii,isii(iii)))
1430 IF(phirel.GT.paru(1)) phirel=paru(2)-phirel
1431 IF((1d0-hazii)*(1d0-hazii*cad)/(1d0+hazii**2-2d0*hazii*cad)
1432 & .LT.
pyr(0)) goto 560
1437 IF(igm.GE.0)
k(im,1)=14
1440 IF(
n.GT.mstu(4)-mstu(32)-10)
THEN
1441 CALL
pyerrm(11,
'(PYSHOW:) no more memory left in PYJETS')
1442 IF(mstu(21).GE.1)
n=ns
1443 IF(mstu(21).GE.1)
RETURN
1448 600
IF(npa.GE.2)
THEN
1452 IF(ip2.GT.0.AND.ip2.LT.ip1)
k(ns+1,3)=ip2
1463 IF(
k(
i,1).LE.10.AND.
k(
i,2).EQ.22)
THEN
1465 ELSEIF(
k(
i,1).LE.10.AND.iabs(
k(
i,2)).GE.11.AND.
1466 & iabs(
k(
i,2)).LE.18)
THEN
1468 ELSEIF(
k(
i,1).LE.10)
THEN
1469 k(
i,4)=mstu(5)*(
k(
i,4)/mstu(5))
1470 k(
i,5)=mstu(5)*(
k(
i,5)/mstu(5))
1471 ELSEIF(
k(mod(
k(
i,4),mstu(5))+1,2).NE.22)
THEN
1472 id1=mod(
k(
i,4),mstu(5))
1473 IF(kq.EQ.1.AND.
k(
i,2).GT.0) id1=mod(
k(
i,4),mstu(5))+1
1474 IF(kq.EQ.2.AND.(
k(id1,2).EQ.21.OR.
k(id1+1,2).EQ.21).AND.
1475 &
pyr(0).GT.0.5d0) id1=mod(
k(
i,4),mstu(5))+1
1476 id2=2*mod(
k(
i,4),mstu(5))+1-id1
1477 k(
i,4)=mstu(5)*(
k(
i,4)/mstu(5))+id1
1478 k(
i,5)=mstu(5)*(
k(
i,5)/mstu(5))+id2
1479 k(id1,4)=
k(id1,4)+mstu(5)*
i
1480 k(id1,5)=
k(id1,5)+mstu(5)*id2
1481 k(id2,4)=
k(id2,4)+mstu(5)*id1
1482 k(id2,5)=
k(id2,5)+mstu(5)*
i
1484 id1=mod(
k(
i,4),mstu(5))
1486 k(
i,4)=mstu(5)*(
k(
i,4)/mstu(5))+id1
1487 k(
i,5)=mstu(5)*(
k(
i,5)/mstu(5))+id1
1488 IF(kq.EQ.1.OR.
k(id1,1).GE.11)
THEN
1489 k(id1,4)=
k(id1,4)+mstu(5)*
i
1490 k(id1,5)=
k(id1,5)+mstu(5)*
i
1502 the=
pyangl(
p(ipa(1),3),sqrt(
p(ipa(1),1)**2+
p(ipa(1),2)**2))
1506 ELSEIF(npa.EQ.2)
THEN
1511 gabep=ga*(ga*(bex*
p(ipa(1),1)+bey*
p(ipa(1),2)+bez*
p(ipa(1),3))
1512 & /(1d0+ga)-
p(ipa(1),4))
1513 the=
pyangl(
p(ipa(1),3)+gabep*bez,sqrt((
p(ipa(1),1)
1514 & +gabep*bex)**2+(
p(ipa(1),2)+gabep*bey)**2))
1515 phi=
pyangl(
p(ipa(1),1)+gabep*bex,
p(ipa(1),2)+gabep*bey)
1533 IF(
n.LE.ns+npa+iim)
THEN
1538 k(ipa(ip),4)=
k(ipa(ip),4)+ns+iim+ip
1539 k(ipa(ip),5)=
k(ipa(ip),5)+ns+iim+ip
1540 k(ns+iim+ip,3)=ipa(ip)
1541 IF(iim.EQ.1.AND.mstu(16).NE.2)
k(ns+iim+ip,3)=ns+1
1542 IF(
k(ns+iim+ip,1).NE.1)
THEN
1543 k(ns+iim+ip,4)=mstu(5)*ipa(ip)+
k(ns+iim+ip,4)
1544 k(ns+iim+ip,5)=mstu(5)*ipa(ip)+
k(ns+iim+ip,5)