13 IMPLICIT DOUBLE PRECISION(a-
h, o-
z)
18 common/
pydat1/mstu(200),paru(200),mstj(200),parj(200)
19 common/
pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
22 common/
pyint8/xpvmd(-6:6),xpanl(-6:6),xpanh(-6:6),xpbeh(-6:6),
24 common/
pyint9/vxpvmd(-6:6),vxpanl(-6:6),vxpanh(-6:6),vxpdgm(-6:6)
25 common/pyintm/kfival(2,3),nmi(2),imi(2,800,2),nvc(2,-6:6),
26 & xassoc(2,-6:6,240),xpsvc(-6:6,-1:240),pvctot(2,-1:1),
27 & xmi(2,240),pt2mi(240),imisep(0:240)
31 dimension xpq(-25:25),xpel(-25:25),xpga(-6:6),vxpga(-6:6),
32 &xppi(-6:6),xppr(-6:6),xpval(-6:6),ppar(6,2)
38 DOUBLE PRECISION xx,qq,upv,dnv,usea,dsea,
str,chm,bot,
top,glu,
41 DATA value/20*0d0/,parm/20*
' '/
45 common/pynucl/inumod,chanum,order
47 DOUBLE PRECISION inumod,chanum
54 DATA alamga/0.2d0/, pmcga/1.3d0/, pmbga/4.6d0/
57 DATA (ppar(1,ipar),ipar=1,2) /0.385d0,1.60d0/
58 DATA (ppar(2,ipar),ipar=1,2) /0.480d0,1.56d0/
59 pavg(ifl,q2)=ppar(ifl,1)/(1d0+ppar(ifl,2)*
72 IF(
x.LE.0d0.OR.
x.GE.1d0)
THEN
73 WRITE(mstu(11),5000)
x
77 IF(kfa.NE.11.AND.kfa.NE.13.AND.kfa.NE.15.AND.kfa.NE.22.AND.
78 &kfa.NE.211.AND.kfa.NE.2112.AND.kfa.NE.2212.AND.kfa.NE.3122.AND.
79 &kfa.NE.3112.AND.kfa.NE.3212.AND.kfa.NE.3222.AND.kfa.NE.3312.AND.
80 &kfa.NE.3322.AND.kfa.NE.3334.AND.kfa.NE.111.AND.kfa.NE.321.AND.
81 &kfa.NE.310.AND.kfa.NE.130)
THEN
82 WRITE(mstu(11),5100) kf
87 IF(kfa.EQ.11.OR.kfa.EQ.13.OR.kfa.EQ.15)
THEN
94 ELSEIF(kfa.EQ.22.AND.
mint(109).LE.1)
THEN
95 IF(
mstp(56).EQ.1.AND.
mstp(55).EQ.1)
THEN
100 xpvu=4d0*(xpq(2)-xpq(1))/3d0
103 xpval(3)=min(xpq(3),xpvu/4d0)
104 xpval(4)=min(xpq(4),xpvu)
105 xpval(5)=min(xpq(5),xpvu/4d0)
111 ELSEIF(
mstp(56).EQ.1.AND.
mstp(55).GE.5.AND.
mstp(55).LE.8)
THEN
114 IF(
mstp(55).GE.7) p2mx=4.0d0
115 IF(
mstp(57).EQ.0) q2mx=p2mx
117 IF(
vint(120).LT.0d0) p2=
vint(120)**2
121 xpval(kfl)=vxpdgm(kfl)
124 ELSEIF(
mstp(56).EQ.1.AND.
mstp(55).GE.9.AND.
mstp(55).LE.12)
THEN
127 IF(
mstp(55).GE.11) p2mx=4.0d0
128 IF(
mstp(57).EQ.0) q2mx=p2mx
130 IF(
vint(120).LT.0d0) p2=
vint(120)**2
133 xpq(kfl)=xpvmd(kfl)+xpanl(kfl)+xpbeh(kfl)+xpdir(kfl)
134 xpval(kfl)=vxpvmd(kfl)+vxpanl(kfl)+xpbeh(kfl)+xpdir(kfl)
137 ELSEIF(
mstp(56).EQ.2)
THEN
145 IF(
mint(93).NE.3000000+
mstp(55))
THEN
150 qq2=
max(0d0,q2min,q2)
151 IF(
mstp(57).EQ.0) qq2=q2min
153 IF(
vint(120).LT.0d0) p2=
vint(120)**2
155 IF(
mstp(55).EQ.5004)
THEN
156 IF(5d0*p2.LT.qq2.AND.
157 & qq2.GT.0.6d0.AND.qq2.LT.5d4.AND.
158 & p2.GE.0d0.AND.p2.LT.10d0.AND.
159 & xx.GT.1d-4.AND.xx.LT.1d0)
THEN
160 CALL
structp(xx,qq2,p2,ip2,upv,dnv,usea,dsea,
str,chm,
175 CALL
structp(xx,qq2,p2,ip2,upv,dnv,usea,dsea,
str,chm,
203 xpvu=4d0*(xpq(2)-xpq(1))/3d0
206 xpval(3)=min(xpq(3),xpvu/4d0)
207 xpval(4)=min(xpq(4),xpvu)
208 xpval(5)=min(xpq(5),xpvu/4d0)
215 WRITE(mstu(11),5200) kf,
mstp(56),
mstp(55)
219 ELSEIF(kfa.EQ.211.OR.kfa.EQ.111.OR.kfa.EQ.321.OR.kfa.EQ.130.OR.
220 &kfa.EQ.310.OR.(kfa.EQ.22.AND.
mint(109).EQ.2))
THEN
221 IF(kfa.EQ.22.AND.
mstp(56).EQ.1.AND.
mstp(55).GE.5.AND.
222 &
mstp(55).LE.12)
THEN
223 iset=1+mod(
mstp(55)-1,4)
226 IF(iset.GE.3) p2mx=4.0d0
227 IF(
mstp(57).EQ.0) q2mx=p2mx
229 IF(
vint(120).LT.0d0) p2=
vint(120)**2
233 xpval(kfl)=vxpvmd(kfl)
236 ELSEIF(
mstp(54).EQ.1.AND.
mstp(53).GE.1.AND.
mstp(53).LE.3)
THEN
241 xpval(2)=xpq(2)-xpq(-2)
242 xpval(-1)=xpq(-1)-xpq(1)
243 ELSEIF(
mstp(54).EQ.2)
THEN
251 IF(
mint(93).NE.2000000+
mstp(53))
THEN
256 qq=sqrt(
max(0d0,q2min,q2))
257 IF(
mstp(57).EQ.0) qq=sqrt(q2min)
258 CALL
structm(xx,qq,upv,dnv,usea,dsea,
str,chm,bot,
top,glu)
276 WRITE(mstu(11),5200) kf,
mstp(54),
mstp(53)
280 ELSEIF(kfa.EQ.22.AND.
mint(109).EQ.3)
THEN
283 IF(
mstp(56).EQ.1.AND.
mstp(55).LE.8)
THEN
284 IF(
mstp(55).EQ.5.OR.
mstp(55).EQ.6) p2mx=0.36d0
285 IF(
mstp(55).EQ.7.OR.
mstp(55).EQ.8) p2mx=4.0d0
286 IF(
mstp(57).EQ.0) q2mx=p2mx
288 IF(
vint(120).LT.0d0) p2=
vint(120)**2
291 xpq(kfl)=xpanl(kfl)+xpanh(kfl)
292 xpval(kfl)=vxpanl(kfl)+vxpanh(kfl)
295 ELSEIF(
mstp(56).EQ.1)
THEN
296 IF(
mstp(55).EQ.9.OR.
mstp(55).EQ.10) p2mx=0.36d0
297 IF(
mstp(55).EQ.11.OR.
mstp(55).EQ.12) p2mx=4.0d0
298 IF(
mstp(57).EQ.0) q2mx=p2mx
300 IF(
vint(120).LT.0d0) p2=
vint(120)**2
303 xpq(kfl)=
max(0d0,xpanl(kfl)+xpbeh(kfl)+xpdir(kfl))
304 xpval(kfl)=
max(0d0,vxpanl(kfl)+xpbeh(kfl)+xpdir(kfl))
307 ELSEIF(
mstp(56).EQ.2)
THEN
308 IF(
mstp(57).EQ.0) q2mx=p2mx
309 CALL
pygano(0,
x,q2mx,p2mx,alamga,xpga,vxpga)
312 xpval(kfl)=vxpga(kfl)
315 ELSEIF(
mstp(55).GE.1.AND.
mstp(55).LE.5)
THEN
316 IF(
mstp(57).EQ.0) q2mx=p2mx
320 xpval(kfl)=vxpga(kfl)
329 IF(rkf.GT.10d0) kfr=5
330 IF(kfr.EQ.4.AND.q2.LT.pmcga**2) goto 220
331 IF(kfr.EQ.5.AND.q2.LT.pmbga**2) goto 220
332 IF(
mstp(57).EQ.0) q2mx=p2mx
336 xpval(kfl)=vxpga(kfl)
343 IF(
mstp(52).EQ.1.AND.
mstp(51).GE.1.AND.
mstp(51).LE.20)
THEN
348 xpval(1)=xpq(1)-xpq(-1)
349 xpval(2)=xpq(2)-xpq(-2)
350 ELSEIF(
mstp(52).EQ.2)
THEN
358 IF(
mint(93).NE.1000000+
mstp(51))
THEN
364 qq=sqrt(
max(0d0,q2min,q2))
365 IF(
mstp(57).EQ.0) qq=sqrt(q2min)
367 IF(inumod.GT.1d0)
THEN
368 IF(order/100.EQ.1)
THEN
371 write(chpset,
'(I2)') ipset
372 nprm=
'EPS09LO,'//chpset
373 CALL setlhaparm(nprm)
374 ELSE IF(order/100.EQ.2)
THEN
377 write(chpset,
'(I2)') ipset
378 nprm=
'EPS09NLO,'//chpset
379 CALL setlhaparm(nprm)
381 CALL structa(xx,qq,inumod,upv,dnv,usea,dsea,
str,chm,bot,
394 upv=
z/a*upvp+(a-
z)/a*upvn
395 usea=
z/a*useap+(a-
z)/a*usean
396 dnv=
z/a*dnvp+(a-
z)/a*dnvn
397 dsea=
z/a*dseap+(a-
z)/a*dsean
399 CALL
structm(xx,qq,upv,dnv,usea,dsea,
str,chm,bot,
top,glu)
419 WRITE(mstu(11),5200) kf,
mstp(52),
mstp(51)
424 IF(kfa.EQ.111.OR.(kfa.EQ.22.AND.
mint(109).EQ.2))
THEN
425 IF(kfa.EQ.22.AND.
mstp(55).GE.5.AND.
mstp(55).LE.12)
THEN
430 xps=0.5d0*(xpq(1)+xpq(-2))
431 xpv=0.5d0*(xpq(2)+xpq(-1))-xps
435 xpvl=0.5d0*(xpval(1)+xpval(2)+xpval(-1)+xpval(-2))+
436 & xpval(3)+xpval(4)+xpval(5)
440 IF(kfa.EQ.22.AND.
mint(105).LE.223)
THEN
441 xpq(1)=xpq(1)+0.2d0*xpv
442 xpq(2)=xpq(2)+0.8d0*xpv
445 ELSEIF(kfa.EQ.22.AND.
mint(105).EQ.333)
THEN
448 ELSEIF(kfa.EQ.22.AND.
mint(105).EQ.443)
THEN
451 IF(
mstp(55).GE.9)
THEN
457 xpq(1)=xpq(1)+0.5d0*xpv
458 xpq(2)=xpq(2)+0.5d0*xpv
464 xpval(-kfl)=xpval(kfl)
469 IF(kfa.EQ.22.AND.
mint(109).EQ.2.AND..NOT.(
mstp(56).EQ.1
470 & .AND.
mstp(55).GE.5.AND.
mstp(55).LE.12))
THEN
472 xpq(kfl)=
vint(281)*xpq(kfl)
473 xpval(kfl)=
vint(281)*xpval(kfl)
479 ELSEIF(kfa.EQ.321)
THEN
480 xpq(-3)=xpq(-3)+xpq(-1)-xpq(1)
484 ELSEIF(kfa.EQ.130.OR.kfa.EQ.310)
THEN
485 xps=0.5d0*(xpq(1)+xpq(-2))
486 xpv=0.5d0*(xpq(2)+xpq(-1))-xps
489 xpq(1)=xpq(1)+0.5d0*xpv
490 xpq(-1)=xpq(-1)+0.5d0*xpv
491 xpq(3)=xpq(3)+0.5d0*xpv
492 xpq(-3)=xpq(-3)+0.5d0*xpv
493 xpv=0.5d0*(xpval(2)+xpval(-1))
502 ELSEIF(kfa.EQ.2112)
THEN
514 ELSEIF(kfa.EQ.3122.OR.kfa.EQ.3112.OR.kfa.EQ.3212.OR.kfa.EQ.3222
515 & .OR.kfa.EQ.3312.OR.kfa.EQ.3322.OR.kfa.EQ.3334)
THEN
516 xpv=(xpq(1)+xpq(2)-xpq(-1)-xpq(-2))/3d0
517 xps=0.5d0*(xpq(-1)+xpq(-2))
522 xpq(kfa/1000)=xpq(kfa/1000)+xpv
523 xpq(mod(kfa/100,10))=xpq(mod(kfa/100,10))+xpv
524 xpq(mod(kfa/10,10))=xpq(mod(kfa/10,10))+xpv
525 xpv=(xpval(1)+xpval(2))/3d0
528 xpval(kfa/1000)=xpval(kfa/1000)+xpv
529 xpval(mod(kfa/100,10))=xpval(mod(kfa/100,10))+xpv
530 xpval(mod(kfa/10,10))=xpval(mod(kfa/10,10))+xpv
536 IF(kfl.EQ.21.OR.kfl.EQ.22.OR.kfl.EQ.23.OR.kfl.EQ.25) goto 290
543 xpval(kfl)=xpval(-kfl)
553 IF ((js.EQ.1.OR.js.EQ.2).AND.
mint(35).GE.2)
THEN
557 IF(kfvsel.NE.0.AND.(kfa.EQ.111.OR.kfa.EQ.22))
THEN
561 xpq(kfl)=
max(0d0,xpq(kfl)-xpval(kfl))
564 xpq(iabs(kfvsel))=xpq(iabs(kfvsel))+xpvl
565 xpval(iabs(kfvsel))=xpvl
568 xpval(-kfl)=xpval(kfl)
573 ELSEIF(kfvsel.NE.0.AND.(kfa.EQ.130.OR.kfa.EQ.310))
THEN
575 IF(kfvsel.EQ.-1.OR.kfvsel.EQ.3) kfs=-1
576 xpq(kfs)=xpq(kfs)+xpval(-kfs)
577 xpval(kfs)=xpval(kfs)+xpval(-kfs)
578 xpq(-kfs)=
max(0d0,xpq(-kfs)-xpval(-kfs))
581 xpq(kfs)=xpq(kfs)+xpval(-kfs)
582 xpval(kfs)=xpval(kfs)+xpval(-kfs)
583 xpq(-kfs)=
max(0d0,xpq(-kfs)-xpval(-kfs))
596 IF(ifl.EQ.0) goto 350
600 IF(kfival(js,1).EQ.ifl) ivorg=ivorg+1
601 IF(kfival(js,2).EQ.ifl) ivorg=ivorg+1
602 IF(kfival(js,3).EQ.ifl) ivorg=ivorg+1
605 IF(kfival(js,1).EQ.0.AND.iabs(ifl).EQ.1) ivorg=1
610 IF (
i1.EQ.
mint(36)) goto 330
611 IF (
k(imi(js,
i1,1),2).EQ.ifl.AND.imi(js,
i1,2).EQ.0)
617 sea=
max(0d0,xpq(ifl)-val)
622 IF (ivorg.NE.0.AND.ivrem.NE.ivorg) xpsvc(ifl,0)=
629 IF(kfa.EQ.2212.OR.kfa.EQ.2112)
THEN
631 IF (kfa.EQ.2112) iaflp=3-iaflp
634 ELSEIF(kfa.GT.1000)
THEN
635 vpavg=(pavg(1,q2)+2d0*pavg(2,q2))/3d0
639 vpavg=0.5d0*(pavg(1,q2)+2d0*pavg(2,q2))
641 pvctot(js,-1)=pvctot(js,-1)+ivorg*vpavg
642 pvctot(js, 0)=pvctot(js, 0)+(ivorg-ivrem)*vpavg
649 DO 340 ivc=1,nvc(js,ifl)
651 IF (xassoc(js,ifl,ivc).LE.0d0)
THEN
657 xs=xassoc(js,ifl,ivc)
665 xcfac=(xrem+xs)/xrem*cmpfac
666 pvctot(js,1)=pvctot(js,1)+xcfac*
pypcmp(ys/cmpfac,
mstp(87))
674 rsfac=1d0+(pvctot(js,0)-pvctot(js,1))/(1d0-pvctot(js,-1))
675 IF (rsfac.LE.0d0)
THEN
677 IF (nresc.EQ.1) cmpfac =
678 & (1d0-(pvctot(js,-1)-pvctot(js,0)))/pvctot(js,1)
682 IF (nresc.LE.10) goto 345
684 &
'(PYPDFU:) Negative reshaping factor persists!')
685 WRITE(mstu(11),5300) (pvctot(js,itmp),itmp=-1,1), rsfac
689 xpsvc(ifl,-1)=rsfac*xpsvc(ifl,-1)
692 DO 360 isvc=-1,nvc(js,ifl)
693 xpq(ifl)=xpq(ifl)+xpsvc(ifl,isvc)
706 xpq(kfl)=
max(0d0,xpq(kfl))
707 IF(iabs(kfl).GT.
mstp(58).AND.iabs(kfl).LE.8) xpq(kfl)=0d0
711 5000
FORMAT(
' Error: x value outside physical range; x =',1
p,d12.3)
712 5100
FORMAT(
' Error: illegal particle code for parton distribution;',
714 5200
FORMAT(
' Error: unknown parton distribution; KF, library, set =',
716 5300
FORMAT(
' Original valence momentum fraction : ',f6.3/
717 &
' Removed valence momentum fraction : ',f6.3/
718 &
' Added companion momentum fraction : ',f6.3/
719 &
' Resulting rescale factor : ',f6.3)