Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pypdfu.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pypdfu.f
1 
2 C*********************************************************************
3 
4 C...PYPDFU
5 C...Gives electron, muon, tau, photon, pi+, neutron, proton and hyperon
6 C...parton distributions according to a few different parametrizations.
7 C...Note that what is coded is x times the probability distribution,
8 C...i.e. xq(x,Q2) etc.
9 
10  SUBROUTINE pypdfu(KF,X,Q2,XPQ)
11 
12 C...Double precision and integer declarations.
13  IMPLICIT DOUBLE PRECISION(a-h, o-z)
14  IMPLICIT INTEGER(i-n)
15  INTEGER pyk,pychge,pycomp
16 C...Commonblocks.
17  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
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)
20  common/pypars/mstp(200),parp(200),msti(200),pari(200)
21  common/pyint1/mint(400),vint(400)
22  common/pyint8/xpvmd(-6:6),xpanl(-6:6),xpanh(-6:6),xpbeh(-6:6),
23  &xpdir(-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)
28  SAVE /pyjets/,/pydat1/,/pydat2/,/pypars/,/pyint1/,/pyint8/,
29  &/pyint9/,/pyintm/
30 C...Local arrays.
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)
33  SAVE ppar
34 
35 C...Interface to PDFLIB.
36  common/w50513/xmin,xmax,q2min,q2max
37  SAVE /w50513/
38  DOUBLE PRECISION xx,qq,upv,dnv,usea,dsea,str,chm,bot,top,glu,
39  &value(20),xmin,xmax,q2min,q2max
40  CHARACTER*20 parm(20)
41  DATA value/20*0d0/,parm/20*' '/
42 
43 C...Added by liang 1/6/12
44 C...Switches for nuclear PDF correction
45  common/pynucl/inumod,chanum,order
46  SAVE /pynucl/
47  DOUBLE PRECISION inumod,chanum
48  INTEGER order
49  CHARACTER*20 chpset
50  CHARACTER*20 nprm
51 
52 
53 C...Data related to Schuler-Sjostrand photon distributions.
54  DATA alamga/0.2d0/, pmcga/1.3d0/, pmbga/4.6d0/
55 
56 C...Valence PDF momentum integral parametrizations PER PARTON!
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)*
60  &log(log(max(q2,1d0)/0.04d0)))
61 
62 C...Reset parton distributions.
63  mint(92)=0
64  DO 100 kfl=-25,25
65  xpq(kfl)=0d0
66  100 CONTINUE
67  DO 110 kfl=-6,6
68  xpval(kfl)=0d0
69  110 CONTINUE
70 
71 C...Check x and particle species.
72  IF(x.LE.0d0.OR.x.GE.1d0) THEN
73  WRITE(mstu(11),5000) x
74  goto 9999
75  ENDIF
76  kfa=iabs(kf)
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
83  goto 9999
84  ENDIF
85 
86 C...Electron (or muon or tau) parton distribution call.
87  IF(kfa.EQ.11.OR.kfa.EQ.13.OR.kfa.EQ.15) THEN
88  CALL pypdel(kfa,x,q2,xpel)
89  DO 120 kfl=-25,25
90  xpq(kfl)=xpel(kfl)
91  120 CONTINUE
92 
93 C...Photon parton distribution call (VDM+anomalous).
94  ELSEIF(kfa.EQ.22.AND.mint(109).LE.1) THEN
95  IF(mstp(56).EQ.1.AND.mstp(55).EQ.1) THEN
96  CALL pypdga(x,q2,xpga)
97  DO 130 kfl=-6,6
98  xpq(kfl)=xpga(kfl)
99  130 CONTINUE
100  xpvu=4d0*(xpq(2)-xpq(1))/3d0
101  xpval(1)=xpvu/4d0
102  xpval(2)=xpvu
103  xpval(3)=min(xpq(3),xpvu/4d0)
104  xpval(4)=min(xpq(4),xpvu)
105  xpval(5)=min(xpq(5),xpvu/4d0)
106  xpval(-1)=xpval(1)
107  xpval(-2)=xpval(2)
108  xpval(-3)=xpval(3)
109  xpval(-4)=xpval(4)
110  xpval(-5)=xpval(5)
111  ELSEIF(mstp(56).EQ.1.AND.mstp(55).GE.5.AND.mstp(55).LE.8) THEN
112  q2mx=q2
113  p2mx=0.36d0
114  IF(mstp(55).GE.7) p2mx=4.0d0
115  IF(mstp(57).EQ.0) q2mx=p2mx
116  p2=0d0
117  IF(vint(120).LT.0d0) p2=vint(120)**2
118  CALL pyggam(mstp(55)-4,x,q2mx,p2,mstp(60),f2gam,xpga)
119  DO 140 kfl=-6,6
120  xpq(kfl)=xpga(kfl)
121  xpval(kfl)=vxpdgm(kfl)
122  140 CONTINUE
123  vint(231)=p2mx
124  ELSEIF(mstp(56).EQ.1.AND.mstp(55).GE.9.AND.mstp(55).LE.12) THEN
125  q2mx=q2
126  p2mx=0.36d0
127  IF(mstp(55).GE.11) p2mx=4.0d0
128  IF(mstp(57).EQ.0) q2mx=p2mx
129  p2=0d0
130  IF(vint(120).LT.0d0) p2=vint(120)**2
131  CALL pyggam(mstp(55)-8,x,q2mx,p2,mstp(60),f2gam,xpga)
132  DO 150 kfl=-6,6
133  xpq(kfl)=xpvmd(kfl)+xpanl(kfl)+xpbeh(kfl)+xpdir(kfl)
134  xpval(kfl)=vxpvmd(kfl)+vxpanl(kfl)+xpbeh(kfl)+xpdir(kfl)
135  150 CONTINUE
136  vint(231)=p2mx
137  ELSEIF(mstp(56).EQ.2) THEN
138 C...Call PDFLIB parton distributions.
139  parm(1)='NPTYPE'
140  value(1)=3
141  parm(2)='NGROUP'
142  value(2)=mstp(55)/1000
143  parm(3)='NSET'
144  value(3)=mod(mstp(55),1000)
145  IF(mint(93).NE.3000000+mstp(55)) THEN
146  CALL pdfset(parm,value)
147  mint(93)=3000000+mstp(55)
148  ENDIF
149  xx=x
150  qq2=max(0d0,q2min,q2)
151  IF(mstp(57).EQ.0) qq2=q2min
152  p2=0d0
153  IF(vint(120).LT.0d0) p2=vint(120)**2
154  ip2=mstp(60)
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,
161  & bot,top,glu)
162  ELSE
163  upv=0d0
164  dnv=0d0
165  usea=0d0
166  dsea=0d0
167  str=0d0
168  chm=0d0
169  bot=0d0
170  top=0d0
171  glu=0d0
172  ENDIF
173  ELSE
174  IF(p2.LT.qq2) THEN
175  CALL structp(xx,qq2,p2,ip2,upv,dnv,usea,dsea,str,chm,
176  & bot,top,glu)
177  ELSE
178  upv=0d0
179  dnv=0d0
180  usea=0d0
181  dsea=0d0
182  str=0d0
183  chm=0d0
184  bot=0d0
185  top=0d0
186  glu=0d0
187  ENDIF
188  ENDIF
189  vint(231)=q2min
190  xpq(0)=glu
191  xpq(1)=dnv
192  xpq(-1)=dnv
193  xpq(2)=upv
194  xpq(-2)=upv
195  xpq(3)=str
196  xpq(-3)=str
197  xpq(4)=chm
198  xpq(-4)=chm
199  xpq(5)=bot
200  xpq(-5)=bot
201  xpq(6)=top
202  xpq(-6)=top
203  xpvu=4d0*(xpq(2)-xpq(1))/3d0
204  xpval(1)=xpvu/4d0
205  xpval(2)=xpvu
206  xpval(3)=min(xpq(3),xpvu/4d0)
207  xpval(4)=min(xpq(4),xpvu)
208  xpval(5)=min(xpq(5),xpvu/4d0)
209  xpval(-1)=xpval(1)
210  xpval(-2)=xpval(2)
211  xpval(-3)=xpval(3)
212  xpval(-4)=xpval(4)
213  xpval(-5)=xpval(5)
214  ELSE
215  WRITE(mstu(11),5200) kf,mstp(56),mstp(55)
216  ENDIF
217 
218 C...Pion/gammaVDM parton distribution call.
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)
224  q2mx=q2
225  p2mx=0.36d0
226  IF(iset.GE.3) p2mx=4.0d0
227  IF(mstp(57).EQ.0) q2mx=p2mx
228  p2=0d0
229  IF(vint(120).LT.0d0) p2=vint(120)**2
230  CALL pyggam(iset,x,q2mx,p2,mstp(60),f2gam,xpga)
231  DO 160 kfl=-6,6
232  xpq(kfl)=xpvmd(kfl)
233  xpval(kfl)=vxpvmd(kfl)
234  160 CONTINUE
235  vint(231)=p2mx
236  ELSEIF(mstp(54).EQ.1.AND.mstp(53).GE.1.AND.mstp(53).LE.3) THEN
237  CALL pypdpi(x,q2,xppi)
238  DO 170 kfl=-6,6
239  xpq(kfl)=xppi(kfl)
240  170 CONTINUE
241  xpval(2)=xpq(2)-xpq(-2)
242  xpval(-1)=xpq(-1)-xpq(1)
243  ELSEIF(mstp(54).EQ.2) THEN
244 C...Call PDFLIB parton distributions.
245  parm(1)='NPTYPE'
246  value(1)=2
247  parm(2)='NGROUP'
248  value(2)=mstp(53)/1000
249  parm(3)='NSET'
250  value(3)=mod(mstp(53),1000)
251  IF(mint(93).NE.2000000+mstp(53)) THEN
252  CALL pdfset(parm,value)
253  mint(93)=2000000+mstp(53)
254  ENDIF
255  xx=x
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)
259  vint(231)=q2min
260  xpq(0)=glu
261  xpq(1)=dsea
262  xpq(-1)=upv+dsea
263  xpq(2)=upv+usea
264  xpq(-2)=usea
265  xpq(3)=str
266  xpq(-3)=str
267  xpq(4)=chm
268  xpq(-4)=chm
269  xpq(5)=bot
270  xpq(-5)=bot
271  xpq(6)=top
272  xpq(-6)=top
273  xpval(2)=upv
274  xpval(-1)=upv
275  ELSE
276  WRITE(mstu(11),5200) kf,mstp(54),mstp(53)
277  ENDIF
278 
279 C...Anomalous photon parton distribution call.
280  ELSEIF(kfa.EQ.22.AND.mint(109).EQ.3) THEN
281  q2mx=q2
282  p2mx=parp(15)**2
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
287  p2=0d0
288  IF(vint(120).LT.0d0) p2=vint(120)**2
289  CALL pyggam(mstp(55)-4,x,q2mx,p2,mstp(60),f2gm,xpga)
290  DO 180 kfl=-6,6
291  xpq(kfl)=xpanl(kfl)+xpanh(kfl)
292  xpval(kfl)=vxpanl(kfl)+vxpanh(kfl)
293  180 CONTINUE
294  vint(231)=p2mx
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
299  p2=0d0
300  IF(vint(120).LT.0d0) p2=vint(120)**2
301  CALL pyggam(mstp(55)-8,x,q2mx,p2,mstp(60),f2gm,xpga)
302  DO 190 kfl=-6,6
303  xpq(kfl)=max(0d0,xpanl(kfl)+xpbeh(kfl)+xpdir(kfl))
304  xpval(kfl)=max(0d0,vxpanl(kfl)+xpbeh(kfl)+xpdir(kfl))
305  190 CONTINUE
306  vint(231)=p2mx
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)
310  DO 200 kfl=-6,6
311  xpq(kfl)=xpga(kfl)
312  xpval(kfl)=vxpga(kfl)
313  200 CONTINUE
314  vint(231)=p2mx
315  ELSEIF(mstp(55).GE.1.AND.mstp(55).LE.5) THEN
316  IF(mstp(57).EQ.0) q2mx=p2mx
317  CALL pygvmd(0,mstp(55),x,q2mx,p2mx,parp(1),xpga,vxpga)
318  DO 210 kfl=-6,6
319  xpq(kfl)=xpga(kfl)
320  xpval(kfl)=vxpga(kfl)
321  210 CONTINUE
322  vint(231)=p2mx
323  ELSE
324  220 rkf=11d0*pyr(0)
325  kfr=1
326  IF(rkf.GT.1d0) kfr=2
327  IF(rkf.GT.5d0) kfr=3
328  IF(rkf.GT.6d0) kfr=4
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
333  CALL pygvmd(0,kfr,x,q2mx,p2mx,parp(1),xpga,vxpga)
334  DO 230 kfl=-6,6
335  xpq(kfl)=xpga(kfl)
336  xpval(kfl)=vxpga(kfl)
337  230 CONTINUE
338  vint(231)=p2mx
339  ENDIF
340 
341 C...Proton parton distribution call.
342  ELSE
343  IF(mstp(52).EQ.1.AND.mstp(51).GE.1.AND.mstp(51).LE.20) THEN
344  CALL pypdpr(x,q2,xppr)
345  DO 240 kfl=-6,6
346  xpq(kfl)=xppr(kfl)
347  240 CONTINUE
348  xpval(1)=xpq(1)-xpq(-1)
349  xpval(2)=xpq(2)-xpq(-2)
350  ELSEIF(mstp(52).EQ.2) THEN
351 C...Call PDFLIB parton distributions.
352  parm(1)='NPTYPE'
353  value(1)=1
354  parm(2)='NGROUP'
355  value(2)=mstp(51)/1000
356  parm(3)='NSET'
357  value(3)=mod(mstp(51),1000)
358  IF(mint(93).NE.1000000+mstp(51)) THEN
359 C...Uncommented by liang to link LHAPDF 1/6/12
360  CALL pdfset(parm,value)
361  mint(93)=1000000+mstp(51)
362  ENDIF
363  xx=x
364  qq=sqrt(max(0d0,q2min,q2))
365  IF(mstp(57).EQ.0) qq=sqrt(q2min)
366 C...Modified by liang to add nuclear correction 1/6/12*************
367  IF(inumod.GT.1d0) THEN
368  IF(order/100.EQ.1) THEN
369 C print*,'Now run LO nuclear pdf'
370  ipset=order-100
371  write(chpset,'(I2)') ipset
372  nprm='EPS09LO,'//chpset
373  CALL setlhaparm(nprm)
374  ELSE IF(order/100.EQ.2) THEN
375 C print*,'Now run NLO nuclear pdf'
376  ipset=order-200
377  write(chpset,'(I2)') ipset
378  nprm='EPS09NLO,'//chpset
379  CALL setlhaparm(nprm)
380  ENDIF
381  CALL structa(xx,qq,inumod,upv,dnv,usea,dsea,str,chm,bot,
382  & top,glu)
383 C...Isospin transformation for nucleus A with Z protons************
384  a=inumod
385  z=chanum
386  upvp=upv
387  upvn=dnv
388  useap=usea
389  usean=dsea
390  dnvp=dnv
391  dnvn=upv
392  dseap=dsea
393  dsean=usea
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
398  ELSE
399  CALL structm(xx,qq,upv,dnv,usea,dsea,str,chm,bot,top,glu)
400  ENDIF
401 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
402  vint(231)=q2min
403  xpq(0)=glu
404  xpq(1)=dnv+dsea
405  xpq(-1)=dsea
406  xpq(2)=upv+usea
407  xpq(-2)=usea
408  xpq(3)=str
409  xpq(-3)=str
410  xpq(4)=chm
411  xpq(-4)=chm
412  xpq(5)=bot
413  xpq(-5)=bot
414  xpq(6)=top
415  xpq(-6)=top
416  xpval(1)=dnv
417  xpval(2)=upv
418  ELSE
419  WRITE(mstu(11),5200) kf,mstp(52),mstp(51)
420  ENDIF
421  ENDIF
422 
423 C...Isospin average for pi0/gammaVDM.
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
426  xpv=xpq(2)-xpq(1)
427  xpq(2)=xpq(1)
428  xpq(-2)=xpq(-1)
429  ELSE
430  xps=0.5d0*(xpq(1)+xpq(-2))
431  xpv=0.5d0*(xpq(2)+xpq(-1))-xps
432  xpq(2)=xps
433  xpq(-1)=xps
434  ENDIF
435  xpvl=0.5d0*(xpval(1)+xpval(2)+xpval(-1)+xpval(-2))+
436  & xpval(3)+xpval(4)+xpval(5)
437  DO 250 kfl=-6,6
438  xpval(kfl)=0d0
439  250 CONTINUE
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
443  xpval(1)=0.2d0*xpvl
444  xpval(2)=0.8d0*xpvl
445  ELSEIF(kfa.EQ.22.AND.mint(105).EQ.333) THEN
446  xpq(3)=xpq(3)+xpv
447  xpval(3)=xpvl
448  ELSEIF(kfa.EQ.22.AND.mint(105).EQ.443) THEN
449  xpq(4)=xpq(4)+xpv
450  xpval(4)=xpvl
451  IF(mstp(55).GE.9) THEN
452  DO 260 kfl=-6,6
453  xpq(kfl)=0d0
454  260 CONTINUE
455  ENDIF
456  ELSE
457  xpq(1)=xpq(1)+0.5d0*xpv
458  xpq(2)=xpq(2)+0.5d0*xpv
459  xpval(1)=0.5d0*xpvl
460  xpval(2)=0.5d0*xpvl
461  ENDIF
462  DO 270 kfl=1,6
463  xpq(-kfl)=xpq(kfl)
464  xpval(-kfl)=xpval(kfl)
465  270 CONTINUE
466 
467 C...Rescale for gammaVDM by effective gamma -> rho coupling.
468 C+++Do not rescale?
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
471  DO 280 kfl=-6,6
472  xpq(kfl)=vint(281)*xpq(kfl)
473  xpval(kfl)=vint(281)*xpval(kfl)
474  280 CONTINUE
475  vint(232)=vint(281)*xpv
476  ENDIF
477 
478 C...Simple recipes for kaons.
479  ELSEIF(kfa.EQ.321) THEN
480  xpq(-3)=xpq(-3)+xpq(-1)-xpq(1)
481  xpq(-1)=xpq(1)
482  xpval(-3)=xpval(-1)
483  xpval(-1)=0d0
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
487  xpq(2)=xps
488  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))
494  xpval(2)=0d0
495  xpval(-1)=0d0
496  xpval(1)=0.5d0*xpv
497  xpval(-1)=0.5d0*xpv
498  xpval(3)=0.5d0*xpv
499  xpval(-3)=0.5d0*xpv
500 
501 C...Isospin conjugation for neutron.
502  ELSEIF(kfa.EQ.2112) THEN
503  xpsv=xpq(1)
504  xpq(1)=xpq(2)
505  xpq(2)=xpsv
506  xpsv=xpq(-1)
507  xpq(-1)=xpq(-2)
508  xpq(-2)=xpsv
509  xpsv=xpval(1)
510  xpval(1)=xpval(2)
511  xpval(2)=xpsv
512 
513 C...Simple recipes for hyperon (average valence parton distribution).
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))
518  xpq(1)=xps
519  xpq(2)=xps
520  xpq(-1)=xps
521  xpq(-2)=xps
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
526  xpval(1)=0d0
527  xpval(2)=0d0
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
531  ENDIF
532 
533 C...Charge conjugation for antiparticle.
534  IF(kf.LT.0) THEN
535  DO 290 kfl=1,25
536  IF(kfl.EQ.21.OR.kfl.EQ.22.OR.kfl.EQ.23.OR.kfl.EQ.25) goto 290
537  xpsv=xpq(kfl)
538  xpq(kfl)=xpq(-kfl)
539  xpq(-kfl)=xpsv
540  290 CONTINUE
541  DO 300 kfl=1,6
542  xpsv=xpval(kfl)
543  xpval(kfl)=xpval(-kfl)
544  xpval(-kfl)=xpsv
545  300 CONTINUE
546  ENDIF
547 
548 C...MULTIPLE INTERACTIONS - PDF RESHAPING.
549 C...Set side.
550  js=mint(30)
551 C...Only reshape PDFs for the non-first interactions;
552 C...But need valence/sea separation already from first interaction.
553  IF ((js.EQ.1.OR.js.EQ.2).AND.mint(35).GE.2) THEN
554  kfvsel=kfival(js,1)
555 C...If valence quark kicked out of pi0 or gamma then that decides
556 C...whether we should consider state as d dbar, u ubar, s sbar, etc.
557  IF(kfvsel.NE.0.AND.(kfa.EQ.111.OR.kfa.EQ.22)) THEN
558  xpvl=0d0
559  DO 310 kfl=1,6
560  xpvl=xpvl+xpval(kfl)
561  xpq(kfl)=max(0d0,xpq(kfl)-xpval(kfl))
562  xpval(kfl)=0d0
563  310 CONTINUE
564  xpq(iabs(kfvsel))=xpq(iabs(kfvsel))+xpvl
565  xpval(iabs(kfvsel))=xpvl
566  DO 320 kfl=1,6
567  xpq(-kfl)=xpq(kfl)
568  xpval(-kfl)=xpval(kfl)
569  320 CONTINUE
570 
571 C...If valence quark kicked out of K0S or K0S then that decides whether
572 C...we should consider state as d sbar or s dbar.
573  ELSEIF(kfvsel.NE.0.AND.(kfa.EQ.130.OR.kfa.EQ.310)) THEN
574  kfs=1
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))
579  xpval(-kfs)=0d0
580  kfs=-3*kfs
581  xpq(kfs)=xpq(kfs)+xpval(-kfs)
582  xpval(kfs)=xpval(kfs)+xpval(-kfs)
583  xpq(-kfs)=max(0d0,xpq(-kfs)-xpval(-kfs))
584  xpval(-kfs)=0d0
585  ENDIF
586 
587 C...XPQ distributions are nominal for a (signed) beam particle
588 C...of KF type, with 1-Sum(x_prev) rescaled to 1.
589  cmpfac=1d0
590  nresc=0
591  345 nresc=nresc+1
592  pvctot(js,-1)=0d0
593  pvctot(js, 0)=0d0
594  pvctot(js, 1)=0d0
595  DO 350 ifl=-6,6
596  IF(ifl.EQ.0) goto 350
597 
598 C...Count up number of original IFL valence quarks.
599  ivorg=0
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
603 C...For pi0/gamma/K0S/K0L without valence flavour decided yet, here
604 C...bookkeep as if d dbar (for total momentum sum in valence sector).
605  IF(kfival(js,1).EQ.0.AND.iabs(ifl).EQ.1) ivorg=1
606 C...Count down number of remaining IFL valence quarks. Skip current
607 C...interaction initiator.
608  ivrem=ivorg
609  DO 330 i1=1,nmi(js)
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)
612  & ivrem=ivrem-1
613  330 CONTINUE
614 
615 C...Separate out original VALENCE and SEA content.
616  val=xpval(ifl)
617  sea=max(0d0,xpq(ifl)-val)
618  xpsvc(ifl,0)=val
619  xpsvc(ifl,-1)=sea
620 
621 C...Rescale valence content if changed.
622  IF (ivorg.NE.0.AND.ivrem.NE.ivorg) xpsvc(ifl,0)=
623  & (val*ivrem)/ivorg
624 
625 C...Momentum integrals of original and removed valence quarks.
626  IF(ivorg.NE.0) THEN
627 C...For p/n/pbar/nbar beams can split into d_val and u_val.
628 C...Isospin conjugation for neutrons
629  IF(kfa.EQ.2212.OR.kfa.EQ.2112) THEN
630  iaflp=iabs(ifl)
631  IF (kfa.EQ.2112) iaflp=3-iaflp
632  vpavg=pavg(iaflp,q2)
633 C...For other baryons average d_val and u_val, like for PDFs.
634  ELSEIF(kfa.GT.1000) THEN
635  vpavg=(pavg(1,q2)+2d0*pavg(2,q2))/3d0
636 C...For mesons and photon average d_val and u_val and scale by 3/2.
637 C...Very crude, especially for photon.
638  ELSE
639  vpavg=0.5d0*(pavg(1,q2)+2d0*pavg(2,q2))
640  ENDIF
641  pvctot(js,-1)=pvctot(js,-1)+ivorg*vpavg
642  pvctot(js, 0)=pvctot(js, 0)+(ivorg-ivrem)*vpavg
643  ENDIF
644 
645 C...Now add companions (at X with partner having been at Z=XASSOC).
646 C...NOTE: due to the assumed simple x scaling, the partner was at what
647 C...corresponds to a higher Z than XASSOC, if there were intermediate
648 C...scatterings. Nothing done about that for the moment.
649  DO 340 ivc=1,nvc(js,ifl)
650 C...Skip companions that have been kicked out
651  IF (xassoc(js,ifl,ivc).LE.0d0) THEN
652  xpsvc(ifl,ivc)=0d0
653  goto 340
654  ELSE
655 C...Momentum fraction of the partner quark.
656 C...Use rescaled YS = XS/(1-Sum_rest) where X and XS are not in "rest".
657  xs=xassoc(js,ifl,ivc)
658  xrem=vint(142+js)
659  ys=xs/(xrem+xs)
660 C...Momentum fraction of the companion quark.
661 C...Rescale from X = x/XREM to Y = x/(1-Sum_rest) -> factor (1-YS).
662  y=x*(1d0-ys)
663  xpsvc(ifl,ivc)=pyfcmp(y/cmpfac,ys/cmpfac,mstp(87))
664 C...Add to momentum sum, with rescaling compensation factor.
665  xcfac=(xrem+xs)/xrem*cmpfac
666  pvctot(js,1)=pvctot(js,1)+xcfac*pypcmp(ys/cmpfac,mstp(87))
667  ENDIF
668  340 CONTINUE
669  350 CONTINUE
670 
671 C...Wait until all flavours treated, then rescale seas and gluon.
672  xpsvc(0,-1)=xpq(0)
673  xpsvc(0,0)=0d0
674  rsfac=1d0+(pvctot(js,0)-pvctot(js,1))/(1d0-pvctot(js,-1))
675  IF (rsfac.LE.0d0) THEN
676 C...First calculate factor needed to exactly restore pz cons.
677  IF (nresc.EQ.1) cmpfac =
678  & (1d0-(pvctot(js,-1)-pvctot(js,0)))/pvctot(js,1)
679 C...Add a bit of headroom
680  cmpfac=0.99*cmpfac
681 C...Try a few times if more headroom is needed, then print error message.
682  IF (nresc.LE.10) goto 345
683  CALL pyerrm(15,
684  & '(PYPDFU:) Negative reshaping factor persists!')
685  WRITE(mstu(11),5300) (pvctot(js,itmp),itmp=-1,1), rsfac
686  rsfac=0d0
687  ENDIF
688  DO 370 ifl=-6,6
689  xpsvc(ifl,-1)=rsfac*xpsvc(ifl,-1)
690 C...Also store resulting distributions in XPQ
691  xpq(ifl)=0d0
692  DO 360 isvc=-1,nvc(js,ifl)
693  xpq(ifl)=xpq(ifl)+xpsvc(ifl,isvc)
694  360 CONTINUE
695  370 CONTINUE
696 C...Save companion reweighting factor for PYPTIS.
697  vint(140)=cmpfac
698  ENDIF
699 
700 
701 C...Allow gluon also in position 21.
702  xpq(21)=xpq(0)
703 
704 C...Check positivity and reset above maximum allowed flavour.
705  DO 380 kfl=-25,25
706  xpq(kfl)=max(0d0,xpq(kfl))
707  IF(iabs(kfl).GT.mstp(58).AND.iabs(kfl).LE.8) xpq(kfl)=0d0
708  380 CONTINUE
709 
710 C...Formats for error printouts.
711  5000 FORMAT(' Error: x value outside physical range; x =',1p,d12.3)
712  5100 FORMAT(' Error: illegal particle code for parton distribution;',
713  &' KF =',i5)
714  5200 FORMAT(' Error: unknown parton distribution; KF, library, set =',
715  &3i5)
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)
720 
721 C...Reset side pointer and return
722  9999 mint(30)=0
723 
724  RETURN
725  END