EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pystfu.F
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pystfu.F
1  SUBROUTINE pystfu(KF,X,Q2,XPQ)
2 
3 C...Gives electron, photon, pi+, neutron, proton and hyperon
4 C...structure functions according to a few different parametrizations.
5 C...Note that what is coded is x times the probability distribution,
6 C...i.e. xq(x,Q2) etc.
7  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
8  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
9  common/pypars/mstp(200),parp(200),msti(200),pari(200)
10  common/pyint1/mint(400),vint(400)
11  common/pyint8/xpvmd(-6:6),xpanl(-6:6),xpanh(-6:6),xpbeh(-6:6),
12  &xpdir(-6:6)
13  SAVE /ludat1/,/ludat2/
14  SAVE /pypars/,/pyint1/,/pyint8/
15  dimension xpq(-25:25),xpel(-25:25),xpga(-6:6),xppi(-6:6),
16  &xppr(-6:6)
17 
18 C...Interface to PDFLIB.
19  common/w50513/xmin,xmax,q2min,q2max
20  SAVE /w50513/
21  DOUBLE PRECISION xx,qq,upv,dnv,usea,dsea,str,chm,bot,top,glu,
22  &value(20),xmin,xmax,q2min,q2max
23  CHARACTER*20 parm(20)
24  DATA value/20*0d0/,parm/20*' '/
25 
26 C...Data related to Schuler-Sjostrand photon distributions.
27  DATA alamga/0.2/, pmcga/1.3/, pmbga/4.6/
28 
29 C...Reset structure functions.
30  mint(92)=0
31  DO 100 kfl=-25,25
32  xpq(kfl)=0.
33  100 CONTINUE
34 
35 C...Check x and particle species.
36  IF(x.LE.0..OR.x.GE.1.) THEN
37  WRITE(mstu(11),5000) x
38  RETURN
39  ENDIF
40  kfa=iabs(kf)
41  IF(kfa.NE.11.AND.kfa.NE.22.AND.kfa.NE.211.AND.kfa.NE.2112.AND.
42  &kfa.NE.2212.AND.kfa.NE.3122.AND.kfa.NE.3112.AND.kfa.NE.3212
43  &.AND.kfa.NE.3222.AND.kfa.NE.3312.AND.kfa.NE.3322.AND.
44  &kfa.NE.3334.AND.kfa.NE.111) THEN
45  WRITE(mstu(11),5100) kf
46  RETURN
47  ENDIF
48 
49 C...Electron structure function call.
50  IF(kfa.EQ.11) THEN
51  CALL pystel(x,q2,xpel)
52  DO 110 kfl=-25,25
53  xpq(kfl)=xpel(kfl)
54  110 CONTINUE
55 
56 C...Photon structure function call (VDM+anomalous).
57  ELSEIF(kfa.EQ.22.AND.mint(109).LE.1) THEN
58  IF(mstp(56).EQ.1.AND.mstp(55).EQ.1) THEN
59  CALL pystga(x,q2,xpga)
60  DO 120 kfl=-6,6
61  xpq(kfl)=xpga(kfl)
62  120 CONTINUE
63  ELSEIF(mstp(56).EQ.1.AND.mstp(55).GE.5.AND.mstp(55).LE.8) THEN
64  q2mx=q2
65  p2mx=0.36
66  IF(mstp(55).GE.7) p2mx=4.0
67  IF(mstp(57).EQ.0) q2mx=p2mx
68  CALL pyggam(mstp(55)-4,x,q2mx,0.,f2gam,xpga)
69  DO 130 kfl=-6,6
70  xpq(kfl)=xpga(kfl)
71  130 CONTINUE
72  vint(231)=p2mx
73  ELSEIF(mstp(56).EQ.1.AND.mstp(55).GE.9.AND.mstp(55).LE.12) THEN
74  q2mx=q2
75  p2mx=0.36
76  IF(mstp(55).GE.11) p2mx=4.0
77  IF(mstp(57).EQ.0) q2mx=p2mx
78  CALL pyggam(mstp(55)-8,x,q2mx,0.,f2gam,xpga)
79  DO 140 kfl=-6,6
80  xpq(kfl)=xpvmd(kfl)+xpanl(kfl)+xpbeh(kfl)+xpdir(kfl)
81  140 CONTINUE
82  vint(231)=p2mx
83  ELSEIF(mstp(56).EQ.2) THEN
84 C...Call PDFLIB structure functions.
85  parm(1)='NPTYPE'
86  value(1)=3
87  parm(2)='NGROUP'
88  value(2)=mstp(55)/1000
89  parm(3)='NSET'
90  value(3)=mod(mstp(55),1000)
91  IF(mint(93).NE.3000000+mstp(55)) THEN
92  CALL pdfset(parm,value)
93  mint(93)=3000000+mstp(55)
94  ENDIF
95  xx=x
96  qq=sqrt(max(0.,sngl(q2min),q2))
97  IF(mstp(57).EQ.0) qq=sqrt(q2min)
98  CALL structm(xx,qq,upv,dnv,usea,dsea,str,chm,bot,top,glu)
99  vint(231)=q2min
100  xpq(0)=glu
101  xpq(1)=dnv
102  xpq(-1)=dnv
103  xpq(2)=upv
104  xpq(-2)=upv
105  xpq(3)=str
106  xpq(-3)=str
107  xpq(4)=chm
108  xpq(-4)=chm
109  xpq(5)=bot
110  xpq(-5)=bot
111  xpq(6)=top
112  xpq(-6)=top
113  ELSE
114  WRITE(mstu(11),5200) kf,mstp(56),mstp(55)
115  ENDIF
116 
117 C...Pion/gammaVDM structure function call.
118  ELSEIF(kfa.EQ.211.OR.kfa.EQ.111.OR.(kfa.EQ.22.AND.
119  &mint(109).EQ.2)) THEN
120  IF(kfa.EQ.22.AND.mstp(56).EQ.1.AND.mstp(55).GE.5.AND.
121  & mstp(55).LE.12) THEN
122  iset=1+mod(mstp(55)-1,4)
123  q2mx=q2
124  p2mx=0.36
125  IF(iset.GE.3) p2mx=4.0
126  IF(mstp(57).EQ.0) q2mx=p2mx
127  CALL pygvmd(iset,2,x,q2mx,p2mx,alamga,xpga)
128  DO 150 kfl=-6,6
129  xpq(kfl)=xpga(kfl)
130  150 CONTINUE
131  vint(231)=p2mx
132  ELSEIF(mstp(54).EQ.1.AND.mstp(53).GE.1.AND.mstp(53).LE.3) THEN
133  CALL pystpi(x,q2,xppi)
134  DO 160 kfl=-6,6
135  xpq(kfl)=xppi(kfl)
136  160 CONTINUE
137  ELSEIF(mstp(54).EQ.2) THEN
138 C...Call PDFLIB structure functions.
139  parm(1)='NPTYPE'
140  value(1)=2
141  parm(2)='NGROUP'
142  value(2)=mstp(53)/1000
143  parm(3)='NSET'
144  value(3)=mod(mstp(53),1000)
145  IF(mint(93).NE.2000000+mstp(53)) THEN
146  CALL pdfset(parm,value)
147  mint(93)=2000000+mstp(53)
148  ENDIF
149  xx=x
150  qq=sqrt(max(0.,sngl(q2min),q2))
151  IF(mstp(57).EQ.0) qq=sqrt(q2min)
152  CALL structm(xx,qq,upv,dnv,usea,dsea,str,chm,bot,top,glu)
153  vint(231)=q2min
154  xpq(0)=glu
155  xpq(1)=dsea
156  xpq(-1)=upv+dsea
157  xpq(2)=upv+usea
158  xpq(-2)=usea
159  xpq(3)=str
160  xpq(-3)=str
161  xpq(4)=chm
162  xpq(-4)=chm
163  xpq(5)=bot
164  xpq(-5)=bot
165  xpq(6)=top
166  xpq(-6)=top
167  ELSE
168  WRITE(mstu(11),5200) kf,mstp(54),mstp(53)
169  ENDIF
170 
171 C...Anomalous photon structure function call.
172  ELSEIF(kfa.EQ.22.AND.mint(109).EQ.3) THEN
173  q2mx=q2
174  p2mx=parp(15)**2
175  IF(mstp(56).EQ.1.AND.mstp(55).LE.8) THEN
176  IF(mstp(55).EQ.5.OR.mstp(55).EQ.6) p2mx=0.36
177  IF(mstp(55).EQ.7.OR.mstp(55).EQ.8) p2mx=4.0
178  IF(mstp(57).EQ.0) q2mx=p2mx
179  CALL pygano(0,x,q2mx,p2mx,alamga,xpga)
180  DO 170 kfl=-6,6
181  xpq(kfl)=xpga(kfl)
182  170 CONTINUE
183  vint(231)=p2mx
184  ELSEIF(mstp(56).EQ.1) THEN
185  IF(mstp(55).EQ.9.OR.mstp(55).EQ.10) p2mx=0.36
186  IF(mstp(55).EQ.11.OR.mstp(55).EQ.12) p2mx=4.0
187  IF(mstp(57).EQ.0) q2mx=p2mx
188  CALL pyggam(mstp(55)-8,x,q2mx,0.,f2gm,xpga)
189  DO 180 kfl=-6,6
190  xpq(kfl)=max(0.,xpanl(kfl)+xpbeh(kfl)+xpdir(kfl))
191  180 CONTINUE
192  vint(231)=p2mx
193  ELSEIF(mstp(56).EQ.2) THEN
194  IF(mstp(57).EQ.0) q2mx=p2mx
195  CALL pygano(0,x,q2mx,p2mx,alamga,xpga)
196  DO 185 kfl=-6,6
197  xpq(kfl)=xpga(kfl)
198  185 CONTINUE
199  vint(231)=p2mx
200  ELSEIF(mstp(55).GE.1.AND.mstp(55).LE.5) THEN
201  IF(mstp(57).EQ.0) q2mx=p2mx
202  CALL pygvmd(0,mstp(55),x,q2mx,p2mx,parp(1),xpga)
203  DO 190 kfl=-6,6
204  xpq(kfl)=xpga(kfl)
205  190 CONTINUE
206  vint(231)=p2mx
207  ELSE
208  200 rkf=11.*rlu(0)
209  kfr=1
210  IF(rkf.GT.1.) kfr=2
211  IF(rkf.GT.5.) kfr=3
212  IF(rkf.GT.6.) kfr=4
213  IF(rkf.GT.10.) kfr=5
214  IF(kfr.EQ.4.AND.q2.LT.pmcga**2) goto 200
215  IF(kfr.EQ.5.AND.q2.LT.pmbga**2) goto 200
216  IF(mstp(57).EQ.0) q2mx=p2mx
217  CALL pygvmd(0,kfr,x,q2mx,p2mx,parp(1),xpga)
218  DO 210 kfl=-6,6
219  xpq(kfl)=xpga(kfl)
220  210 CONTINUE
221  vint(231)=p2mx
222  ENDIF
223 
224 C...Proton structure function call.
225  ELSE
226  IF(mstp(52).EQ.1.AND.mstp(51).GE.1.AND.mstp(51).LE.11) THEN
227  CALL pystpr(x,q2,xppr)
228  DO 220 kfl=-6,6
229  xpq(kfl)=xppr(kfl)
230  220 CONTINUE
231  ELSEIF(mstp(52).EQ.2) THEN
232 C...Call PDFLIB structure functions.
233  parm(1)='NPTYPE'
234  value(1)=1
235  parm(2)='NGROUP'
236  value(2)=mstp(51)/1000
237  parm(3)='NSET'
238  value(3)=mod(mstp(51),1000)
239  IF(mint(93).NE.1000000+mstp(51)) THEN
240  CALL pdfset(parm,value)
241  mint(93)=1000000+mstp(51)
242  ENDIF
243  xx=x
244  qq=sqrt(max(0.,sngl(q2min),q2))
245  IF(mstp(57).EQ.0) qq=sqrt(q2min)
246  CALL structm(xx,qq,upv,dnv,usea,dsea,str,chm,bot,top,glu)
247  vint(231)=q2min
248  xpq(0)=glu
249  xpq(1)=dnv+dsea
250  xpq(-1)=dsea
251  xpq(2)=upv+usea
252  xpq(-2)=usea
253  xpq(3)=str
254  xpq(-3)=str
255  xpq(4)=chm
256  xpq(-4)=chm
257  xpq(5)=bot
258  xpq(-5)=bot
259  xpq(6)=top
260  xpq(-6)=top
261  ELSE
262  WRITE(mstu(11),5200) kf,mstp(52),mstp(51)
263  ENDIF
264  ENDIF
265 
266 C...Isospin average for pi0/gammaVDM.
267  IF(kfa.EQ.111.OR.(kfa.EQ.22.AND.mint(109).EQ.2)) THEN
268  IF(kfa.EQ.22.AND.mstp(55).GE.5.AND.mstp(55).LE.12) THEN
269  xpv=xpq(2)-xpq(1)
270  xpq(2)=xpq(1)
271  xpq(-2)=xpq(-1)
272  ELSE
273  xps=0.5*(xpq(1)+xpq(-2))
274  xpv=0.5*(xpq(2)+xpq(-1))-xps
275  xpq(2)=xps
276  xpq(-1)=xps
277  ENDIF
278  IF(kfa.EQ.22.AND.mint(105).LE.223) THEN
279  xpq(1)=xpq(1)+0.2*xpv
280  xpq(-1)=xpq(-1)+0.2*xpv
281  xpq(2)=xpq(2)+0.8*xpv
282  xpq(-2)=xpq(-2)+0.8*xpv
283  ELSEIF(kfa.EQ.22.AND.mint(105).EQ.333) THEN
284  xpq(3)=xpq(3)+xpv
285  xpq(-3)=xpq(-3)+xpv
286  ELSEIF(kfa.EQ.22.AND.mint(105).EQ.443) THEN
287  xpq(4)=xpq(4)+xpv
288  xpq(-4)=xpq(-4)+xpv
289  IF(mstp(55).GE.9) THEN
290  DO 230 kfl=-6,6
291  xpq(kfl)=0.
292  230 CONTINUE
293  ENDIF
294  ELSE
295  xpq(1)=xpq(1)+0.5*xpv
296  xpq(-1)=xpq(-1)+0.5*xpv
297  xpq(2)=xpq(2)+0.5*xpv
298  xpq(-2)=xpq(-2)+0.5*xpv
299  ENDIF
300 
301 C...Rescale for gammaVDM by effective gamma -> rho coupling.
302  IF(kfa.EQ.22.AND.mint(109).EQ.2) THEN
303  DO 240 kfl=-6,6
304  xpq(kfl)=vint(281)*xpq(kfl)
305  240 CONTINUE
306  vint(232)=vint(281)*xpv
307  ENDIF
308 
309 C...Isospin conjugation for neutron.
310  ELSEIF(kfa.EQ.2112) THEN
311  xps=xpq(1)
312  xpq(1)=xpq(2)
313  xpq(2)=xps
314  xps=xpq(-1)
315  xpq(-1)=xpq(-2)
316  xpq(-2)=xps
317 
318 C...Simple recipes for hyperon (average valence structure function).
319  ELSEIF(kfa.EQ.3122.OR.kfa.EQ.3112.OR.kfa.EQ.3212.OR.kfa.EQ.3222
320  &.OR.kfa.EQ.3312.OR.kfa.EQ.3322.OR.kfa.EQ.3334) THEN
321  xpval=(xpq(1)+xpq(2)-xpq(-1)-xpq(-2))/3.
322  xpsea=0.5*(xpq(-1)+xpq(-2))
323  xpq(1)=xpsea
324  xpq(2)=xpsea
325  xpq(-1)=xpsea
326  xpq(-2)=xpsea
327  xpq(kfa/1000)=xpq(kfa/1000)+xpval
328  xpq(mod(kfa/100,10))=xpq(mod(kfa/100,10))+xpval
329  xpq(mod(kfa/10,10))=xpq(mod(kfa/10,10))+xpval
330  ENDIF
331 
332 C...Charge conjugation for antiparticle.
333  IF(kf.LT.0) THEN
334  DO 250 kfl=1,25
335  IF(kfl.EQ.21.OR.kfl.EQ.22.OR.kfl.EQ.23.OR.kfl.EQ.25) goto 250
336  xps=xpq(kfl)
337  xpq(kfl)=xpq(-kfl)
338  xpq(-kfl)=xps
339  250 CONTINUE
340  ENDIF
341 
342 C...Allow gluon also in position 21.
343  xpq(21)=xpq(0)
344 
345 C...Check positivity and reset above maximum allowed flavour.
346  DO 260 kfl=-25,25
347  xpq(kfl)=max(0.,xpq(kfl))
348  IF(iabs(kfl).GT.mstp(58).AND.iabs(kfl).LE.8) xpq(kfl)=0.
349  260 CONTINUE
350 
351 C...Formats for error printouts.
352  5000 FORMAT(' Error: x value outside physical range; x =',1p,e12.3)
353  5100 FORMAT(' Error: illegal particle code for structure function;',
354  &' KF =',i5)
355  5200 FORMAT(' Error: unknown structure function; KF, library, set =',
356  &3i5)
357 
358  RETURN
359  END