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