12 IMPLICIT DOUBLE PRECISION(a-
h, o-
z)
19 common/
pydat1/mstu(200),paru(200),mstj(200),parj(200)
20 common/
pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
21 common/
pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
22 common/
pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
25 common/
pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
26 common/
pyint3/xsfx(2,-40:40),isig(1000,3),sigh(1000)
28 common/
pyint7/sigt(0:6,0:6,0:5)
29 common/pyintm/kfival(2,3),nmi(2),imi(2,800,2),nvc(2,-6:6),
30 & xassoc(2,-6:6,240),xpsvc(-6:6,-1:240),pvctot(2,-1:1),
31 & xmi(2,240),pt2mi(240),imisep(0:240)
35 dimension nmul(20),sigm(20),kstr(500,2),vintsv(80),
36 &wdtp(0:400),wdte(0:400,0:5),xpq(-25:25),ksav(4,5),psav(4,5)
37 SAVE xt2,xt2fac,xc2,xts,irbin,rbin,nmul,sigm,p83a,p83b,p83c,
38 &cq2i,cq2r,pik,bdiv,b,plowb,phighb,pallb,s4a,s4b,s4c,powip,
39 &rpwip,b2rpdv,b2rpmx,bavg,vnt145,vnt146,vnt147
43 IF(
mstp(122).GE.1)
WRITE(mstu(11),5000)
mstp(82)
56 DO 110 itry=1,
mstp(83)
57 rsca=0.05d0*((21-ixt2)-
pyr(0))
63 IF(
pyr(0).LE.coef(isub,1))
THEN
64 taut=(2d0*(1d0+sqrt(1d0-xt2))/xt2-1d0)**
pyr(0)
65 tau=xt2*(1d0+taut)**2/(4d0*taut)
67 tau=xt2*(1d0+tan(
pyr(0)*atan(sqrt(1d0/xt2-1d0)))**2)
73 IF(ryst.GT.coef(isub,8)) myst=2
74 IF(ryst.GT.coef(isub,8)+coef(isub,9)) myst=3
81 sigm(ixt2)=sigm(ixt2)+sigs
83 sigsum=sigsum+sigm(ixt2)
85 sigsum=sigsum/(20d0*
mstp(83))
88 IF(sigsum.LT.1.1d0*sigt(0,0,5))
THEN
89 IF(
mstp(122).GE.1)
WRITE(mstu(11),5100)
96 IF(
mstp(122).GE.1)
WRITE(mstu(11),5200)
100 yke=sigsum/
max(1d-10,sigt(0,0,5))
101 p83a=(1d0-
parp(83))**2
105 cq2r=2d0/(1d0+
parp(84)**2)
113 130
IF(iit.EQ.0)
THEN
115 ELSEIF(iit.EQ.1)
THEN
118 xk=xi+(yke-yi)*(xf-xi)/(yf-yi)
122 IF(
mstp(82).EQ.2)
THEN
123 sp=0.5d0*paru(1)*(1d0-exp(-xk))
126 IF(
mstp(82).EQ.3)
THEN
128 ELSEIF(
mstp(82).EQ.4)
THEN
129 deltab=min(0.01d0,0.05d0*
parp(84))
133 deltab=
max(0.02d0,0.02d0*(2d0/powip)**(1d0/powip))
143 IF(
mstp(82).EQ.3)
THEN
144 ov=exp(-b**2)/paru(2)
145 ELSEIF(
mstp(82).EQ.4)
THEN
146 ov=(p83a*exp(-min(50d0,b**2))+
147 & p83b*cq2r*exp(-min(50d0,b**2*cq2r))+
148 & p83c*cq2i*exp(-min(50d0,b**2*cq2i)))/paru(2)
150 ov=exp(-b**powip)/paru(2)
151 so=so+paru(2)*b*deltab*ov
153 IF(ibdiv.EQ.1) sohigh=sohigh+paru(2)*b*deltab*ov
154 pacc=1d0-exp(-min(50d0,paru(1)*xk*ov))
155 sp=
sp+paru(2)*b*deltab*pacc
156 sop=sop+paru(2)*b*deltab*ov*pacc
157 bsp=bsp+b*paru(2)*b*deltab*pacc
158 IF(ibdiv.EQ.0.AND.paru(1)*xk*ov.LT.1d0)
THEN
162 IF(b.LT.1d0.OR.b*pacc.GT.1d-6) goto 140
176 IF(abs(yk-yke).GE.1d-5*yke) goto 130
187 pik=(vnt146/vnt147)*yke
190 plowb=paru(1)*bdiv**2
191 IF(
mstp(82).EQ.3)
THEN
192 phighb=pik*0.5*exp(-bdiv**2)
193 ELSEIF(
mstp(82).EQ.4)
THEN
194 s4a=p83a*exp(-bdiv**2)
195 s4b=p83b*exp(-bdiv**2*cq2r)
196 s4c=p83c*exp(-bdiv**2*cq2i)
197 phighb=pik*0.5*(s4a+s4b+s4c)
198 ELSEIF(
parp(83).GE.1.999d0)
THEN
204 b2rpmx=
max(2d0*rpwip,b2rpdv)
209 ELSEIF(mmul.EQ.2)
THEN
213 IF(
mstp(82).LE.0)
THEN
214 ELSEIF(
mstp(82).EQ.1)
THEN
217 IF(
mint(141).NE.0.OR.
mint(142).NE.0) sigrat=sigrat*
219 xt2fac=sigrat*
vint(149)/(1d0-
vint(149))
220 ELSEIF(
mstp(82).EQ.2)
THEN
222 xt2fac=vnt146*
xsec(96,1)/
max(1d-10,sigt(0,0,5))*
225 xc2=4d0*ckin(3)**2/
vint(2)
226 IF(ckin(3).LE.ckin(5).OR.
mint(82).GE.2) xc2=0d0
230 IF(
mstp(82).LE.2)
RETURN
231 142
IF(
pyr(0)*pallb.LT.plowb)
THEN
235 IF(
mstp(82).EQ.3)
THEN
236 ov=exp(-b**2)/paru(2)
237 ELSEIF(
mstp(82).EQ.4)
THEN
238 ov=(p83a*exp(-min(50d0,b**2))+
239 & p83b*cq2r*exp(-min(50d0,b**2*cq2r))+
240 & p83c*cq2i*exp(-min(50d0,b**2*cq2i)))/paru(2)
242 ov=exp(-b**powip)/paru(2)
245 pacc=1d0-exp(-min(50d0,pik*ov))
247 xt2fac=vnt146*
vint(148)*
xsec(96,1)/
max(1d-10,sigt(0,0,5))*
252 IF(
mstp(82).EQ.3)
THEN
253 b=sqrt(bdiv**2-
log(
pyr(0)))
254 ov=exp(-b**2)/paru(2)
255 ELSEIF(
mstp(82).EQ.4)
THEN
256 s4rndm=
pyr(0)*(s4a+s4b+s4c)
257 IF(s4rndm.LT.s4a)
THEN
258 b=sqrt(bdiv**2-
log(
pyr(0)))
259 ELSEIF(s4rndm.LT.s4a+s4b)
THEN
260 b=sqrt(bdiv**2-
log(
pyr(0))/cq2r)
262 b=sqrt(bdiv**2-
log(
pyr(0))/cq2i)
264 ov=(p83a*exp(-min(50d0,b**2))+
265 & p83b*cq2r*exp(-min(50d0,b**2*cq2r))+
266 & p83c*cq2i*exp(-min(50d0,b**2*cq2i)))/paru(2)
267 ELSEIF(
parp(83).GE.1.999d0)
THEN
268 144 b2rpw=b2rpdv-
log(
pyr(0))
269 accip=(b2rpw/b2rpdv)**rpwip
270 IF(accip.LT.
pyr(0)) goto 144
271 ov=exp(-b2rpw)/paru(2)
274 146 b2rpw=b2rpdv-2d0*
log(
pyr(0))
275 accip=(b2rpw/b2rpmx)**rpwip*exp(-0.5d0*(b2rpw-b2rpmx))
276 IF(accip.LT.
pyr(0)) goto 146
277 ov=exp(-b2rpw)/paru(2)
281 pacc=(1d0-exp(-min(50d0,pik*ov)))/(pik*ov)
283 IF(pacc.LT.
pyr(0)) goto 142
286 ELSEIF(mmul.EQ.3)
THEN
294 IF(
mstp(82).LE.0)
THEN
296 ELSEIF(
mstp(82).EQ.1)
THEN
297 xt2=xt2fac*xt2/(xt2fac-xt2*
log(
pyr(0)))
299 ELSEIF(
mstp(82).EQ.2.OR.
mint(39).EQ.1)
THEN
300 IF(xt2.LT.1d0.AND.exp(-xt2fac*xt2/(
vint(149)*(xt2+
301 &
vint(149)))).GT.
pyr(0)) xt2=1d0
303 xt2=(1d0+
vint(149))*xt2fac/(xt2fac-(1d0+
vint(149))*
log(1d0-
304 &
pyr(0)*(1d0-exp(-xt2fac/(
vint(149)*(1d0+
vint(149)))))))-
307 xt2=-xt2fac/
log(exp(-xt2fac/(xt2+
vint(149)))+
pyr(0)*
308 & (exp(-xt2fac/
vint(149))-exp(-xt2fac/(xt2+
vint(149)))))-
321 IF(
mstp(82).LE.1.AND.xt2.LT.
vint(149))
THEN
334 IF(
pyr(0).LE.coef(isub,1))
THEN
335 taut=(2d0*(1d0+sqrt(1d0-xt2))/xt2-1d0)**
pyr(0)
336 tau=xt2*(1d0+taut)**2/(4d0*taut)
338 tau=xt2*(1d0+tan(
pyr(0)*atan(sqrt(1d0/xt2-1d0)))**2)
344 IF(ryst.GT.coef(isub,8)) myst=2
345 IF(ryst.GT.coef(isub,8)+coef(isub,9)) myst=3
352 ELSEIF(mmul.EQ.4)
THEN
358 IF(iset(isub).EQ.1) xts=
vint(21)
361 IF(iset(isub).GE.3.AND.iset(isub).LE.5) xts=
vint(26)
362 rbin=
max(0.000001d0,min(0.999999d0,xts*(1d0+
vint(149))/
364 irbin=int(1d0+20d0*rbin)
365 IF(isub.EQ.96.AND.
mstp(171).EQ.0)
THEN
366 nmul(irbin)=nmul(irbin)+1
367 sigm(irbin)=sigm(irbin)+
vint(153)
371 ELSEIF(mmul.EQ.5)
THEN
376 150
IF(
mint(39).GT.0)
THEN
377 ELSEIF(
mstp(82).EQ.3)
THEN
380 vint(148)=expb2/(paru(2)*vnt147)
381 vint(139)=sqrt(b2)/bavg
382 ELSEIF(
mstp(82).EQ.4)
THEN
384 IF(rtype.LT.p83a)
THEN
386 ELSEIF(rtype.LT.p83a+p83b)
THEN
391 vint(148)=(p83a*exp(-min(50d0,b2))+
392 & p83b*cq2r*exp(-min(50d0,b2*cq2r))+
393 & p83c*cq2i*exp(-min(50d0,b2*cq2i)))/(paru(2)*vnt147)
394 vint(139)=sqrt(b2)/bavg
395 ELSEIF(
parp(83).GE.1.999d0)
THEN
398 prob1=powip/(2d0*exp(-1d0)+powip)
399 160
IF(
pyr(0).LT.prob1)
THEN
400 b2rpw=
pyr(0)**(0.5d0*powip)
406 IF(accip.LT.
pyr(0)) goto 160
407 vint(148)=exp(-b2rpw)/(paru(2)*vnt147)
408 vint(139)=b2rpw**(1d0/powip)/bavg
412 prob1=rpwip/(rpwip+2d0**rpwip*exp(-rpwip))
413 170
IF(
pyr(0).LT.prob1)
THEN
414 b2rpw=2d0*rpwip*
pyr(0)
415 accip=(b2rpw/rpwip)**rpwip*exp(rpwip-b2rpw)
417 b2rpw=2d0*(rpwip-
log(
pyr(0)))
418 accip=(0.5d0*b2rpw/rpwip)**rpwip*exp(rpwip-0.5d0*b2rpw)
420 IF(accip.lt .
pyr(0)) goto 170
421 vint(148)=exp(-b2rpw)/(paru(2)*vnt147)
422 vint(139)=b2rpw**(1d0/powip)/bavg
429 IF(
mint(39).NE.1)
THEN
430 rncor=(irbin-20d0*rbin)*nmul(irbin)
431 sigcor=(irbin-20d0*rbin)*sigm(irbin)
432 DO 180 ibin=irbin+1,20
433 rncor=rncor+nmul(ibin)
434 sigcor=sigcor+sigm(ibin)
436 sigabv=(sigcor/rncor)*
vint(149)*(1d0-xts)/(xts+
vint(149))
438 vint(150)=exp(-min(50d0,vnt146*
vint(148)*
439 & sigabv/
max(1d-10,sigt(0,0,5))))
441 IF(
mstp(86).EQ.3.OR.(
mstp(86).EQ.2.AND.isub.NE.11.AND.
442 & isub.NE.12.AND.isub.NE.13.AND.isub.NE.28.AND.isub.NE.53
443 & .AND.isub.NE.68.AND.isub.NE.95.AND.isub.NE.96))
THEN
444 IF(
vint(150).LT.
pyr(0)) goto 150
449 ELSEIF(mmul.EQ.6)
THEN
487 kfsbm=isign(1,kfbeam)
494 IF(kfabm.GT.1000)
THEN
495 kfival(js,1)=kfsbm*mod(kfabm/1000,10)
496 kfival(js,2)=kfsbm*mod(kfabm/100,10)
497 kfival(js,3)=kfsbm*mod(kfabm/10,10)
499 ELSEIF(kfabm.EQ.211)
THEN
502 ELSEIF(kfabm.EQ.321)
THEN
503 kfival(js,1)=-kfsbm*3
516 IF(
k(
i,3).EQ.
mint(83)+2+js)
THEN
524 IF (iabs(ifl).GT.6) goto 230
529 IF(
mstp(61).GE.1)
THEN
541 rvcs=
pyr(0)*(sea+val)
543 220
IF (rvcs.LE.val.AND.ivnow.GE.1)
THEN
546 IF(kfival(js,1).EQ.ifl) ivnow=ivnow+1
547 IF(kfival(js,2).EQ.ifl) ivnow=ivnow+1
548 IF(kfival(js,3).EQ.ifl) ivnow=ivnow+1
549 IF(kfival(js,1).EQ.0)
THEN
550 IF(kfbeam.EQ.111.AND.iabs(ifl).LE.2) ivnow=1
551 IF(kfbeam.EQ.22.AND.iabs(ifl).LE.5) ivnow=1
552 IF((kfbeam.EQ.130.OR.kfbeam.EQ.310).AND.
553 & (iabs(ifl).EQ.1.OR.iabs(ifl).EQ.3)) ivnow=1
555 IF(ivnow.EQ.0) goto 220
559 IF(kfival(js,1).EQ.0)
THEN
560 IF(kfbeam.EQ.111.OR.kfbeam.EQ.22)
THEN
563 ELSEIF(kfbeam.EQ.130.OR.kfbeam.EQ.310)
THEN
565 IF(iabs(ifl).EQ.1) kfival(js,2)=isign(3,-ifl)
566 IF(iabs(ifl).NE.1) kfival(js,2)=isign(1,-ifl)
572 nvc(js,-ifl)=nvc(js,-ifl)+1
573 xassoc(js,-ifl,nvc(js,-ifl))=
x
575 imi(js,1,2)=-nvc(js,-ifl)
584 IF(
mstp(86).EQ.3.OR.(
mstp(86).EQ.2.AND.isubsv.NE.11.AND.
585 & isubsv.NE.12.AND.isubsv.NE.13.AND.isubsv.NE.28.AND.
586 & isubsv.NE.53.AND.isubsv.NE.68.AND.isubsv.NE.95.AND.
587 & isubsv.NE.96))
THEN
591 IF(iset(isubsv).EQ.1) xt2=
vint(21)
592 IF(iset(isubsv).EQ.2)
594 IF(iset(isubsv).GE.3.AND.iset(isubsv).LE.5) xt2=
vint(26)
596 IF(
mstp(82).LE.1)
THEN
598 IF(
mint(141).NE.0.OR.
mint(142).NE.0) sigrat=sigrat*
600 xt2fac=sigrat*
vint(149)/(1d0-
vint(149))
602 xt2fac=vnt146*
vint(148)*
xsec(isub,1)/
609 240
IF((
mint(35).EQ.2.AND.
mstp(81).EQ.10).OR.isubsv.EQ.95)
THEN
612 ELSEIF(
mstp(82).LE.1)
THEN
613 xt2=xt2fac*xt2/(xt2fac-xt2*
log(
pyr(0)))
614 IF(xt2.LT.
vint(149)) goto 440
616 IF(xt2.LE.0.01001d0*
vint(149)) goto 440
617 xt2=xt2fac*(xt2+
vint(149))/(xt2fac-(xt2+
vint(149))*
619 IF(xt2.LE.0d0) goto 440
625 IF(
pyr(0).LE.coef(isub,1))
THEN
626 taut=(2d0*(1d0+sqrt(1d0-xt2))/xt2-1d0)**
pyr(0)
627 tau=xt2*(1d0+taut)**2/(4d0*taut)
629 tau=xt2*(1d0+tan(
pyr(0)*atan(sqrt(1d0/xt2-1d0)))**2)
633 IF(
tau*
vint(2).LT.1d0) goto 240
637 IF(ryst.GT.coef(isub,8)) myst=2
638 IF(ryst.GT.coef(isub,8)+coef(isub,9)) myst=3
645 IF(
vint(143)-x1m.LT.0.01d0.OR.
vint(144)-x2m.LT.0.01d0) goto 240
648 IF(
mint(141).NE.0.OR.
mint(142).NE.0) sigs=sigs*
vint(320)
649 IF(sigs.LT.
xsec(isub,1)*
pyr(0)) goto 240
650 IF(
mint(141).NE.0.OR.
mint(142).NE.0) sigs=sigs/
vint(320)
660 pt=0.5d0*
vint(1)*sqrt(xt2)
668 rsigs=rsigs-sigh(ichn)
669 IF(rsigs.LE.0d0) goto 280
674 iconmi=mod(iconmi,10)
677 IF(isubmi.EQ.12.OR.isubmi.EQ.53)
THEN
679 CALL
pywidt(21,sh,wdtp,wdte)
680 290 rkfl=(wdte(0,1)+wdte(0,2)+wdte(0,4))*
pyr(0)
681 DO 300
i=1,mdcy(21,3)
682 kflf=kfdp(
i+mdcy(21,2)-1,1)
683 rkfl=rkfl-(wdte(
i,1)+wdte(
i,2)+wdte(
i,4))
684 IF(rkfl.LE.0d0) goto 310
686 310
IF(isubmi.EQ.53.AND.iconmi.LE.2)
THEN
687 IF(kflf.GE.4) goto 290
688 ELSEIF(isubmi.EQ.53.AND.iconmi.LE.4)
THEN
691 ELSEIF(isubmi.EQ.53)
THEN
704 IF(isubmi.EQ.11)
THEN
707 IF(kfl1*kfl2.LT.0) kcc=kcc+2
709 ELSEIF(isubmi.EQ.12)
THEN
711 kfl3=isign(kflf,kfl1)
715 ELSEIF(isubmi.EQ.13)
THEN
721 ELSEIF(isubmi.EQ.28)
THEN
725 IF(kfl1.EQ.21) kcc=kcc+2
726 IF(kfl1.NE.21) kcs=isign(1,kfl1)
727 IF(kfl2.NE.21) kcs=isign(1,kfl2)
729 ELSEIF(isubmi.EQ.53)
THEN
731 kcs=(-1)**int(1.5d0+
pyr(0))
736 ELSEIF(isubmi.EQ.68)
THEN
739 kcs=(-1)**int(1.5d0+
pyr(0))
768 IF(icol(kcc,1,jc).NE.0)
k(
n+1,
j+3)=
n+icol(kcc,1,jc)
769 IF(icol(kcc,2,jc).NE.0)
k(
n+2,
j+3)=
n+icol(kcc,2,jc)
770 IF(icol(kcc,3,jc).NE.0)
k(
n+3,
j+3)=mstu(5)*(
n+icol(kcc,3,jc))
771 IF(icol(kcc,4,jc).NE.0)
k(
n+4,
j+3)=mstu(5)*(
n+icol(kcc,4,jc))
782 IF(
p(
n+3,5)+
p(
n+4,5).GE.shr) goto 240
783 p(
n+3,4)=0.5d0*(shr+(
p(
n+3,5)**2-
p(
n+4,5)**2)/shr)
784 p(
n+3,3)=sqrt(
max(0d0,
p(
n+3,4)**2-
p(
n+3,5)**2))
785 p(
n+4,4)=shr-
p(
n+3,4)
804 IF(
mstp(84).GE.1.AND.
mstp(61).GE.1)
THEN
818 IF(
mint(51).GE.1)
THEN
834 370 imi(1,
mint(31),1)=ipu1
835 imi(2,
mint(31),1)=ipu2
847 kfsbm=isign(1,
mint(10+js))
848 ifl=
k(imi(js,
mint(31),1),2)
850 IF (iabs(ifl).GT.6) goto 430
856 IF(
mstp(84).GE.1.AND.
mstp(61).GE.1)
THEN
867 DO 380 ivc=1,nvc(js,ifl)
868 cmp=cmp+xpsvc(ifl,ivc)
872 rvcs=
pyr(0)*(sea+val+cmp)
874 390
IF (rvcs.LE.val.AND.ivnow.GE.1)
THEN
877 IF(kfival(js,1).EQ.ifl) ivnow=ivnow+1
878 IF(kfival(js,2).EQ.ifl) ivnow=ivnow+1
879 IF(kfival(js,3).EQ.ifl) ivnow=ivnow+1
880 IF(kfival(js,1).EQ.0)
THEN
881 IF(kfbeam.EQ.111.AND.iabs(ifl).LE.2) ivnow=1
882 IF(kfbeam.EQ.22.AND.iabs(ifl).LE.5) ivnow=1
883 IF((kfbeam.EQ.130.OR.kfbeam.EQ.310).AND.
884 & (iabs(ifl).EQ.1.OR.iabs(ifl).EQ.3)) ivnow=1
887 IF (
k(imi(js,
i1,1),2).EQ.ifl.AND.imi(js,
i1,2).EQ.0)
891 IF(ivnow.EQ.0) goto 390
895 IF(kfival(js,1).EQ.0)
THEN
896 IF(kfbeam.EQ.111.OR.kfbeam.EQ.22)
THEN
899 ELSEIF(kfbeam.EQ.130.OR.kfbeam.EQ.310)
THEN
901 IF(iabs(ifl).EQ.1) kfival(js,2)=isign(3,-ifl)
902 IF(iabs(ifl).NE.1) kfival(js,2)=isign(1,-ifl)
906 ELSEIF (rvcs.LE.val+sea.OR.nvc(js,ifl).EQ.0)
THEN
908 nvc(js,-ifl)=nvc(js,-ifl)+1
909 xassoc(js,-ifl,nvc(js,-ifl))=
x
911 imi(js,
mint(31),2)=-nvc(js,-ifl)
917 cmpsum=cmpsum+xpsvc(ifl,isel)
918 IF (rvcs.GT.cmpsum.AND.isel.LT.nvc(js,ifl)) goto 410
922 IF (
k(imi(js,
i1,1),2).NE.-ifl) goto 420
923 IF (-imi(js,
i1,2).EQ.isel)
THEN
924 imi(js,
mint(31),2)=imi(js,
i1,1)
925 imi(js,
i1,2)=imi(js,
mint(31),1)
930 x=xassoc(js,ifl,isel)
932 xassoc(js,ifl,isel)=0d0
942 IF(
n.GT.mstu(4)-mstu(32)-10)
THEN
943 CALL
pyerrm(11,
'(PYMIGN:) no more memory left in PYJETS')
955 IF(
mint(31).LT.240) goto 240
975 5000
FORMAT(/1
x,
'****** PYMIGN: initialization of multiple inter',
976 &
'actions for MSTP(82) =',
i2,
' ******')
977 5100
FORMAT(8
x,
'pT0 =',f5.2,
' GeV gives sigma(parton-parton) =',1
p,
978 &d9.2,
' mb: rejected')
979 5200
FORMAT(8
x,
'pT0 =',f5.2,
' GeV gives sigma(parton-parton) =',1
p,
980 &d9.2,
' mb: accepted')