EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pygvmd.F
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pygvmd.F
1  SUBROUTINE pygvmd(ISET,KF,X,Q2,P2,ALAM,XPGA)
2 C...Purpose: to evaluate the VMD parton distributions of a photon,
3 C...evolved homogeneously from an initial scale P2 to Q2.
4 C...Does not include dipole suppression factor.
5 C...ISET is parton distribution set, see above;
6 C...additionally ISET=0 is used for the evolution of an anomalous photon
7 C...which branched at a scale P2 and then evolved homogeneously to Q2.
8 C...ALAM is the 4-flavour Lambda, which is automatically converted
9 C...to 3- and 5-flavour equivalents as needed.
10  dimension xpga(-6:6)
11  DATA pmc/1.3/, pmb/4.6/, aem/0.007297/, aem2pi/0.0011614/
12 
13 C...Reset output.
14  DO 100 kfl=-6,6
15  xpga(kfl)=0.
16  100 CONTINUE
17  kfa=iabs(kf)
18 
19 C...Calculate Lambda; protect against unphysical Q2 and P2 input.
20  alam3=alam*(pmc/alam)**(2./27.)
21  alam5=alam*(alam/pmb)**(2./23.)
22  p2eff=max(p2,1.2*alam3**2)
23  IF(kfa.EQ.4) p2eff=max(p2eff,pmc**2)
24  IF(kfa.EQ.5) p2eff=max(p2eff,pmb**2)
25  q2eff=max(q2,p2eff)
26 
27 C...Find number of flavours at lower and upper scale.
28  nfp=4
29  IF(p2eff.LT.pmc**2) nfp=3
30  IF(p2eff.GT.pmb**2) nfp=5
31  nfq=4
32  IF(q2eff.LT.pmc**2) nfq=3
33  IF(q2eff.GT.pmb**2) nfq=5
34 
35 C...Find s as sum of 3-, 4- and 5-flavour parts.
36  s=0.
37  IF(nfp.EQ.3) THEN
38  q2div=pmc**2
39  IF(nfq.EQ.3) q2div=q2eff
40  s=s+(6./27.)*log(log(q2div/alam3**2)/log(p2eff/alam3**2))
41  ENDIF
42  IF(nfp.LE.4.AND.nfq.GE.4) THEN
43  p2div=p2eff
44  IF(nfp.EQ.3) p2div=pmc**2
45  q2div=q2eff
46  IF(nfq.EQ.5) q2div=pmb**2
47  s=s+(6./25.)*log(log(q2div/alam**2)/log(p2div/alam**2))
48  ENDIF
49  IF(nfq.EQ.5) THEN
50  p2div=pmb**2
51  IF(nfp.EQ.5) p2div=p2eff
52  s=s+(6./23.)*log(log(q2eff/alam5**2)/log(p2div/alam5**2))
53  ENDIF
54 
55 C...Calculate frequent combinations of x and s.
56  x1=1.-x
57  xl=-log(x)
58  s2=s**2
59  s3=s**3
60  s4=s**4
61 
62 C...Evaluate homogeneous anomalous parton distributions below or
63 C...above threshold.
64  IF(iset.EQ.0) THEN
65  IF(q2.LE.p2.OR.(kfa.EQ.4.AND.q2.LT.pmc**2).OR.
66  &(kfa.EQ.5.AND.q2.LT.pmb**2)) THEN
67  xval = x * 1.5 * (x**2+x1**2)
68  xglu = 0.
69  xsea = 0.
70  ELSE
71  xval = (1.5/(1.-0.197*s+4.33*s2)*x**2 + (1.5+2.10*s)/
72  & (1.+3.29*s)*x1**2 + 5.23*s/(1.+1.17*s+19.9*s3)*x*x1) *
73  & x**(1./(1.+1.5*s)) * (1.-x**2)**(2.667*s)
74  xglu = 4.*s/(1.+4.76*s+15.2*s2+29.3*s4) *
75  & x**(-2.03*s/(1.+2.44*s)) * (x1*xl)**(1.333*s) *
76  & ((4.*x**2+7.*x+4.)*x1/3. - 2.*x*(1.+x)*xl)
77  xsea = s2/(1.+4.54*s+8.19*s2+8.05*s3) *
78  & x**(-1.54*s/(1.+1.29*s)) * x1**(2.667*s) *
79  & ((8.-73.*x+62.*x**2)*x1/9. + (3.-8.*x**2/3.)*x*xl +
80  & (2.*x-1.)*x*xl**2)
81  ENDIF
82 
83 C...Evaluate set 1D parton distributions below or above threshold.
84  ELSEIF(iset.EQ.1) THEN
85  IF(q2.LE.p2.OR.(kfa.EQ.4.AND.q2.LT.pmc**2).OR.
86  &(kfa.EQ.5.AND.q2.LT.pmb**2)) THEN
87  xval = 1.294 * x**0.80 * x1**0.76
88  xglu = 1.273 * x**0.40 * x1**1.76
89  xsea = 0.100 * x1**3.76
90  ELSE
91  xval = 1.294/(1.+0.252*s+3.079*s2) * x**(0.80-0.13*s) *
92  & x1**(0.76+0.667*s) * xl**(2.*s)
93  xglu = 7.90*s/(1.+5.50*s) * exp(-5.16*s) *
94  & x**(-1.90*s/(1.+3.60*s)) * x1**1.30 * xl**(0.50+3.*s) +
95  & 1.273 * exp(-10.*s) * x**0.40 * x1**(1.76+3.*s)
96  xsea = (0.1-0.397*s2+1.121*s3)/(1.+5.61*s2+5.26*s3) *
97  & x**(-7.32*s2/(1.+10.3*s2)) *
98  & x1**((3.76+15.*s+12.*s2)/(1.+4.*s))
99  xsea0 = 0.100 * x1**3.76
100  ENDIF
101 
102 C...Evaluate set 1M parton distributions below or above threshold.
103  ELSEIF(iset.EQ.2) THEN
104  IF(q2.LE.p2.OR.(kfa.EQ.4.AND.q2.LT.pmc**2).OR.
105  &(kfa.EQ.5.AND.q2.LT.pmb**2)) THEN
106  xval = 0.8477 * x**0.51 * x1**1.37
107  xglu = 3.42 * x**0.255 * x1**2.37
108  xsea = 0.
109  ELSE
110  xval = 0.8477/(1.+1.37*s+2.18*s2+3.73*s3) * x**(0.51+0.21*s)
111  & * x1**1.37 * xl**(2.667*s)
112  xglu = 24.*s/(1.+9.6*s+0.92*s2+14.34*s3) * exp(-5.94*s) *
113  & x**((-0.013-1.80*s)/(1.+3.14*s)) * x1**(2.37+0.4*s) *
114  & xl**(0.32+3.6*s) + 3.42 * exp(-12.*s) * x**0.255 *
115  & x1**(2.37+3.*s)
116  xsea = 0.842*s/(1.+21.3*s-33.2*s2+229.*s3) *
117  & x**((0.13-2.90*s)/(1.+5.44*s)) * x1**(3.45+0.5*s) *
118  & xl**(2.8*s)
119  xsea0 = 0.
120  ENDIF
121 
122 C...Evaluate set 2D parton distributions below or above threshold.
123  ELSEIF(iset.EQ.3) THEN
124  IF(q2.LE.p2.OR.(kfa.EQ.4.AND.q2.LT.pmc**2).OR.
125  &(kfa.EQ.5.AND.q2.LT.pmb**2)) THEN
126  xval = x**0.46 * x1**0.64 + 0.76 * x
127  xglu = 1.925 * x1**2
128  xsea = 0.242 * x1**4
129  ELSE
130  xval = (1.+0.186*s)/(1.-0.209*s+1.495*s2) * x**(0.46+0.25*s)
131  & * x1**((0.64+0.14*s+5.*s2)/(1.+s)) * xl**(1.9*s) +
132  & (0.76+0.4*s) * x * x1**(2.667*s)
133  xglu = (1.925+5.55*s+147.*s2)/(1.-3.59*s+3.32*s2) *
134  & exp(-18.67*s) * x**((-5.81*s-5.34*s2)/(1.+29.*s-4.26*s2))
135  & * x1**((2.-5.9*s)/(1.+1.7*s)) * xl**(9.3*s/(1.+1.7*s))
136  xsea = (0.242-0.252*s+1.19*s2)/(1.-0.607*s+21.95*s2) *
137  & x**(-12.1*s2/(1.+2.62*s+16.7*s2)) * x1**4 * xl**s
138  xsea0 = 0.242 * x1**4
139  ENDIF
140 
141 C...Evaluate set 2M parton distributions below or above threshold.
142  ELSEIF(iset.EQ.4) THEN
143  IF(q2.LE.p2.OR.(kfa.EQ.4.AND.q2.LT.pmc**2).OR.
144  &(kfa.EQ.5.AND.q2.LT.pmb**2)) THEN
145  xval = 1.168 * x**0.50 * x1**2.60 + 0.965 * x
146  xglu = 1.808 * x1**2
147  xsea = 0.209 * x1**4
148  ELSE
149  xval = (1.168+1.771*s+29.35*s2) * exp(-5.776*s) *
150  & x**((0.5+0.208*s)/(1.-0.794*s+1.516*s2)) *
151  & x1**((2.6+7.6*s)/(1.+5.*s)) * xl**(5.15*s/(1.+2.*s)) +
152  & (0.965+22.35*s)/(1.+18.4*s) * x * x1**(2.667*s)
153  xglu = (1.808+29.9*s)/(1.+26.4*s) * exp(-5.28*s) *
154  & x**((-5.35*s-10.11*s2)/(1.+31.71*s)) *
155  & x1**((2.-7.3*s+4.*s2)/(1.+2.5*s)) *
156  & xl**(10.9*s/(1.+2.5*s))
157  xsea = (0.209+0.644*s2)/(1.+0.319*s+17.6*s2) *
158  & x**((-0.373*s-7.71*s2)/(1.+0.815*s+11.0*s2)) *
159  & x1**(4.+s) * xl**(0.45*s)
160  xsea0 = 0.209 * x1**4
161  ENDIF
162  ENDIF
163 
164 C...Threshold factors for c and b sea.
165  sll=log(log(q2eff/alam**2)/log(p2eff/alam**2))
166  xchm=0.
167  IF(q2.GT.pmc**2.AND.q2.GT.1.001*p2eff) THEN
168  sch=max(0.,log(log(pmc**2/alam**2)/log(p2eff/alam**2)))
169  IF(iset.EQ.0) THEN
170  xchm=xsea*(1.-(sch/sll)**2)
171  ELSE
172  xchm=max(0.,xsea-xsea0*x1**(2.667*s))*(1.-sch/sll)
173  ENDIF
174  ENDIF
175  xbot=0.
176  IF(q2.GT.pmb**2.AND.q2.GT.1.001*p2eff) THEN
177  sbt=max(0.,log(log(pmb**2/alam**2)/log(p2eff/alam**2)))
178  IF(iset.EQ.0) THEN
179  xbot=xsea*(1.-(sbt/sll)**2)
180  ELSE
181  xbot=max(0.,xsea-xsea0*x1**(2.667*s))*(1.-sbt/sll)
182  ENDIF
183  ENDIF
184 
185 C...Fill parton distributions.
186  xpga(0)=xglu
187  xpga(1)=xsea
188  xpga(2)=xsea
189  xpga(3)=xsea
190  xpga(4)=xchm
191  xpga(5)=xbot
192  xpga(kfa)=xpga(kfa)+xval
193  DO 110 kfl=1,5
194  xpga(-kfl)=xpga(kfl)
195  110 CONTINUE
196 
197  RETURN
198  END