]> git.uio.no Git - u/mrichter/AliRoot.git/blame - LHAPDF/lhapdf5.3.1/wrapwhitg.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix for IsTriggerInputFired,GetFiredTriggerInputs
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapwhitg.f
CommitLineData
4e9e3152 1 subroutine WHITGevolvep(xin,qin,p2in,ip2in,pdf)
2 include 'parmsetup.inc'
3 real*8 xin,qin,q2in,p2in,pdf(-6:6),xval(45),qcdl4,qcdl5
4 real*8 upv,dnv,usea,dsea,str,chm,bot,top,glu
5 character*16 name(nmxset)
6 integer nmem(nmxset),ndef(nmxset),mmem
7 common/NAME/name,nmem,ndef,mmem
8 integer nset
9
10 save
11 call getnset(iset)
12 call getnmem(iset,imem)
13
14 if(imem.eq.1.or.imem.eq.0) then
15 call SFWHI1(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
16
17 elseif(imem.eq.2) then
18 call SFWHI2(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
19
20 elseif(imem.eq.3) then
21 call SFWHI3(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
22
23 elseif(imem.eq.4) then
24 call SFWHI4(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
25
26 elseif(imem.eq.5) then
27 call SFWHI5(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
28
29 elseif(imem.eq.6) then
30 call SFWHI6(xin,qin,upv,dnv,usea,dsea,str,chm,glu)
31
32 else
33 CONTINUE
34 endif
35
36 pdf(-6)= 0.0d0
37 pdf(6)= 0.0d0
38 pdf(-5)= 0.0d0
39 pdf(5 )= 0.0d0
40 pdf(-4)= chm
41 pdf(4 )= chm
42 pdf(-3)= str
43 pdf(3 )= str
44 pdf(-2)= usea
45 pdf(2 )= upv
46 pdf(-1)= dsea
47 pdf(1 )= dnv
48 pdf(0 )= glu
49
50 return
51ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
52 entry WHITGread(nset)
53 read(1,*)nmem(nset),ndef(nset)
54 return
55c
56ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
57 entry WHITGalfa(alfas,qalfa)
58 call getnset(iset)
59 call getnmem(iset,imem)
60 call GetOrderAsM(iset,iord)
61 call Getlam4M(iset,imem,qcdl4)
62 call Getlam5M(iset,imem,qcdl5)
63 call aspdflib(alfas,Qalfa,iord,qcdl5)
64
65 return
66c
67ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
68 entry WHITGinit(Eorder,Q2fit)
69 return
70c
71ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
72 entry WHITGpdf(mem)
73 call getnset(iset)
74 call setnmem(iset,mem)
75c imem = mem
76 return
77c
78 1000 format(5e13.5)
79 end
80ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
81c-------------------------------------------------------
82 subroutine SFWHI1(ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL)
83c-------------------------------------------------------
84c WHIT1 parton distribution in the photon
85c
86c INPUT: integer ic : if ic=0 then qc=0
87c else qc is calculated
88c DOUBLE PRECISION Q2 : energy scale Q^2 (GeV^2)
89c DOUBLE PRECISION x : energy fraction
90c
91c OUTPUT: DOUBLE PRECISION qu : up-quark dist.
92c DOUBLE PRECISION qd : down- or strange-quark dist.
93c DOUBLE PRECISION qc : charm-quark dist.
94c DOUBLE PRECISION g : gluon dist.
95c-------------------------------------------------------
96c Modified by M.Tanaka on July 22, 1994.
97c The bug pointed out by M.Drees is fixed.
98c-------------------------------------------------------
99c Modified by I.Watanabe on July 22, 1994.
100c-------------------------------------------------------
101 implicit none
102 external WHIT1G
103 double precision
104 + ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL
105c arg
106 integer ic
107 DOUBLE PRECISION Q2,x
108 DOUBLE PRECISION qu,qd,qc,g
109c const
110 DOUBLE PRECISION q42it,q52it,lam42,lam52
111 DOUBLE PRECISION alinv,mc,PI
112c local
113 DOUBLE PRECISION qv,qsea,cv,cs,dcv,dcs
114 DOUBLE PRECISION A0val,A1val,A2val,Bval,Cval,
115 $ A0sea,B0sea,BB0sea,C0sea
116 DOUBLE PRECISION A0dcv,A1dcv,A2dcv,A3dcv,Bdcv,Cdcv
117 DOUBLE PRECISION Adcs, B0dcs, B1dcs, Cdcs
118 DOUBLE PRECISION x1,x2,mc2q2
119 DOUBLE PRECISION s,s2,s3,s4,prsccf,alstpi
120 DOUBLE PRECISION WHIT1G
121c parameters
122 parameter(lam42=0.16d0, lam52=0.091411319d0)
123 parameter(Q42IT=4.0d0, Q52IT=100.0d0)
124 parameter(alinv=137.036d0, mc=1.5d0)
125 parameter(pi=3.14159265358979323846d0)
126 common /scale/ s,s2,s3,s4,prsccf
127c
128c begin
129 x=ZX
130 Q2=ZQ*ZQ
131 ic=1
132c
133 x1=1.0d0-x
134 x2=x**2
135 mc2q2=mc**2/Q2
136c
137 if(Q2.lt.100.0d0) then
138c under 100 GeV^2
139c
140c set scale s
141 if(Q2.lt.4.0d0) then
142cccc for under 4GeV^2 prescription
143 s= 0.0d0
144 prsccf = log(Q2/LAM42)/ log(Q42IT/LAM42)
145 alstpi = 6.0d0/25.0d0/ log(Q42IT/LAM42)
146 else
147 s= log( log(Q2/LAM42)/ log(Q42IT/LAM42))
148 prsccf = 1.0d0
149 alstpi = 6.0d0/25.0d0/ log(Q2/LAM42)
150 endif
151 s2=s**2
152 s3=s2*s
153 s4=s2**2
154c
155cccccc WHIT1 quark (U100)
156c
157 A0val= 1.882000d+00+s*( 1.213000d+00)+s2*( 6.970000d-01)
158 A1val= s*(-2.361000d+00)+s2*(-1.136000d+00)
159 A2val= s*( 5.280000d-01)+s2*( 2.406000d+00)
160 Bval = 5.000000d-01+s*( 2.107000d-02)+s2*( 4.130000d-03)
161 Cval = 2.500000d-01+s*(-2.376000d-01)+s2*( 2.018000d-01)
162 $ +s3*(-5.040000d-02)
163 A0sea= 6.510000d-01+s*( 1.291000d+00)+s2*(-4.470000d+00)
164 $ +s3*( 5.140000d+00)+s4*(-2.091000d+00)
165 B0sea=-3.820000d-02+s*( 9.010000d-02)+s2*(-1.356000d+00)
166 $ +s3*( 1.582000d+00)+s4*(-6.440000d-01)
167 BB0sea=2.084000d+00+s*( 7.740000d+00)+s2*(-2.970000d+01)
168 $ +s3*( 3.860000d+01)+s4*(-1.705000d+01)
169 C0sea= 7.000000d+00+s*(-1.608000d+01)+s2*( 4.670000d+01)
170 $ +s3*(-5.710000d+01)+s4*( 2.386000d+01)
171c
172 qv = prsccf/alinv/x*
173 $ (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
174 qsea= prsccf/alinv/x*
175 $ A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
176c
177 qu = qv/3.0d0 + qsea/6.0d0
178 qu = qu*x
179 ZUV=qu
180 ZUB=qu
181 qd = qv/12.0d0 + qsea/6.0d0
182 qd = qd*x
183 ZDV=qd
184 ZDB=qd
185 ZSB=qd
186c
187 if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
188 call WHIT1Q(x,mc2q2,cv,cs)
189 qc = cv/alinv/2.0d0/PI + cs*alstpi
190 qc = qc*x
191 ZCB=qc
192 else
193 qc = 0.0d0
194 ZCB=qc
195 endif
196c
197 g = WHIT1G(x,Q2)
198 g = g*x
199 ZGL=g
200c
201 else
202c over 100 GeV^2
203c
204c set scale s
205 s= log( log(Q2/LAM52)/ log(Q52IT/LAM52))
206 prsccf = 1.0d0
207 alstpi = 6.0d0/23.0d0/ log(Q2/LAM52)
208 s2=s**2
209 s3=s2*s
210 s4=s2**2
211c
212cccccc WHIT1 quark (O100)
213c
214 A0val= 3.058000d+00+s*( 2.474000d+00)+s2*( 1.002000d+00)
215 A1val=-2.182000d+00+s*(-4.480000d+00)+s2*(-2.251000d-01)
216 A2val= 1.522000d+00+s*( 4.310000d+00)+s2*( 1.314000d+00)
217 Bval = 5.170000d-01+s*( 4.040000d-02)+s2*(-2.100000d-02)
218 Cval = 1.655000d-01+s*(-2.062000d-02)+s2*( 5.360000d-02)
219 A0sea= 6.250000d-01+s*(-5.890000d-01)+s2*( 4.180000d+00)
220 $ +s3*(-1.206000d+01)+s4*( 1.257000d+01)
221 B0sea=-2.492000d-01+s*(-4.110000d-01)+s2*( 9.660000d-01)
222 $ +s3*(-2.584000d+00)+s4*( 2.670000d+00)
223 BB0sea=2.100000d+00+s*(-5.750000d+00)+s2*( 4.780000d+01)
224 $ +s3*(-1.407000d+02)+s4*( 1.476000d+02)
225 C0sea= 4.780000d+00+s*( 4.860000d+00)+s2*(-4.890000d+01)
226 $ +s3*( 1.477000d+02)+s4*(-1.602000d+02)
227c
228 qv = 1.0d0/alinv/x*
229 $ (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
230 qsea= 1.0d0/alinv/x*
231 $ A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
232c
233 qu = qv/3.0d0 + qsea/6.0d0
234 qu = qu*x
235 ZUV=qu
236 ZUB=qu
237 qd = qv/12.0d0 + qsea/6.0d0
238 qd = qd*x
239 ZDV=qd
240 ZDB=qd
241 ZSB=qd
242 g = WHIT1G(x,Q2)
243 g = g*x
244 ZGL=g
245c
246 if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
247 A0dcv= s*( 1.219000d-01)+s2*( 6.200000d+00)
248 $ +s3*(-2.504000d+01)+s4*( 3.098000d+01)
249 A1dcv= s*( 1.913000d+00)+s2*(-7.690000d+01)
250 $ +s3*( 3.180000d+02)+s4*(-3.920000d+02)
251 A2dcv= s*(-7.160000d+00)+s2*( 2.503000d+02)
252 $ +s3*(-1.062000d+03)+s4*( 1.308000d+03)
253 A3dcv= s*( 3.190000d+00)+s2*(-2.301000d+02)
254 $ +s3*( 1.012000d+03)+s4*(-1.250000d+03)
255 Bdcv = 4.990000d-01+s*( 3.470000d+00)+s2*(-1.526000d+01)
256 $ +s3*( 1.967000d+01)
257 Cdcv = 3.290000d-01+s*( 8.240000d+00)+s2*(-3.800000d+01)
258 $ +s3*( 4.630000d+01)
259 Adcs = s*(-1.815000d-02)+s2*( 2.043000d-03)
260 $ +s3*(-4.130000d-03)
261 B0dcs=-3.086000d-01+s*(-2.565000d-01)+s2*( 9.840000d-02)
262 B1dcs= 1.376000d+00+s*(-4.630000d-01)+s2*( 1.232000d+00)
263 Cdcs = 3.650000d+00+s*( 7.290000d-01)+s2*(-7.570000d+00)
264 $ +s3*( 7.790000d+00)
265c
266 dcv = 1.0d0/alinv/x*
267 $ (A0dcv+x*A1dcv+x2*A2dcv+x2*x*A3dcv) * x**Bdcv * x1**Cdcv
268 dcs = 1.0d0/alinv/x*
269 $ Adcs * x**(B0dcs+B1dcs*x) * x1**Cdcs
270c
271 call WHIT1Q(x,mc*mc/Q2,cv,cs)
272 qc = cv/alinv/2.0d0/PI + cs*alstpi + dcs + dcv
273 qc = qc*x
274 ZCB=qc
275 else
276 qc = 0.0d0
277 ZCB=qc
278 endif
279 endif
280c
281 return
282 end
283c-------------------------------------------------------
284 subroutine SFWHI2(ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL)
285c-------------------------------------------------------
286c WHIT2 parton distribution in the photon
287c
288c INPUT: integer ic : if ic=0 then qc=0
289c else qc is calculated
290c DOUBLE PRECISION Q2 : energy scale Q^2 (GeV^2)
291c DOUBLE PRECISION x : energy fraction
292c
293c OUTPUT: DOUBLE PRECISION qu : up-quark dist.
294c DOUBLE PRECISION qd : down- or strange-quark dist.
295c DOUBLE PRECISION qc : charm-quark dist.
296c DOUBLE PRECISION g : gluon dist.
297c-------------------------------------------------------
298c Modified by M.Tanaka on July 22, 1994.
299c The bug pointed out by M.Drees is fixed.
300c-------------------------------------------------------
301c Modified by I.Watanabe on July 22, 1994.
302c-------------------------------------------------------
303 implicit none
304 external WHIT2G
305 double precision
306 + ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL
307c arg
308 integer ic
309 DOUBLE PRECISION Q2,x
310 DOUBLE PRECISION qu,qd,qc,g
311c const
312 DOUBLE PRECISION q42it,q52it,lam42,lam52
313 DOUBLE PRECISION alinv,mc,PI
314c local
315 DOUBLE PRECISION qv,qsea,cv,cs,dcv,dcs
316 DOUBLE PRECISION A0val,A1val,A2val,Bval,Cval,
317 $ A0sea,B0sea,BB0sea,C0sea
318 DOUBLE PRECISION A0dcv,A1dcv,A2dcv,A3dcv,Bdcv,Cdcv
319 DOUBLE PRECISION Adcs, B0dcs, B1dcs, Cdcs
320 DOUBLE PRECISION x1,x2,mc2q2
321 DOUBLE PRECISION s,s2,s3,s4,prsccf,alstpi
322 DOUBLE PRECISION WHIT2G
323c parameters
324 parameter(lam42=0.16d0, lam52=0.091411319d0)
325 parameter(Q42IT=4.0d0, Q52IT=100.0d0)
326 parameter(alinv=137.036d0, mc=1.5d0)
327 parameter(pi=3.14159265358979323846d0)
328 common /scale/ s,s2,s3,s4,prsccf
329c
330c begin
331 x=ZX
332 Q2=ZQ*ZQ
333 ic=1
334c
335 x1=1.0d0-x
336 x2=x**2
337 mc2q2=mc**2/Q2
338c
339 if(Q2.lt.100.0d0) then
340c under 100 GeV^2
341c
342c set scale s
343 if(Q2.lt.4.0d0) then
344cccc for under 4GeV^2 prescription
345 s= 0.0d0
346 prsccf = log(Q2/LAM42)/ log(Q42IT/LAM42)
347 alstpi = 6.0d0/25.0d0/ log(Q42IT/LAM42)
348 else
349 s= log( log(Q2/LAM42)/ log(Q42IT/LAM42))
350 prsccf = 1.0d0
351 alstpi = 6.0d0/25.0d0/ log(Q2/LAM42)
352 endif
353 s2=s**2
354 s3=s2*s
355 s4=s2**2
356c
357cccccc WHIT2 quark (U100)
358c
359 A0val= 1.882000d+00+s*( 1.213000d+00)+s2*( 6.970000d-01)
360 A1val= s*(-2.361000d+00)+s2*(-1.136000d+00)
361 A2val= s*( 5.280000d-01)+s2*( 2.406000d+00)
362 Bval= 5.000000d-01+s*( 2.107000d-02)+s2*( 4.130000d-03)
363 Cval= 2.500000d-01+s*(-2.376000d-01)+s2*( 2.018000d-01)
364 $ +s3*(-5.040000d-02)
365 A0sea= 1.237000d+00+s*( 3.390000d+00)+s2*(-1.075000d+01)
366 $ +s3*( 1.246000d+01)+s4*(-5.580000d+00)
367 B0sea=-7.270000d-02+s*( 1.748000d-01)+s2*(-1.392000d+00)
368 $ +s3*( 1.711000d+00)+s4*(-7.960000d-01)
369 BB0sea=4.290000d+00+s*( 1.787000d+01)+s2*(-5.810000d+01)
370 $ +s3*( 8.190000d+01)+s4*(-4.140000d+01)
371 C0sea= 1.434000d+01+s*(-4.490000d+01)+s2*( 1.197000d+02)
372 $ +s3*(-1.585000d+02)+s4*( 7.530000d+01)
373c
374 qv = prsccf/alinv/x*
375 $ (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
376 qsea= prsccf/alinv/x*
377 $ A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
378c
379 qu = qv/3.0d0 + qsea/6.0d0
380 qu = qu*x
381 ZUV=qu
382 ZUB=qu
383 qd = qv/12.0d0 + qsea/6.0d0
384 qd = qd*x
385 ZDV=qd
386 ZDB=qd
387 ZSB=qd
388c
389 if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
390 call WHIT2Q(x,mc2q2,cv,cs)
391 qc = cv/alinv/2.0d0/PI + cs*alstpi
392 qc = qc*x
393 ZCB=qc
394 else
395 qc = 0.0d0
396 ZCB=qc
397 endif
398c
399 g = WHIT2G(x,Q2)
400 g = g*x
401 ZGL=g
402c
403 else
404c over 100 GeV^2
405c
406c set scale s
407 s= log( log(Q2/LAM52)/ log(Q52IT/LAM52))
408 prsccf = 1.0d0
409 alstpi = 6.0d0/23.0d0/ log(Q2/LAM52)
410 s2=s**2
411 s3=s2*s
412 s4=s2**2
413c
414cccccc WHIT2 quark (O100)
415c
416 A0val= 3.058000d+00+s*( 2.474000d+00)+s2*( 1.002000d+00)
417 A1val=-2.182000d+00+s*(-4.480000d+00)+s2*(-2.259000d-01)
418 A2val= 1.522000d+00+s*( 4.300000d+00)+s2*( 1.315000d+00)
419 Bval = 5.170000d-01+s*( 4.030000d-02)+s2*(-2.098000d-02)
420 Cval = 1.655000d-01+s*(-2.063000d-02)+s2*( 5.370000d-02)
421 A0sea= 1.287000d+00+s*(-2.069000d+00)+s2*( 1.157000d+01)
422 $ +s3*(-3.570000d+01)+s4*( 3.740000d+01)
423 B0sea=-2.340000d-01+s*(-4.430000d-01)+s2*( 1.235000d+00)
424 $ +s3*(-3.720000d+00)+s4*( 3.840000d+00)
425 BB0sea=6.460000d+00+s*(-1.048000d+01)+s2*( 8.980000d+01)
426 $ +s3*(-2.847000d+02)+s4*( 2.998000d+02)
427 C0sea= 5.350000d+00+s*( 1.011000d+01)+s2*(-1.337000d+02)
428 $ +s3*( 4.270000d+02)+s4*(-4.570000d+02)
429c
430 qv = 1.0d0/alinv/x*
431 $ (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
432 qsea= 1.0d0/alinv/x*
433 $ A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
434c
435 qu = qv/3.0d0 + qsea/6.0d0
436 qu = qu*x
437 ZUV=qu
438 ZUB=qu
439 qd = qv/12.0d0 + qsea/6.0d0
440 qd = qd*x
441 ZDV=qd
442 ZDB=qd
443 ZSB=qd
444 g = WHIT2G(x,Q2)
445 g = g*x
446 ZGL=g
447c
448 if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
449 A0dcv= s*( 1.219000d-01)+s2*( 6.200000d+00)
450 $ +s3*(-2.504000d+01)+s4*( 3.098000d+01)
451 A1dcv= s*( 1.913000d+00)+s2*(-7.690000d+01)
452 $ +s3*( 3.180000d+02)+s4*(-3.920000d+02)
453 A2dcv= s*(-7.160000d+00)+s2*( 2.503000d+02)
454 $ +s3*(-1.062000d+03)+s4*( 1.308000d+03)
455 A3dcv= s*( 3.190000d+00)+s2*(-2.301000d+02)
456 $ +s3*( 1.012000d+03)+s4*(-1.250000d+03)
457 Bdcv = 4.990000d-01+s*( 3.470000d+00)+s2*(-1.526000d+01)
458 $ +s3*( 1.967000d+01)
459 Cdcv = 3.290000d-01+s*( 8.240000d+00)+s2*(-3.800000d+01)
460 $ +s3*( 4.630000d+01)
461 Adcs = s*(-2.786000d-02)+s2*( 3.490000d-02)
462 $ +s3*(-2.223000d-02)
463 B0dcs=-3.141000d-01+s*(-4.250000d-01)+s2*( 1.564000d-01)
464 B1dcs= 4.720000d+00+s*(-5.480000d+00)+s2*( 2.686000d+00)
465 Cdcs = 2.961000d+00+s*( 7.760000d-01)+s2*(-8.280000d+00)
466 $ +s3*( 9.780000d+00)
467c
468 dcv = 1.0d0/alinv/x*
469 $ (A0dcv+x*A1dcv+x2*A2dcv+x2*x*A3dcv) * x**Bdcv * x1**Cdcv
470 dcs = 1.0d0/alinv/x*
471 $ Adcs * x**(B0dcs+B1dcs*x) * x1**Cdcs
472c
473 call WHIT2Q(x,mc*mc/Q2,cv,cs)
474 qc = cv/alinv/2.0d0/PI + cs*alstpi + dcs + dcv
475 qc = qc*x
476 ZCB=qc
477 else
478 qc = 0.0d0
479 ZCB=qc
480 endif
481 endif
482c
483 return
484 end
485c-------------------------------------------------------
486 subroutine SFWHI3(ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL)
487c-------------------------------------------------------
488c WHIT3 parton distribution in the photon
489c
490c INPUT: integer ic : if ic=0 then qc=0
491c else qc is calculated
492c DOUBLE PRECISION Q2 : energy scale Q^2 (GeV^2)
493c DOUBLE PRECISION x : energy fraction
494c
495c OUTPUT: DOUBLE PRECISION qu : up-quark dist.
496c DOUBLE PRECISION qd : down- or strange-quark dist.
497c DOUBLE PRECISION qc : charm-quark dist.
498c DOUBLE PRECISION g : gluon dist.
499c-------------------------------------------------------
500c Modified by M.Tanaka on July 22, 1994.
501c The bug pointed out by M.Drees is fixed.
502c-------------------------------------------------------
503c Modified by I.Watanabe on July 22, 1994.
504c-------------------------------------------------------
505 implicit none
506 external whit3g
507 double precision
508 + ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL
509c arg
510 integer ic
511 DOUBLE PRECISION Q2,x
512 DOUBLE PRECISION qu,qd,qc,g
513c const
514 DOUBLE PRECISION q42it,q52it,lam42,lam52
515 DOUBLE PRECISION alinv,mc,PI
516c local
517 DOUBLE PRECISION qv,qsea,cv,cs,dcv,dcs
518 DOUBLE PRECISION A0val,A1val,A2val,Bval,Cval,
519 $ A0sea,B0sea,BB0sea,C0sea
520 DOUBLE PRECISION A0dcv,A1dcv,A2dcv,A3dcv,Bdcv,Cdcv
521 DOUBLE PRECISION Adcs, B0dcs, B1dcs, Cdcs
522 DOUBLE PRECISION x1,x2,mc2q2
523 DOUBLE PRECISION s,s2,s3,s4,prsccf,alstpi
524 DOUBLE PRECISION WHIT3g
525c parameters
526 parameter(lam42=0.16d0, lam52=0.091411319d0)
527 parameter(Q42IT=4.0d0, Q52IT=100.0d0)
528 parameter(alinv=137.036d0, mc=1.5d0)
529 parameter(pi=3.14159265358979323846d0)
530 common /scale/ s,s2,s3,s4,prsccf
531c
532c begin
533 x=ZX
534 Q2=ZQ*ZQ
535 ic=1
536c
537 x1=1.0d0-x
538 x2=x**2
539 mc2q2=mc**2/Q2
540c
541 if(Q2.lt.100.0d0) then
542c under 100 GeV^2
543c
544c set scale s
545 if(Q2.lt.4.0d0) then
546cccc for under 4GeV^2 prescription
547 s= 0.0d0
548 prsccf = log(Q2/LAM42)/ log(Q42IT/LAM42)
549 alstpi = 6.0d0/25.0d0/ log(Q42IT/LAM42)
550 else
551 s= log( log(Q2/LAM42)/ log(Q42IT/LAM42))
552 prsccf = 1.0d0
553 alstpi = 6.0d0/25.0d0/ log(Q2/LAM42)
554 endif
555 s2=s**2
556 s3=s2*s
557 s4=s2**2
558c
559cccccc WHIT3 quark (U100)
560c
561 A0val= 1.882000d+00+s*( 1.213000d+00)+s2*( 6.970000d-01)
562 A1val= s*(-2.361000d+00)+s2*(-1.136000d+00)
563 A2val= s*( 5.280000d-01)+s2*( 2.406000d+00)
564 Bval = 5.000000d-01+s*( 2.107000d-02)+s2*( 4.130000d-03)
565 Cval = 2.500000d-01+s*(-2.376000d-01)+s2*( 2.018000d-01)
566 $ +s3*(-5.040000d-02)
567 A0sea= 1.587000d+00+s*( 5.050000d+00)+s2*(-1.126000d+01)
568 $ +s3*( 7.560000d+00)+s4*(-1.471000d+00)
569 B0sea=-1.006000d-01+s*( 2.259000d-01)+s2*(-1.195000d+00)
570 $ +s3*( 1.175000d+00)+s4*(-4.460000d-01)
571 BB0sea=5.730000d+00+s*( 2.564000d+01)+s2*(-5.870000d+01)
572 $ +s3*( 6.320000d+01)+s4*(-2.577000d+01)
573 C0sea= 2.136000d+01+s*(-7.290000d+01)+s2*( 1.532000d+02)
574 $ +s3*(-1.679000d+02)+s4*( 6.740000d+01)
575c
576 qv = prsccf/alinv/x*
577 $ (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
578 qsea= prsccf/alinv/x*
579 $ A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
580c
581 qu = qv/3.0d0 + qsea/6.0d0
582 qu = qu*x
583 ZUV=qu
584 ZUB=qu
585 qd = qv/12.0d0 + qsea/6.0d0
586 qd = qd*x
587 ZDV=qd
588 ZDB=qd
589 ZSB=qd
590c
591 if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
592 call WHIT3Q(x,mc2q2,cv,cs)
593 qc = cv/alinv/2.0d0/PI + cs*alstpi
594 qc = qc*x
595 ZCB=qc
596 else
597 qc = 0.0d0
598 ZCB=qc
599 endif
600c
601 g = WHIT3G(x,Q2)
602 g = g*x
603 ZGL=g
604c
605 else
606c over 100 GeV^2
607c
608c set scale s
609 s= log( log(Q2/LAM52)/ log(Q52IT/LAM52))
610 prsccf = 1.0d0
611 alstpi = 6.0d0/23.0d0/ log(Q2/LAM52)
612 s2=s**2
613 s3=s2*s
614 s4=s2**2
615c
616cccccc WHIT3 quark (O100)
617c
618 A0val= 3.058000d+00+s*( 2.474000d+00)+s2*( 1.002000d+00)
619 A1val=-2.182000d+00+s*(-4.480000d+00)+s2*(-2.264000d-01)
620 A2val= 1.522000d+00+s*( 4.300000d+00)+s2*( 1.315000d+00)
621 Bval = 5.170000d-01+s*( 4.030000d-02)+s2*(-2.097000d-02)
622 Cval = 1.655000d-01+s*(-2.064000d-02)+s2*( 5.370000d-02)
623 A0sea= 1.850000d+00+s*(-3.670000d+00)+s2*( 2.714000d+01)
624 $ +s3*(-1.066000d+02)+s4*( 1.309000d+02)
625 B0sea=-2.299000d-01+s*(-4.970000d-01)+s2*( 2.464000d+00)
626 $ +s3*(-9.950000d+00)+s4*( 1.232000d+01)
627 BB0sea=1.042000d+01+s*(-1.074000d+01)+s2*( 1.327000d+02)
628 $ +s3*(-5.390000d+02)+s4*( 6.560000d+02)
629 C0sea= 4.070000d+00+s*( 4.110000d+00)+s2*(-1.719000d+02)
630 $ +s3*( 7.070000d+02)+s4*(-8.590000d+02)
631c
632 qv = 1.0d0/alinv/x*
633 $ (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
634 qsea= 1.0d0/alinv/x*
635 $ A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
636c
637 qu = qv/3.0d0 + qsea/6.0d0
638 qu = qu*x
639 ZUV=qu
640 ZUB=qu
641 qd = qv/12.0d0 + qsea/6.0d0
642 qd = qd*x
643 ZDV=qd
644 ZDB=qd
645 ZSB=qd
646 g = WHIT3G(x,Q2)
647 g = g*x
648 ZGL=g
649c
650 if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
651 A0dcv= s*( 1.219000d-01)+s2*( 6.200000d+00)
652 $ +s3*(-2.504000d+01)+s4*( 3.098000d+01)
653 A1dcv= s*( 1.913000d+00)+s2*(-7.690000d+01)
654 $ +s3*( 3.180000d+02)+s4*(-3.920000d+02)
655 A2dcv= s*(-7.160000d+00)+s2*( 2.503000d+02)
656 $ +s3*(-1.062000d+03)+s4*( 1.308000d+03)
657 A3dcv= s*( 3.190000d+00)+s2*(-2.301000d+02)
658 $ +s3*( 1.012000d+03)+s4*(-1.250000d+03)
659 Bdcv = 4.990000d-01+s*( 3.470000d+00)+s2*(-1.526000d+01)
660 $ +s3*( 1.967000d+01)
661 Cdcv = 3.290000d-01+s*( 8.240000d+00)+s2*(-3.800000d+01)
662 $ +s3*( 4.630000d+01)
663 Adcs = s*(-1.948000d-02)+s2*( 2.861000d-02)
664 $ +s3*(-2.036000d-02)
665 B0dcs=-4.130000d-01+s*(-4.390000d-01)+s2*( 1.810000d-01)
666 B1dcs= 5.190000d+00+s*(-7.400000d+00)+s2*( 3.400000d+00)
667 Cdcs = 2.359000d+00+s*( 9.770000d-01)+s2*(-7.730000d+00)
668 $ +s3*( 9.480000d+00)
669c
670 dcv = 1.0d0/alinv/x*
671 $ (A0dcv+x*A1dcv+x2*A2dcv+x2*x*A3dcv) * x**Bdcv * x1**Cdcv
672 dcs = 1.0d0/alinv/x*
673 $ Adcs * x**(B0dcs+B1dcs*x) * x1**Cdcs
674c
675 call WHIT3Q(x,mc*mc/Q2,cv,cs)
676 qc = cv/alinv/2.0d0/PI + cs*alstpi + dcs + dcv
677 qc = qc*x
678 ZCB=qc
679 else
680 qc = 0.0d0
681 ZCB=qc
682 endif
683 endif
684c
685 return
686 end
687c-------------------------------------------------------
688 subroutine SFWHI4(ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL)
689c-------------------------------------------------------
690c WHIT4 parton distribution in the photon
691c
692c INPUT: integer ic : if ic=0 then qc=0
693c else qc is calculated
694c DOUBLE PRECISION Q2 : energy scale Q^2 (GeV^2)
695c DOUBLE PRECISION x : energy fraction
696c
697c OUTPUT: DOUBLE PRECISION qu : up-quark dist.
698c DOUBLE PRECISION qd : down- or strange-quark dist.
699c DOUBLE PRECISION qc : charm-quark dist.
700c DOUBLE PRECISION g : gluon dist.
701c-------------------------------------------------------
702c Modified by M.Tanaka on July 22, 1994.
703c The bug pointed out by M.Drees is fixed.
704c-------------------------------------------------------
705c Modified by I.Watanabe on July 22, 1994.
706c-------------------------------------------------------
707 implicit none
708 external WHIT4G
709 double precision
710 + ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL
711c arg
712 integer ic
713 DOUBLE PRECISION Q2,x
714 DOUBLE PRECISION qu,qd,qc,g
715c const
716 DOUBLE PRECISION q42it,q52it,lam42,lam52
717 DOUBLE PRECISION alinv,mc,PI
718c local
719 DOUBLE PRECISION qv,qsea,cv,cs,dcv,dcs
720 DOUBLE PRECISION A0val,A1val,A2val,Bval,Cval,
721 $ A0sea,B0sea,BB0sea,C0sea
722 DOUBLE PRECISION A0dcv,A1dcv,A2dcv,A3dcv,Bdcv,Cdcv
723 DOUBLE PRECISION Adcs, B0dcs, B1dcs, Cdcs
724 DOUBLE PRECISION x1,x2,mc2q2
725 DOUBLE PRECISION s,s2,s3,s4,prsccf,alstpi
726 DOUBLE PRECISION WHIT4G
727c parameters
728 parameter(lam42=0.16d0, lam52=0.091411319d0)
729 parameter(Q42IT=4.0d0, Q52IT=100.0d0)
730 parameter(alinv=137.036d0, mc=1.5d0)
731 parameter(pi=3.14159265358979323846d0)
732 common /scale/ s,s2,s3,s4,prsccf
733c
734c begin
735 x=ZX
736 Q2=ZQ*ZQ
737 ic=1
738c
739 x1=1.0d0-x
740 x2=x**2
741 mc2q2=mc**2/Q2
742c
743 if(Q2.lt.100.0d0) then
744c under 100 GeV^2
745c
746c set scale s
747 if(Q2.lt.4.0d0) then
748cccc for under 4GeV^2 prescription
749 s= 0.0d0
750 prsccf = log(Q2/LAM42)/ log(Q42IT/LAM42)
751 alstpi = 6.0d0/25.0d0/ log(Q42IT/LAM42)
752 else
753 s= log( log(Q2/LAM42)/ log(Q42IT/LAM42))
754 prsccf = 1.0d0
755 alstpi = 6.0d0/25.0d0/ log(Q2/LAM42)
756 endif
757 s2=s**2
758 s3=s2*s
759 s4=s2**2
760c
761cccccc WHIT4 quark (U100)
762c
763 A0val= 2.540000d+00+s*( 2.000000d+00)+s2*( 7.180000d-01)
764 A1val= 6.230000d-02+s*(-7.010000d+00)+s2*( 1.251000d-01)
765 A2val=-1.642000d-01+s*(-4.360000d-01)+s2*( 1.048000d+01)
766 $ +s3*(-5.200000d+00)
767 Bval = 6.990000d-01+s*(-2.796000d-02)+s2*(-3.650000d-03)
768 Cval = 4.420000d-01+s*(-1.255000d+00)+s2*( 1.941000d+00)
769 $ +s3*(-9.950000d-01)
770 A0sea= 1.308000d+00+s*( 2.315000d+00)+s2*(-7.880000d+00)
771 $ +s3*( 8.260000d+00)+s4*(-3.004000d+00)
772 B0sea=-3.730000d-02+s*( 5.630000d-02)+s2*(-1.133000d+00)
773 $ +s3*( 1.185000d+00)+s4*(-4.180000d-01)
774 BB0sea=2.103000d+00+s*( 4.850000d+00)+s2*(-1.781000d+01)
775 $ +s3*( 2.062000d+01)+s4*(-7.940000d+00)
776 C0sea= 7.000000d+00+s*(-1.017000d+01)+s2*( 2.600000d+01)
777 $ +s3*(-2.960000d+01)+s4*( 1.227000d+01)
778c
779 qv = prsccf/alinv/x*
780 $ (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
781 qsea= prsccf/alinv/x*
782 $ A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
783c
784 qu = qv/3.0d0 + qsea/6.0d0
785 qu = qu*x
786 ZUV=qu
787 ZUB=qu
788 qd = qv/12.0d0 + qsea/6.0d0
789 qd = qd*x
790 ZDV=qd
791 ZDB=qd
792 ZSB=qd
793c
794 if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
795 call WHIT4Q(x,mc2q2,cv,cs)
796 qc = cv/alinv/2.0d0/PI + cs*alstpi
797 qc = qc*x
798 ZCB=qc
799 else
800 qc = 0.0d0
801 ZCB=qc
802 endif
803c
804 g = WHIT4G(x,Q2)
805 g = g*x
806 ZGL=g
807c
808 else
809c over 100 GeV^2
810c
811c set scale s
812 s= log( log(Q2/LAM52)/ log(Q52IT/LAM52))
813 prsccf = 1.0d0
814 alstpi = 6.0d0/23.0d0/ log(Q2/LAM52)
815 s2=s**2
816 s3=s2*s
817 s4=s2**2
818c
819cccccc WHIT4 quark (O100)
820c
821 A0val= 4.270000d+00+s*( 3.096000d+00)+s2*( 1.619000d+00)
822 A1val=-4.740000d+00+s*(-6.900000d+00)+s2*(-2.430000d+00)
823 A2val= 2.837000d+00+s*( 6.470000d+00)+s2*( 4.090000d+00)
824 Bval = 6.780000d-01+s*(-3.940000d-02)+s2*( 1.756000d-02)
825 Cval = 1.728000d-01+s*(-2.479000d-02)+s2*( 1.446000d-01)
826 A0sea= 1.188000d+00+s*(-1.396000d+00)+s2*( 8.710000d+00)
827 $ +s3*(-2.542000d+01)+s4*( 2.492000d+01)
828 B0sea=-2.448000d-01+s*(-4.190000d-01)+s2*( 1.007000d+00)
829 $ +s3*(-2.689000d+00)+s4*( 2.517000d+00)
830 BB0sea=1.942000d+00+s*(-6.040000d+00)+s2*( 5.030000d+01)
831 $ +s3*(-1.478000d+02)+s4*( 1.481000d+02)
832 C0sea= 5.420000d+00+s*( 6.110000d+00)+s2*(-5.380000d+01)
833 $ +s3*( 1.632000d+02)+s4*(-1.716000d+02)
834c
835 qv = 1.0d0/alinv/x*
836 $ (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
837 qsea= 1.0d0/alinv/x*
838 $ A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
839c
840 qu = qv/3.0d0 + qsea/6.0d0
841 qu = qu*x
842 ZUV=qu
843 ZUB=qu
844 qd = qv/12.0d0 + qsea/6.0d0
845 qd = qd*x
846 ZDV=qd
847 ZDB=qd
848 ZSB=qd
849 g = WHIT4G(x,Q2)
850 g = g*x
851 ZGL=g
852c
853 if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
854 A0dcv= s*( 1.219000d-01)+s2*( 6.200000d+00)
855 $ +s3*(-2.504000d+01)+s4*( 3.098000d+01)
856 A1dcv= s*( 1.913000d+00)+s2*(-7.690000d+01)
857 $ +s3*( 3.180000d+02)+s4*(-3.920000d+02)
858 A2dcv= s*(-7.160000d+00)+s2*( 2.503000d+02)
859 $ +s3*(-1.062000d+03)+s4*( 1.308000d+03)
860 A3dcv= s*( 3.190000d+00)+s2*(-2.301000d+02)
861 $ +s3*( 1.012000d+03)+s4*(-1.250000d+03)
862 Bdcv = 4.990000d-01+s*( 3.470000d+00)+s2*(-1.526000d+01)
863 $ +s3*( 1.967000d+01)
864 Cdcv = 3.290000d-01+s*( 8.240000d+00)+s2*(-3.800000d+01)
865 $ +s3*( 4.630000d+01)
866 Adcs = s*(-2.821000d-02)+s2*(-2.649000d-04)
867 $ +s3*( 7.040000d-03)
868 B0dcs=-3.270000d-01+s*(-2.298000d-01)+s2*( 3.500000d-02)
869 B1dcs= 1.254000d+00+s*( 8.780000d-01)+s2*( 2.086000d-01)
870 Cdcs = 4.170000d+00+s*( 6.400000d-01)+s2*(-7.630000d+00)
871 $ +s3*( 7.170000d+00)
872c
873 dcv = 1.0d0/alinv/x*
874 $ (A0dcv+x*A1dcv+x2*A2dcv+x2*x*A3dcv) * x**Bdcv * x1**Cdcv
875 dcs = 1.0d0/alinv/x*
876 $ Adcs * x**(B0dcs+B1dcs*x) * x1**Cdcs
877c
878 call WHIT4Q(x,mc*mc/Q2,cv,cs)
879 qc = cv/alinv/2.0d0/PI + cs*alstpi + dcs + dcv
880 qc = qc*x
881 ZCB=qc
882 else
883 qc = 0.0d0
884 ZCB=qc
885 endif
886 endif
887c
888 return
889 end
890c-------------------------------------------------------
891 subroutine SFWHI5(ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL)
892c-------------------------------------------------------
893c WHIT5 parton distribution in the photon
894c
895c INPUT: integer ic : if ic=0 then qc=0
896c else qc is calculated
897c DOUBLE PRECISION Q2 : energy scale Q^2 (GeV^2)
898c DOUBLE PRECISION x : energy fraction
899c
900c OUTPUT: DOUBLE PRECISION qu : up-quark dist.
901c DOUBLE PRECISION qd : down- or strange-quark dist.
902c DOUBLE PRECISION qc : charm-quark dist.
903c DOUBLE PRECISION g : gluon dist.
904c-------------------------------------------------------
905c Modified by M.Tanaka on July 22, 1994.
906c The bug pointed out by M.Drees is fixed.
907c-------------------------------------------------------
908c Modified by I.Watanabe on July 22, 1994.
909c-------------------------------------------------------
910 implicit none
911 external WHIT5G
912 double precision
913 + ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL
914c arg
915 integer ic
916 DOUBLE PRECISION Q2,x
917 DOUBLE PRECISION qu,qd,qc,g
918c const
919 DOUBLE PRECISION q42it,q52it,lam42,lam52
920 DOUBLE PRECISION alinv,mc,PI
921c local
922 DOUBLE PRECISION qv,qsea,cv,cs,dcv,dcs
923 DOUBLE PRECISION A0val,A1val,A2val,Bval,Cval,
924 $ A0sea,B0sea,BB0sea,C0sea
925 DOUBLE PRECISION A0dcv,A1dcv,A2dcv,A3dcv,Bdcv,Cdcv
926 DOUBLE PRECISION Adcs, B0dcs, B1dcs, Cdcs
927 DOUBLE PRECISION x1,x2,mc2q2
928 DOUBLE PRECISION s,s2,s3,s4,prsccf,alstpi
929 DOUBLE PRECISION WHIT5G
930c parameters
931 parameter(lam42=0.16d0, lam52=0.091411319d0)
932 parameter(Q42IT=4.0d0, Q52IT=100.0d0)
933 parameter(alinv=137.036d0, mc=1.5d0)
934 parameter(pi=3.14159265358979323846d0)
935 common /scale/ s,s2,s3,s4,prsccf
936c
937c begin
938 x=ZX
939 Q2=ZQ*ZQ
940 ic=1
941c
942 x1=1.0d0-x
943 x2=x**2
944 mc2q2=mc**2/Q2
945c
946 if(Q2.lt.100.0d0) then
947c under 100 GeV^2
948c
949c set scale s
950 if(Q2.lt.4.0d0) then
951cccc for under 4GeV^2 prescription
952 s= 0.0d0
953 prsccf = log(Q2/LAM42)/ log(Q42IT/LAM42)
954 alstpi = 6.0d0/25.0d0/ log(Q42IT/LAM42)
955 else
956 s= log( log(Q2/LAM42)/ log(Q42IT/LAM42))
957 prsccf = 1.0d0
958 alstpi = 6.0d0/25.0d0/ log(Q2/LAM42)
959 endif
960 s2=s**2
961 s3=s2*s
962 s4=s2**2
963c
964cccccc WHIT5 quark (U100)
965c
966 A0val= 2.540000d+00+s*( 2.000000d+00)+s2*( 7.180000d-01)
967 A1val= 6.230000d-02+s*(-7.010000d+00)+s2*( 1.251000d-01)
968 A2val=-1.642000d-01+s*(-4.360000d-01)+s2*( 1.048000d+01)
969 $ +s3*(-5.200000d+00)
970 Bval = 6.990000d-01+s*(-2.796000d-02)+s2*(-3.650000d-03)
971 Cval = 4.420000d-01+s*(-1.255000d+00)+s2*( 1.941000d+00)
972 $ +s3*(-9.950000d-01)
973 A0sea= 2.227000d+00+s*( 5.720000d+00)+s2*(-1.295000d+01)
974 $ +s3*( 7.220000d+00)+s4*(-2.514000d-01)
975 B0sea=-8.810000d-02+s*( 1.465000d-01)+s2*(-9.750000d-01)
976 $ +s3*( 7.820000d-01)+s4*(-2.074000d-01)
977 BB0sea=3.370000d+00+s*( 1.416000d+01)+s2*(-3.150000d+01)
978 $ +s3*( 2.789000d+01)+s4*(-8.710000d+00)
979 C0sea= 1.581000d+01+s*(-3.630000d+01)+s2*( 7.710000d+01)
980 $ +s3*(-7.810000d+01)+s4*( 2.948000d+01)
981c
982 qv = prsccf/alinv/x*
983 $ (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
984 qsea= prsccf/alinv/x*
985 $ A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
986c
987 qu = qv/3.0d0 + qsea/6.0d0
988 qu = qu*x
989 ZUV=qu
990 ZUB=qu
991 qd = qv/12.0d0 + qsea/6.0d0
992 qd = qd*x
993 ZDV=qd
994 ZDB=qd
995 ZSB=qd
996c
997 if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
998 call WHIT5Q(x,mc2q2,cv,cs)
999 qc = cv/alinv/2.0d0/PI + cs*alstpi
1000 qc = qc*x
1001 ZCB=qc
1002 else
1003 qc = 0.0d0
1004 ZCB=qc
1005 endif
1006c
1007 g = WHIT5G(x,Q2)
1008 g = g*x
1009 ZGL=g
1010c
1011 else
1012c over 100 GeV^2
1013c
1014c set scale s
1015 s= log( log(Q2/LAM52)/ log(Q52IT/LAM52))
1016 prsccf = 1.0d0
1017 alstpi = 6.0d0/23.0d0/ log(Q2/LAM52)
1018 s2=s**2
1019 s3=s2*s
1020 s4=s2**2
1021c
1022cccccc WHIT5 quark (O100)
1023c
1024 A0val= 4.270000d+00+s*( 3.096000d+00)+s2*( 1.617000d+00)
1025 A1val=-4.740000d+00+s*(-6.900000d+00)+s2*(-2.417000d+00)
1026 A2val= 2.837000d+00+s*( 6.470000d+00)+s2*( 4.070000d+00)
1027 Bval = 6.780000d-01+s*(-3.940000d-02)+s2*( 1.750000d-02)
1028 Cval = 1.728000d-01+s*(-2.457000d-02)+s2*( 1.440000d-01)
1029 A0sea= 2.318000d+00+s*(-3.760000d+00)+s2*( 2.026000d+01)
1030 $ +s3*(-5.950000d+01)+s4*( 5.900000d+01)
1031 B0sea=-2.425000d-01+s*(-4.360000d-01)+s2*( 1.241000d+00)
1032 $ +s3*(-3.510000d+00)+s4*( 3.360000d+00)
1033 BB0sea=5.330000d+00+s*(-8.680000d+00)+s2*( 7.420000d+01)
1034 $ +s3*(-2.070000d+02)+s4*( 1.967000d+02)
1035 C0sea= 8.480000d+00+s*( 9.310000d+00)+s2*(-1.041000d+02)
1036 $ +s3*( 2.801000d+02)+s4*(-2.663000d+02)
1037c
1038 qv = 1.0d0/alinv/x*
1039 $ (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
1040 qsea= 1.0d0/alinv/x*
1041 $ A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
1042c
1043 qu = qv/3.0d0 + qsea/6.0d0
1044 qu = qu*x
1045 ZUV=qu
1046 ZUB=qu
1047 qd = qv/12.0d0 + qsea/6.0d0
1048 qd = qd*x
1049 ZDV=qd
1050 ZDB=qd
1051 ZSB=qd
1052 g = WHIT5G(x,Q2)
1053 g = g*x
1054 ZGL=g
1055c
1056 if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
1057 A0dcv= s*( 1.219000d-01)+s2*( 6.200000d+00)
1058 $ +s3*(-2.504000d+01)+s4*( 3.098000d+01)
1059 A1dcv= s*( 1.913000d+00)+s2*(-7.690000d+01)
1060 $ +s3*( 3.180000d+02)+s4*(-3.920000d+02)
1061 A2dcv= s*(-7.160000d+00)+s2*( 2.503000d+02)
1062 $ +s3*(-1.062000d+03)+s4*( 1.308000d+03)
1063 A3dcv= s*( 3.190000d+00)+s2*(-2.301000d+02)
1064 $ +s3*( 1.012000d+03)+s4*(-1.250000d+03)
1065 Bdcv = 4.990000d-01+s*( 3.470000d+00)+s2*(-1.526000d+01)
1066 $ +s3*( 1.967000d+01)
1067 Cdcv = 3.290000d-01+s*( 8.240000d+00)+s2*(-3.800000d+01)
1068 $ +s3*( 4.630000d+01)
1069 Adcs = s*(-6.580000d-02)+s2*( 1.059000d-01)
1070 $ +s3*(-6.630000d-02)
1071 B0dcs=-2.750000d-01+s*(-4.760000d-01)+s2*( 1.191000d-01)
1072 B1dcs= 6.370000d+00+s*(-5.320000d+00)+s2*( 1.986000d+00)
1073 Cdcs = 3.400000d+00+s*( 3.750000d-01)+s2*(-8.790000d+00)
1074 $ +s3*( 1.001000d+01)
1075c
1076 dcv = 1.0d0/alinv/x*
1077 $ (A0dcv+x*A1dcv+x2*A2dcv+x2*x*A3dcv) * x**Bdcv * x1**Cdcv
1078 dcs = 1.0d0/alinv/x*
1079 $ Adcs * x**(B0dcs+B1dcs*x) * x1**Cdcs
1080c
1081 call WHIT5Q(x,mc*mc/Q2,cv,cs)
1082 qc = cv/alinv/2.0d0/PI + cs*alstpi + dcs + dcv
1083 qc = qc*x
1084 ZCB=qc
1085 else
1086 qc = 0.0d0
1087 ZCB=qc
1088 endif
1089 endif
1090c
1091 return
1092 end
1093c-------------------------------------------------------
1094 subroutine SFWHI6(ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL)
1095c-------------------------------------------------------
1096c WHIT6 parton distribution in the photon
1097c
1098c INPUT: integer ic : if ic=0 then qc=0
1099c else qc is calculated
1100c DOUBLE PRECISION Q2 : energy scale Q^2 (GeV^2)
1101c DOUBLE PRECISION x : energy fraction
1102c
1103c OUTPUT: DOUBLE PRECISION qu : up-quark dist.
1104c DOUBLE PRECISION qd : down- or strange-quark dist.
1105c DOUBLE PRECISION qc : charm-quark dist.
1106c DOUBLE PRECISION g : gluon dist.
1107c-------------------------------------------------------
1108c Modified by M.Tanaka on July 22, 1994.
1109c The bug pointed out by M.Drees is fixed.
1110c-------------------------------------------------------
1111c Modified by I.Watanabe on July 22, 1994.
1112c-------------------------------------------------------
1113 implicit none
1114 external WHIT6G
1115 double precision
1116 + ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZGL
1117c arg
1118 integer ic
1119 DOUBLE PRECISION Q2,x
1120 DOUBLE PRECISION qu,qd,qc,g
1121c const
1122 DOUBLE PRECISION q42it,q52it,lam42,lam52
1123 DOUBLE PRECISION alinv,mc,PI
1124c local
1125 DOUBLE PRECISION qv,qsea,cv,cs,dcv,dcs
1126 DOUBLE PRECISION A0val,A1val,A2val,Bval,Cval,
1127 $ A0sea,B0sea,BB0sea,C0sea
1128 DOUBLE PRECISION A0dcv,A1dcv,A2dcv,A3dcv,Bdcv,Cdcv
1129 DOUBLE PRECISION Adcs, B0dcs, B1dcs, Cdcs
1130 DOUBLE PRECISION x1,x2,mc2q2
1131 DOUBLE PRECISION s,s2,s3,s4,prsccf,alstpi
1132 DOUBLE PRECISION WHIT6G
1133c parameters
1134 parameter(lam42=0.16d0, lam52=0.091411319d0)
1135 parameter(Q42IT=4.0d0, Q52IT=100.0d0)
1136 parameter(alinv=137.036d0, mc=1.5d0)
1137 parameter(pi=3.14159265358979323846d0)
1138 common /scale/ s,s2,s3,s4,prsccf
1139c
1140c begin
1141 x=ZX
1142 Q2=ZQ*ZQ
1143 ic=1
1144c
1145 x1=1.0d0-x
1146 x2=x**2
1147 mc2q2=mc**2/Q2
1148c
1149 if(Q2.lt.100.0d0) then
1150c under 100 GeV^2
1151c
1152c set scale s
1153 if(Q2.lt.4.0d0) then
1154cccc for under 4GeV^2 prescription
1155 s= 0.0d0
1156 prsccf = log(Q2/LAM42)/ log(Q42IT/LAM42)
1157 alstpi = 6.0d0/25.0d0/ log(Q42IT/LAM42)
1158 else
1159 s= log( log(Q2/LAM42)/ log(Q42IT/LAM42))
1160 prsccf = 1.0d0
1161 alstpi = 6.0d0/25.0d0/ log(Q2/LAM42)
1162 endif
1163 s2=s**2
1164 s3=s2*s
1165 s4=s2**2
1166c
1167cccccc WHIT6 quark (U100)
1168c
1169 A0val= 2.540000d+00+s*( 2.000000d+00)+s2*( 7.180000d-01)
1170 A1val= 6.230000d-02+s*(-7.010000d+00)+s2*( 1.251000d-01)
1171 A2val=-1.642000d-01+s*(-4.360000d-01)+s2*( 1.048000d+01)
1172 $ +s3*(-5.200000d+00)
1173 Bval = 6.990000d-01+s*(-2.796000d-02)+s2*(-3.650000d-03)
1174 Cval = 4.420000d-01+s*(-1.255000d+00)+s2*( 1.941000d+00)
1175 $ +s3*(-9.950000d-01)
1176 A0sea= 3.180000d+00+s*( 8.690000d+00)+s2*(-2.287000d+01)
1177 $ +s3*( 1.896000d+01)+s4*(-5.140000d+00)
1178 B0sea=-1.003000d-01+s*( 1.603000d-01)+s2*(-1.037000d+00)
1179 $ +s3*( 9.440000d-01)+s4*(-2.915000d-01)
1180 BB0sea=5.690000d+00+s*( 1.867000d+01)+s2*(-4.670000d+01)
1181 $ +s3*( 5.050000d+01)+s4*(-1.835000d+01)
1182 C0sea= 2.149000d+01+s*(-5.650000d+01)+s2*( 1.293000d+02)
1183 $ +s3*(-1.459000d+02)+s4*( 5.750000d+01)
1184c
1185 qv = prsccf/alinv/x*
1186 $ (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
1187 qsea= prsccf/alinv/x*
1188 $ A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
1189c
1190 qu = qv/3.0d0 + qsea/6.0d0
1191 qu = qu*x
1192 ZUV=qu
1193 ZUB=qu
1194 qd = qv/12.0d0 + qsea/6.0d0
1195 qd = qd*x
1196 ZDV=qd
1197 ZDB=qd
1198 ZSB=qd
1199c
1200 if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
1201 call WHIT6Q(x,mc2q2,cv,cs)
1202 qc = cv/alinv/2.0d0/PI + cs*alstpi
1203 qc = qc*x
1204 ZCB=qc
1205 else
1206 qc = 0.0d0
1207 ZCB=qc
1208 endif
1209c
1210 g = WHIT6G(x,Q2)
1211 g = g*x
1212 ZGL=g
1213c
1214 else
1215c over 100 GeV^2
1216c
1217c set scale s
1218 s= log( log(Q2/LAM52)/ log(Q52IT/LAM52))
1219 prsccf = 1.0d0
1220 alstpi = 6.0d0/23.0d0/ log(Q2/LAM52)
1221 s2=s**2
1222 s3=s2*s
1223 s4=s2**2
1224c
1225cccccc WHIT6 quark (O100)
1226c
1227 A0val= 4.270000d+00+s*( 3.096000d+00)+s2*( 1.621000d+00)
1228 A1val=-4.740000d+00+s*(-6.900000d+00)+s2*(-2.439000d+00)
1229 A2val= 2.837000d+00+s*( 6.460000d+00)+s2*( 4.100000d+00)
1230 Bval = 6.780000d-01+s*(-3.940000d-02)+s2*( 1.758000d-02)
1231 Cval = 1.728000d-01+s*(-2.493000d-02)+s2*( 1.451000d-01)
1232 A0sea= 3.340000d+00+s*(-5.610000d+00)+s2*( 5.000000d+01)
1233 $ +s3*(-2.207000d+02)+s4*( 3.028000d+02)
1234 B0sea=-2.402000d-01+s*(-4.090000d-01)+s2*( 2.263000d+00)
1235 $ +s3*(-1.050000d+01)+s4*( 1.487000d+01)
1236 BB0sea=8.790000d+00+s*(-8.860000d+00)+s2*( 1.640000d+02)
1237 $ +s3*(-7.120000d+02)+s4*( 9.730000d+02)
1238 C0sea= 9.160000d+00+s*( 9.290000d+00)+s2*(-2.784000d+02)
1239 $ +s3*( 1.175000d+03)+s4*(-1.592000d+03)
1240c
1241 qv = 1.0d0/alinv/x*
1242 $ (A0val+A1val*x+A2val*x2) * x**Bval * x1**Cval
1243 qsea= 1.0d0/alinv/x*
1244 $ A0sea * x**(B0sea+BB0sea*x) * x1**C0sea
1245c
1246 qu = qv/3.0d0 + qsea/6.0d0
1247 qu = qu*x
1248 ZUV=qu
1249 ZUB=qu
1250 qd = qv/12.0d0 + qsea/6.0d0
1251 qd = qd*x
1252 ZDV=qd
1253 ZDB=qd
1254 ZSB=qd
1255 g = WHIT6G(x,Q2)
1256 g = g*x
1257 ZGL=g
1258c
1259 if((ic.ne.0) .and. (x*(1.0d0+4.0d0*mc2q2).lt.1.0d0)) then
1260 A0dcv= s*( 1.219000d-01)+s2*( 6.200000d+00)
1261 $ +s3*(-2.504000d+01)+s4*( 3.098000d+01)
1262 A1dcv= s*( 1.913000d+00)+s2*(-7.690000d+01)
1263 $ +s3*( 3.180000d+02)+s4*(-3.920000d+02)
1264 A2dcv= s*(-7.160000d+00)+s2*( 2.503000d+02)
1265 $ +s3*(-1.062000d+03)+s4*( 1.308000d+03)
1266 A3dcv= s*( 3.190000d+00)+s2*(-2.301000d+02)
1267 $ +s3*( 1.012000d+03)+s4*(-1.250000d+03)
1268 Bdcv = 4.990000d-01+s*( 3.470000d+00)+s2*(-1.526000d+01)
1269 $ +s3*( 1.967000d+01)
1270 Cdcv = 3.290000d-01+s*( 8.240000d+00)+s2*(-3.800000d+01)
1271 $ +s3*( 4.630000d+01)
1272 Adcs = s*(-4.990000d-02)+s2*( 1.026000d-01)
1273 $ +s3*(-7.870000d-02)
1274 B0dcs=-3.610000d-01+s*(-5.760000d-01)+s2*( 2.257000d-01)
1275 B1dcs= 7.680000d+00+s*(-8.830000d+00)+s2*( 3.880000d+00)
1276 Cdcs = 2.548000d+00+s*( 6.910000d-01)+s2*(-8.700000d+00)
1277 $ +s3*( 1.065000d+01)
1278c
1279 dcv = 1.0d0/alinv/x*
1280 $ (A0dcv+x*A1dcv+x2*A2dcv+x2*x*A3dcv) * x**Bdcv * x1**Cdcv
1281 dcs = 1.0d0/alinv/x*
1282 $ Adcs * x**(B0dcs+B1dcs*x) * x1**Cdcs
1283c
1284 call WHIT6Q(x,mc*mc/Q2,cv,cs)
1285 qc = cv/alinv/2.0d0/PI + cs*alstpi + dcs + dcv
1286 qc = qc*x
1287 ZCB=qc
1288 else
1289 qc = 0.0d0
1290 ZCB=qc
1291 endif
1292 endif
1293c
1294 return
1295 end
1296c
1297ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1298 DOUBLE PRECISION function WHIT1G(x,Q2)
1299c input: x,Q2
1300c output: clg
1301c (gluon dist.)
1302ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1303 implicit none
1304c arg
1305 DOUBLE PRECISION Q2,x
1306c const
1307 DOUBLE PRECISION q42it,q52it,lam42,lam52
1308 DOUBLE PRECISION alinv
1309c local
1310 DOUBLE PRECISION A0g,B0g,C0g,A1g,AA1g,B1g,C1g
1311 DOUBLE PRECISION s,s2,s3,s4,prsccf
1312 DOUBLE PRECISION x1
1313c parameters
1314 parameter(lam42=0.16d0, lam52=0.091411319d0)
1315 parameter(Q42IT=4.0d0, Q52IT=100.0d0)
1316 parameter(alinv=137.036d0)
1317 common /scale/ s,s2,s3,s4,prsccf
1318c
1319c begin
1320 x1=1.0d0-x
1321c
1322 if(Q2.lt.100.0d0) then
1323c under 100 GeV^2
1324c
1325cccccc WHIT1 gluon (U100)
1326c
1327 A0g = 2.000000d+00+s*(-3.280000d+00)+s2*( 2.894000d+00)
1328 $ +s3*(-1.561000d+00)+s4*( 8.180000d-01)
1329 B0g = s*(-7.610000d-01)+s2*(-4.900000d-02)
1330 $ +s3*( 4.460000d-01)
1331 C0g = 3.000000d+00+s*( 1.586000d+00)+s2*(-9.490000d-01)
1332 $ +s3*( 2.425000d+00)
1333 A1g = s*( 4.610000d-01)+s2*( 1.041000d-01)
1334 $ +s3*(-1.753000d-02)+s4*(-2.717000d-01)
1335 AA1g= s*( 9.680000d-03)+s2*(-4.170000d-01)
1336 $ +s3*(-3.950000d-01)+s4*( 8.430000d-01)
1337 B1g =-4.140000d-01+s*(-6.060000d-02)+s2*( 2.847000d-01)
1338 $ +s3*(-5.070000d-01)
1339 C1g = 1.244000d+00+s*( 5.880000d-01)+s2*(-1.228000d+00)
1340 $ +s3*( 8.090000d-01)
1341 else
1342c over 100 GeV^2
1343c
1344cccccc WHIT1 gluon (O100)
1345c
1346 A0g = 7.840000d-01+s*(-2.238000d+00)+s2*( 1.617000d+01)
1347 $ +s3*(-6.250000d+01)+s4*( 8.390000d+01)
1348 B0g =-4.030000d-01+s*(-1.307000d+00)+s2*( 8.780000d+00)
1349 $ +s3*(-3.580000d+01)+s4*( 5.350000d+01)
1350 C0g = 4.450000d+00+s*( 1.027000d+00)+s2*( 4.460000d+01)
1351 $ +s3*(-1.600000d+02)+s4*( 1.816000d+02)
1352 A1g = 3.010000d-01+s*( 1.275000d+00)+s2*(-1.563000d+00)
1353 $ +s3*( 4.100000d+00)+s4*(-1.337000d+01)
1354 AA1g=-1.305000d-01+s*(-1.245000d+00)+s2*( 2.438000d+00)
1355 $ +s3*(-2.539000d+00)+s4*( 1.273000d+01)
1356 B1g =-4.890000d-01+s*( 9.550000d-01)+s2*(-4.400000d+00)
1357 $ +s3*( 1.022000d+01)+s4*(-1.713000d+01)
1358 C1g = 1.331000d+00+s*(-2.481000d-01)+s2*( 1.950000d+00)
1359 $ +s3*(-2.072000d+00)
1360 endif
1361c
1362 WHIT1G = prsccf/alinv/x*
1363 $ ( A0g * x**B0g * x1**C0g
1364 $ +(A1g+AA1g*x) * x**B1g * x1**C1g )
1365c
1366 return
1367 end
1368c
1369cccccccccccccccccccccccccccccc
1370c QPM calculation
1371 subroutine WHIT1Q(x,r,cv,cs)
1372ccc INPUTS : x,r=mc^2/Q^2
1373ccc OUTPUTS: cv,cs (valence- and sea- charm quark dist)
1374ccc cv <-- cv / ( alpha / 2PI)
1375ccc cs <-- cs / ( alpha_s/2PI)
1376c
1377 implicit none
1378c arg
1379 DOUBLE PRECISION x,r
1380 DOUBLE PRECISION cv,cs
1381c CONST
1382 DOUBLE PRECISION ec,mc
1383 parameter(ec=2.0d0/3.0d0,mc=1.5d0)
1384c N=15 Gauss int. weights and points
1385 integer GN,i
1386 parameter(GN=15)
1387 DOUBLE PRECISION XG(GN), XW(GN)
1388 DATA (XG(i),i=1,GN)/6.003741d-03,
1389 $ 3.136330d-02, 7.589671d-02, 1.377911d-01, 2.145139d-01,
1390 $ 3.029243d-01, 3.994030d-01, 5.000000d-01, 6.005970d-01,
1391 $ 6.970757d-01, 7.854861d-01, 8.622089d-01, 9.241033d-01,
1392 $ 9.686367d-01, 9.939963d-01/
1393 DATA (XW(i),i=1,GN)/1.537662d-02,
1394 $ 3.518302d-02, 5.357961d-02, 6.978534d-02, 8.313460d-02,
1395 $ 9.308050d-02, 9.921574d-02, 1.012891d-01, 9.921574d-02,
1396 $ 9.308050d-02, 8.313460d-02, 6.978534d-02, 5.357961d-02,
1397 $ 3.518302d-02, 1.537662d-02/
1398c local
1399 DOUBLE PRECISION sum, y,z,beta,w,WHIT1G,L
1400 parameter(L=4.0d0)
1401 DOUBLE PRECISION x1,rx,z1,rz
1402c
1403c begin
1404 x1=1.0d0-x
1405 rx=4.0d0*r*x
1406c
1407c direct
1408 beta=dsqrt(1.0d0-rx/x1)
1409 w=x*( beta*(-1.0d0+8.0d0*x*x1-rx*x1)
1410 $ +(x**2+x1**2+rx*(1.0d0-3.0d0*x)-0.5d0*rx**2)
1411 $ * log( (1.0d0+beta)/(1.0d0-beta) ))
1412 cv = 3.0d0*ec**2 * w / x
1413c
1414c resolved
1415 sum=0.0d0
1416 do 10 i=1,GN
1417 y= x+rx + (x1-rx)*XG(i)**L
1418 z=x/y
1419 z1=1.0d0-z
1420 rz=4.0d0*r*z
1421 beta=dsqrt(1.0d0-rz/z1)
1422 w=z*( beta*(-1.0d0+8.0d0*z*z1-rz*z1)
1423 $ +(z**2+z1**2+rz*(1.0d0-3.0d0*z)-0.5d0*rz**2)
1424 $ * log( (1.0d0+beta)/(1.0d0-beta) ))
1425c
1426 sum= sum + w * WHIT1G(y,mc**2/r)* L*XG(i)**(L-1.0d0)*XW(i)
1427c
1428 10 continue
1429c
1430 cs = 0.5d0/x * (x1-rx) * sum
1431c
1432 return
1433 end
1434c
1435ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1436 DOUBLE PRECISION function WHIT2G(x,Q2)
1437c input: x,Q2
1438c output: clg
1439c (gluon dist.)
1440ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1441 implicit none
1442c arg
1443 DOUBLE PRECISION Q2,x
1444c const
1445 DOUBLE PRECISION q42it,q52it,lam42,lam52
1446 DOUBLE PRECISION alinv
1447c local
1448 DOUBLE PRECISION A0g,B0g,C0g,A1g,AA1g,B1g,C1g
1449 DOUBLE PRECISION s,s2,s3,s4,prsccf
1450 DOUBLE PRECISION x1
1451c parameters
1452 parameter(lam42=0.16d0, lam52=0.091411319d0)
1453 parameter(Q42IT=4.0d0, Q52IT=100.0d0)
1454 parameter(alinv=137.036d0)
1455 common /scale/ s,s2,s3,s4,prsccf
1456c
1457c begin
1458 x1=1.0d0-x
1459c
1460 if(Q2.lt.100.0d0) then
1461c under 100 GeV^2
1462c
1463cccccc WHIT2 gluon (U100)
1464c
1465 A0g = 5.000000d+00+s*(-1.499000d+01)+s2*( 2.617000d+01)
1466 $ +s3*(-2.530000d+01)+s4*( 1.012000d+01)
1467 B0g = s*(-9.370000d-01)+s2*( 4.100000d-01)
1468 $ +s3*( 3.390000d-02)
1469 C0g = 9.000000d+00+s*( 7.090000d-01)+s2*( 3.118000d+00)
1470 $ +s3*(-5.820000d-04)
1471 A1g = s*( 4.610000d-01)+s2*( 1.041000d-01)
1472 $ +s3*(-1.753000d-02)+s4*(-2.717000d-01)
1473 AA1g= s*( 9.680000d-03)+s2*(-4.170000d-01)
1474 $ +s3*(-3.950000d-01)+s4*( 8.430000d-01)
1475 B1g =-4.140000d-01+s*(-6.060000d-02)+s2*( 2.847000d-01)
1476 $ +s3*(-5.070000d-01)
1477 C1g = 1.244000d+00+s*( 5.880000d-01)+s2*(-1.228000d+00)
1478 $ +s3*( 8.090000d-01)
1479 else
1480c over 100 GeV^2
1481c
1482cccccc WHIT2 gluon (O100)
1483c
1484 A0g = 1.095000d+00+s*(-2.388000d+00)+s2*( 9.190000d+00)
1485 $ +s3*(-3.032000d+01)+s4*( 3.480000d+01)
1486 B0g =-4.410000d-01+s*(-9.070000d-01)+s2*( 4.680000d+00)
1487 $ +s3*(-1.866000d+01)+s4*( 2.717000d+01)
1488 C0g = 1.099000d+01+s*( 4.710000d+00)+s2*( 2.801000d+01)
1489 $ +s3*(-1.279000d+02)+s4*( 1.640000d+02)
1490 A1g = 3.010000d-01+s*( 1.275000d+00)+s2*(-1.563000d+00)
1491 $ +s3*( 4.100000d+00)+s4*(-1.337000d+01)
1492 AA1g=-1.305000d-01+s*(-1.245000d+00)+s2*( 2.438000d+00)
1493 $ +s3*(-2.539000d+00)+s4*( 1.273000d+01)
1494 B1g =-4.890000d-01+s*( 9.550000d-01)+s2*(-4.400000d+00)
1495 $ +s3*( 1.022000d+01)+s4*(-1.713000d+01)
1496 C1g = 1.331000d+00+s*(-2.481000d-01)+s2*( 1.950000d+00)
1497 $ +s3*(-2.072000d+00)
1498 endif
1499c
1500 WHIT2G = prsccf/alinv/x*
1501 $ ( A0g * x**B0g * x1**C0g
1502 $ +(A1g+AA1g*x) * x**B1g * x1**C1g )
1503c
1504 return
1505 end
1506c
1507cccccccccccccccccccccccccccccc
1508c QPM calculation
1509 subroutine WHIT2Q(x,r,cv,cs)
1510ccc INPUTS : x,r=mc^2/Q^2
1511ccc OUTPUTS: cv,cs (valence- and sea- charm quark dist)
1512ccc cv <-- cv / ( alpha / 2PI)
1513ccc cs <-- cs / ( alpha_s/2PI)
1514c
1515 implicit none
1516c arg
1517 DOUBLE PRECISION x,r
1518 DOUBLE PRECISION cv,cs
1519c CONST
1520 DOUBLE PRECISION ec,mc
1521 parameter(ec=2.0d0/3.0d0,mc=1.5d0)
1522c N=15 Gauss int. weights and points
1523 integer GN,i
1524 parameter(GN=15)
1525 DOUBLE PRECISION XG(GN), XW(GN)
1526 DATA (XG(i),i=1,GN)/6.003741d-03,
1527 $ 3.136330d-02, 7.589671d-02, 1.377911d-01, 2.145139d-01,
1528 $ 3.029243d-01, 3.994030d-01, 5.000000d-01, 6.005970d-01,
1529 $ 6.970757d-01, 7.854861d-01, 8.622089d-01, 9.241033d-01,
1530 $ 9.686367d-01, 9.939963d-01/
1531 DATA (XW(i),i=1,GN)/1.537662d-02,
1532 $ 3.518302d-02, 5.357961d-02, 6.978534d-02, 8.313460d-02,
1533 $ 9.308050d-02, 9.921574d-02, 1.012891d-01, 9.921574d-02,
1534 $ 9.308050d-02, 8.313460d-02, 6.978534d-02, 5.357961d-02,
1535 $ 3.518302d-02, 1.537662d-02/
1536c local
1537 DOUBLE PRECISION sum, y,z,beta,w,WHIT2G,L
1538 parameter(L=4.0d0)
1539 DOUBLE PRECISION x1,rx,z1,rz
1540c
1541c begin
1542 x1=1.0d0-x
1543 rx=4.0d0*r*x
1544c
1545c direct
1546 beta=dsqrt(1.0d0-rx/x1)
1547 w=x*( beta*(-1.0d0+8.0d0*x*x1-rx*x1)
1548 $ +(x**2+x1**2+rx*(1.0d0-3.0d0*x)-0.5d0*rx**2)
1549 $ * log( (1.0d0+beta)/(1.0d0-beta) ))
1550 cv = 3.0d0*ec**2 * w / x
1551c
1552c resolved
1553 sum=0.0d0
1554 do 10 i=1,GN
1555 y= x+rx + (x1-rx)*XG(i)**L
1556 z=x/y
1557 z1=1.0d0-z
1558 rz=4.0d0*r*z
1559 beta=dsqrt(1.0d0-rz/z1)
1560 w=z*( beta*(-1.0d0+8.0d0*z*z1-rz*z1)
1561 $ +(z**2+z1**2+rz*(1.0d0-3.0d0*z)-0.5d0*rz**2)
1562 $ * log( (1.0d0+beta)/(1.0d0-beta) ))
1563c
1564 sum= sum + w * WHIT2G(y,mc**2/r)* L*XG(i)**(L-1.0d0)*XW(i)
1565c
1566 10 continue
1567c
1568 cs = 0.5d0/x * (x1-rx) * sum
1569c
1570 return
1571 end
1572c
1573ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1574 DOUBLE PRECISION function WHIT3G(x,Q2)
1575c input: x,Q2
1576c output: clg
1577c (gluon dist.)
1578ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1579 implicit none
1580c arg
1581 DOUBLE PRECISION Q2,x
1582c const
1583 DOUBLE PRECISION q42it,q52it,lam42,lam52
1584 DOUBLE PRECISION alinv
1585c local
1586 DOUBLE PRECISION A0g,B0g,C0g,A1g,AA1g,B1g,C1g
1587 DOUBLE PRECISION s,s2,s3,s4,prsccf
1588 DOUBLE PRECISION x1
1589c parameters
1590 parameter(lam42=0.16d0, lam52=0.091411319d0)
1591 parameter(Q42IT=4.0d0, Q52IT=100.0d0)
1592 parameter(alinv=137.036d0)
1593 common /scale/ s,s2,s3,s4,prsccf
1594c
1595c begin
1596 x1=1.0d0-x
1597c
1598 if(Q2.lt.100.0d0) then
1599c under 100 GeV^2
1600c
1601cccccc WHIT3 gluon (U100)
1602c
1603 A0g = 8.000000d+00+s*(-2.864000d+01)+s2*( 5.590000d+01)
1604 $ +s3*(-5.760000d+01)+s4*( 2.366000d+01)
1605 B0g = s*(-9.870000d-01)+s2*( 5.100000d-01)
1606 $ +s3*(-6.670000d-02)
1607 C0g = 1.500000d+01+s*( 3.310000d-01)+s2*( 3.500000d+00)
1608 $ +s3*( 8.920000d-01)
1609 A1g = s*( 4.610000d-01)+s2*( 1.041000d-01)
1610 $ +s3*(-1.753000d-02)+s4*(-2.717000d-01)
1611 AA1g= s*( 9.680000d-03)+s2*(-4.170000d-01)
1612 $ +s3*(-3.950000d-01)+s4*( 8.430000d-01)
1613 B1g =-4.140000d-01+s*(-6.060000d-02)+s2*( 2.847000d-01)
1614 $ +s3*(-5.070000d-01)
1615 C1g = 1.244000d+00+s*( 5.880000d-01)+s2*(-1.228000d+00)
1616 $ +s3*( 8.090000d-01)
1617 else
1618c over 100 GeV^2
1619c
1620cccccc WHIT3 gluon (O100)
1621c
1622 A0g = 1.270000d+00+s*(-2.817000d+00)+s2*( 5.740000d+00)
1623 $ +s3*(-1.327000d+01)+s4*( 1.268000d+01)
1624 B0g =-4.610000d-01+s*(-8.170000d-01)+s2*( 3.320000d+00)
1625 $ +s3*(-1.296000d+01)+s4*( 1.893000d+01)
1626 C0g = 1.721000d+01+s*( 1.257000d+00)+s2*( 5.050000d+01)
1627 $ +s3*(-2.761000d+02)+s4*( 4.900000d+02)
1628 A1g = 3.010000d-01+s*( 1.275000d+00)+s2*(-1.563000d+00)
1629 $ +s3*( 4.100000d+00)+s4*(-1.337000d+01)
1630 AA1g=-1.305000d-01+s*(-1.245000d+00)+s2*( 2.438000d+00)
1631 $ +s3*(-2.539000d+00)+s4*( 1.273000d+01)
1632 B1g =-4.890000d-01+s*( 9.550000d-01)+s2*(-4.400000d+00)
1633 $ +s3*( 1.022000d+01)+s4*(-1.713000d+01)
1634 C1g = 1.331000d+00+s*(-2.481000d-01)+s2*( 1.950000d+00)
1635 $ +s3*(-2.072000d+00)
1636 endif
1637c
1638 WHIT3G = prsccf/alinv/x*
1639 $ ( A0g * x**B0g * x1**C0g
1640 $ +(A1g+AA1g*x) * x**B1g * x1**C1g )
1641c
1642 return
1643 end
1644c
1645cccccccccccccccccccccccccccccc
1646c QPM calculation
1647 subroutine WHIT3Q(x,r,cv,cs)
1648ccc INPUTS : x,r=mc^2/Q^2
1649ccc OUTPUTS: cv,cs (valence- and sea- charm quark dist)
1650ccc cv <-- cv / ( alpha / 2PI)
1651ccc cs <-- cs / ( alpha_s/2PI)
1652c
1653 implicit none
1654c arg
1655 DOUBLE PRECISION x,r
1656 DOUBLE PRECISION cv,cs
1657c CONST
1658 DOUBLE PRECISION ec,mc
1659 parameter(ec=2.0d0/3.0d0,mc=1.5d0)
1660c N=15 Gauss int. weights and points
1661 integer GN,i
1662 parameter(GN=15)
1663 DOUBLE PRECISION XG(GN), XW(GN)
1664 DATA (XG(i),i=1,GN)/6.003741d-03,
1665 $ 3.136330d-02, 7.589671d-02, 1.377911d-01, 2.145139d-01,
1666 $ 3.029243d-01, 3.994030d-01, 5.000000d-01, 6.005970d-01,
1667 $ 6.970757d-01, 7.854861d-01, 8.622089d-01, 9.241033d-01,
1668 $ 9.686367d-01, 9.939963d-01/
1669 DATA (XW(i),i=1,GN)/1.537662d-02,
1670 $ 3.518302d-02, 5.357961d-02, 6.978534d-02, 8.313460d-02,
1671 $ 9.308050d-02, 9.921574d-02, 1.012891d-01, 9.921574d-02,
1672 $ 9.308050d-02, 8.313460d-02, 6.978534d-02, 5.357961d-02,
1673 $ 3.518302d-02, 1.537662d-02/
1674c local
1675 DOUBLE PRECISION sum, y,z,beta,w,WHIT3G,L
1676 parameter(L=4.0d0)
1677 DOUBLE PRECISION x1,rx,z1,rz
1678c
1679c begin
1680 x1=1.0d0-x
1681 rx=4.0d0*r*x
1682c
1683c direct
1684 beta=dsqrt(1.0d0-rx/x1)
1685 w=x*( beta*(-1.0d0+8.0d0*x*x1-rx*x1)
1686 $ +(x**2+x1**2+rx*(1.0d0-3.0d0*x)-0.5d0*rx**2)
1687 $ * log( (1.0d0+beta)/(1.0d0-beta) ))
1688 cv = 3.0d0*ec**2 * w / x
1689c
1690c resolved
1691 sum=0.0d0
1692 do 10 i=1,GN
1693 y= x+rx + (x1-rx)*XG(i)**L
1694 z=x/y
1695 z1=1.0d0-z
1696 rz=4.0d0*r*z
1697 beta=dsqrt(1.0d0-rz/z1)
1698 w=z*( beta*(-1.0d0+8.0d0*z*z1-rz*z1)
1699 $ +(z**2+z1**2+rz*(1.0d0-3.0d0*z)-0.5d0*rz**2)
1700 $ * log( (1.0d0+beta)/(1.0d0-beta) ))
1701c
1702 sum= sum + w * WHIT3G(y,mc**2/r)* L*XG(i)**(L-1.0d0)*XW(i)
1703c
1704 10 continue
1705c
1706 cs = 0.5d0/x * (x1-rx) * sum
1707c
1708 return
1709 end
1710c
1711ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1712 DOUBLE PRECISION function WHIT4G(x,Q2)
1713c input: x,Q2
1714c output: clg
1715c (gluon dist.)
1716ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1717 implicit none
1718c arg
1719 DOUBLE PRECISION Q2,x
1720c const
1721 DOUBLE PRECISION q42it,q52it,lam42,lam52
1722 DOUBLE PRECISION alinv
1723c local
1724 DOUBLE PRECISION A0g,B0g,C0g,A1g,AA1g,B1g,C1g
1725 DOUBLE PRECISION s,s2,s3,s4,prsccf
1726 DOUBLE PRECISION x1
1727c parameters
1728 parameter(lam42=0.16d0, lam52=0.091411319d0)
1729 parameter(Q42IT=4.0d0, Q52IT=100.0d0)
1730 parameter(alinv=137.036d0)
1731 common /scale/ s,s2,s3,s4,prsccf
1732c
1733c begin
1734 x1=1.0d0-x
1735c
1736 if(Q2.lt.100.0d0) then
1737c under 100 GeV^2
1738c
1739cccccc WHIT4 gluon (U100)
1740c
1741 A0g = 4.000000d+00+s*(-9.400000d+00)+s2*( 1.555000d+01)
1742 $ +s3*(-1.450000d+01)+s4*( 5.470000d+00)
1743 B0g = s*(-1.142000d+00)+s2*( 1.034000d+00)
1744 $ +s3*(-4.410000d-01)
1745 C0g = 3.000000d+00+s*( 8.720000d-01)+s2*( 1.006000d+00)
1746 $ +s3*( 3.560000d-01)
1747 A1g = s*( 6.020000d-01)+s2*( 5.090000d-01)
1748 $ +s3*(-2.054000d+00)+s4*( 1.392000d+00)
1749 AA1g= s*(-9.220000d-02)+s2*(-1.899000d+00)
1750 $ +s3*( 4.180000d+00)+s4*(-2.494000d+00)
1751 B1g =-2.895000d-01+s*( 3.760000d-01)+s2*(-1.719000d+00)
1752 $ +s3*( 1.116000d+00)
1753 C1g = 1.439000d+00+s*(-5.570000d-01)+s2*( 3.660000d-01)
1754 $ +s3*( 7.330000d-01)+s4*(-7.620000d-01)
1755 else
1756c over 100 GeV^2
1757c
1758cccccc WHIT4 gluon (O100)
1759c
1760 A0g = 1.384000d+00+s*(-2.455000d+00)+s2*( 8.940000d+00)
1761 $ +s3*(-2.906000d+01)+s4*( 3.710000d+01)
1762 B0g =-4.420000d-01+s*(-7.190000d-01)+s2*( 2.961000d+00)
1763 $ +s3*(-1.209000d+01)+s4*( 1.916000d+01)
1764 C0g = 4.210000d+00+s*( 2.524000d+00)+s2*( 1.003000d+01)
1765 $ +s3*(-1.827000d+01)+s4*( 2.162000d+00)
1766 A1g = 2.992000d-01+s*( 1.179000d+00)+s2*(-1.915000d+00)
1767 $ +s3*( 7.260000d+00)+s4*(-1.839000d+01)
1768 AA1g=-1.600000d-01+s*(-1.114000d+00)+s2*( 2.939000d+00)
1769 $ +s3*(-6.660000d+00)+s4*( 1.923000d+01)
1770 B1g =-4.830000d-01+s*( 7.550000d-01)+s2*(-3.800000d+00)
1771 $ +s3*( 1.075000d+01)+s4*(-1.993000d+01)
1772 C1g = 1.297000d+00+s*(-1.669000d-01)+s2*( 1.906000d+00)
1773 $ +s3*(-2.057000d+00)
1774 endif
1775c
1776 WHIT4G = prsccf/alinv/x*
1777 $ ( A0g * x**B0g * x1**C0g
1778 $ +(A1g+AA1g*x) * x**B1g * x1**C1g )
1779c
1780 return
1781 end
1782c
1783cccccccccccccccccccccccccccccc
1784c QPM calculation
1785 subroutine WHIT4Q(x,r,cv,cs)
1786ccc INPUTS : x,r=mc^2/Q^2
1787ccc OUTPUTS: cv,cs (valence- and sea- charm quark dist)
1788ccc cv <-- cv / ( alpha / 2PI)
1789ccc cs <-- cs / ( alpha_s/2PI)
1790c
1791 implicit none
1792c arg
1793 DOUBLE PRECISION x,r
1794 DOUBLE PRECISION cv,cs
1795c CONST
1796 DOUBLE PRECISION ec,mc
1797 parameter(ec=2.0d0/3.0d0,mc=1.5d0)
1798c N=15 Gauss int. weights and points
1799 integer GN,i
1800 parameter(GN=15)
1801 DOUBLE PRECISION XG(GN), XW(GN)
1802 DATA (XG(i),i=1,GN)/6.003741d-03,
1803 $ 3.136330d-02, 7.589671d-02, 1.377911d-01, 2.145139d-01,
1804 $ 3.029243d-01, 3.994030d-01, 5.000000d-01, 6.005970d-01,
1805 $ 6.970757d-01, 7.854861d-01, 8.622089d-01, 9.241033d-01,
1806 $ 9.686367d-01, 9.939963d-01/
1807 DATA (XW(i),i=1,GN)/1.537662d-02,
1808 $ 3.518302d-02, 5.357961d-02, 6.978534d-02, 8.313460d-02,
1809 $ 9.308050d-02, 9.921574d-02, 1.012891d-01, 9.921574d-02,
1810 $ 9.308050d-02, 8.313460d-02, 6.978534d-02, 5.357961d-02,
1811 $ 3.518302d-02, 1.537662d-02/
1812c local
1813 DOUBLE PRECISION sum, y,z,beta,w,WHIT4G,L
1814 parameter(L=4.0d0)
1815 DOUBLE PRECISION x1,rx,z1,rz
1816c
1817c begin
1818 x1=1.0d0-x
1819 rx=4.0d0*r*x
1820c
1821c direct
1822 beta=dsqrt(1.0d0-rx/x1)
1823 w=x*( beta*(-1.0d0+8.0d0*x*x1-rx*x1)
1824 $ +(x**2+x1**2+rx*(1.0d0-3.0d0*x)-0.5d0*rx**2)
1825 $ * log( (1.0d0+beta)/(1.0d0-beta) ))
1826 cv = 3.0d0*ec**2 * w / x
1827c
1828c resolved
1829 sum=0.0d0
1830 do 10 i=1,GN
1831 y= x+rx + (x1-rx)*XG(i)**L
1832 z=x/y
1833 z1=1.0d0-z
1834 rz=4.0d0*r*z
1835 beta=dsqrt(1.0d0-rz/z1)
1836 w=z*( beta*(-1.0d0+8.0d0*z*z1-rz*z1)
1837 $ +(z**2+z1**2+rz*(1.0d0-3.0d0*z)-0.5d0*rz**2)
1838 $ * log( (1.0d0+beta)/(1.0d0-beta) ))
1839c
1840 sum= sum + w * WHIT4G(y,mc**2/r)* L*XG(i)**(L-1.0d0)*XW(i)
1841c
1842 10 continue
1843c
1844 cs = 0.5d0/x * (x1-rx) * sum
1845c
1846 return
1847 end
1848c
1849ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1850 DOUBLE PRECISION function WHIT5G(x,Q2)
1851c input: x,Q2
1852c output: clg
1853c (gluon dist.)
1854ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1855 implicit none
1856c arg
1857 DOUBLE PRECISION Q2,x
1858c const
1859 DOUBLE PRECISION q42it,q52it,lam42,lam52
1860 DOUBLE PRECISION alinv
1861c local
1862 DOUBLE PRECISION A0g,B0g,C0g,A1g,AA1g,B1g,C1g
1863 DOUBLE PRECISION s,s2,s3,s4,prsccf
1864 DOUBLE PRECISION x1
1865c parameters
1866 parameter(lam42=0.16d0, lam52=0.091411319d0)
1867 parameter(Q42IT=4.0d0, Q52IT=100.0d0)
1868 parameter(alinv=137.036d0)
1869 common /scale/ s,s2,s3,s4,prsccf
1870c
1871c begin
1872 x1=1.0d0-x
1873c
1874 if(Q2.lt.100.0d0) then
1875c under 100 GeV^2
1876c
1877cccccc WHIT5 gluon (U100)
1878c
1879 A0g = 1.000000d+01+s*(-3.400000d+01)+s2*( 6.900000d+01)
1880 $ +s3*(-7.530000d+01)+s4*( 3.230000d+01)
1881 B0g = s*(-1.126000d+00)+s2*( 9.260000d-01)
1882 $ +s3*(-3.930000d-01)
1883 C0g = 9.000000d+00+s*( 4.810000d-01)+s2*( 3.200000d+00)
1884 $ +s3*(-3.470000d-01)
1885 A1g = s*( 6.020000d-01)+s2*( 5.090000d-01)
1886 $ +s3*(-2.054000d+00)+s4*( 1.392000d+00)
1887 AA1g= s*(-9.220000d-02)+s2*(-1.899000d+00)
1888 $ +s3*( 4.180000d+00)+s4*(-2.494000d+00)
1889 B1g =-2.895000d-01+s*( 3.760000d-01)+s2*(-1.719000d+00)
1890 $ +s3*( 1.116000d+00)
1891 C1g = 1.439000d+00+s*(-5.570000d-01)+s2*( 3.660000d-01)
1892 $ +s3*( 7.330000d-01)+s4*(-7.620000d-01)
1893 else
1894c over 100 GeV^2
1895c
1896cccccc WHIT5 gluon (O100)
1897c
1898 A0g = 1.995000d+00+s*(-3.260000d+00)+s2*( 1.818000d+00)
1899 $ +s3*( 1.711000d+00)+s4*(-4.990000d+00)
1900 B0g =-4.660000d-01+s*(-6.100000d-01)+s2*( 1.691000d+00)
1901 $ +s3*(-6.680000d+00)+s4*( 1.019000d+01)
1902 C0g = 1.075000d+01+s*( 5.420000d+00)+s2*( 6.550000d+00)
1903 $ +s3*(-2.297000d+01)+s4*( 1.867000d+01)
1904 A1g = 2.992000d-01+s*( 1.179000d+00)+s2*(-1.915000d+00)
1905 $ +s3*( 7.260000d+00)+s4*(-1.839000d+01)
1906 AA1g=-1.600000d-01+s*(-1.114000d+00)+s2*( 2.939000d+00)
1907 $ +s3*(-6.660000d+00)+s4*( 1.923000d+01)
1908 B1g =-4.830000d-01+s*( 7.550000d-01)+s2*(-3.800000d+00)
1909 $ +s3*( 1.075000d+01)+s4*(-1.993000d+01)
1910 C1g = 1.297000d+00+s*(-1.669000d-01)+s2*( 1.906000d+00)
1911 $ +s3*(-2.057000d+00)
1912 endif
1913c
1914 WHIT5G = prsccf/alinv/x*
1915 $ ( A0g * x**B0g * x1**C0g
1916 $ +(A1g+AA1g*x) * x**B1g * x1**C1g )
1917c
1918 return
1919 end
1920c
1921cccccccccccccccccccccccccccccc
1922c QPM calculation
1923 subroutine WHIT5Q(x,r,cv,cs)
1924ccc INPUTS : x,r=mc^2/Q^2
1925ccc OUTPUTS: cv,cs (valence- and sea- charm quark dist)
1926ccc cv <-- cv / ( alpha / 2PI)
1927ccc cs <-- cs / ( alpha_s/2PI)
1928c
1929 implicit none
1930c arg
1931 DOUBLE PRECISION x,r
1932 DOUBLE PRECISION cv,cs
1933c CONST
1934 DOUBLE PRECISION ec,mc
1935 parameter(ec=2.0d0/3.0d0,mc=1.5d0)
1936c N=15 Gauss int. weights and points
1937 integer GN,i
1938 parameter(GN=15)
1939 DOUBLE PRECISION XG(GN), XW(GN)
1940 DATA (XG(i),i=1,GN)/6.003741d-03,
1941 $ 3.136330d-02, 7.589671d-02, 1.377911d-01, 2.145139d-01,
1942 $ 3.029243d-01, 3.994030d-01, 5.000000d-01, 6.005970d-01,
1943 $ 6.970757d-01, 7.854861d-01, 8.622089d-01, 9.241033d-01,
1944 $ 9.686367d-01, 9.939963d-01/
1945 DATA (XW(i),i=1,GN)/1.537662d-02,
1946 $ 3.518302d-02, 5.357961d-02, 6.978534d-02, 8.313460d-02,
1947 $ 9.308050d-02, 9.921574d-02, 1.012891d-01, 9.921574d-02,
1948 $ 9.308050d-02, 8.313460d-02, 6.978534d-02, 5.357961d-02,
1949 $ 3.518302d-02, 1.537662d-02/
1950c local
1951 DOUBLE PRECISION sum, y,z,beta,w,WHIT5G,L
1952 parameter(L=4.0d0)
1953 DOUBLE PRECISION x1,rx,z1,rz
1954c
1955c begin
1956 x1=1.0d0-x
1957 rx=4.0d0*r*x
1958c
1959c direct
1960 beta=dsqrt(1.0d0-rx/x1)
1961 w=x*( beta*(-1.0d0+8.0d0*x*x1-rx*x1)
1962 $ +(x**2+x1**2+rx*(1.0d0-3.0d0*x)-0.5d0*rx**2)
1963 $ * log( (1.0d0+beta)/(1.0d0-beta) ))
1964 cv = 3.0d0*ec**2 * w / x
1965c
1966c resolved
1967 sum=0.0d0
1968 do 10 i=1,GN
1969 y= x+rx + (x1-rx)*XG(i)**L
1970 z=x/y
1971 z1=1.0d0-z
1972 rz=4.0d0*r*z
1973 beta=dsqrt(1.0d0-rz/z1)
1974 w=z*( beta*(-1.0d0+8.0d0*z*z1-rz*z1)
1975 $ +(z**2+z1**2+rz*(1.0d0-3.0d0*z)-0.5d0*rz**2)
1976 $ * log( (1.0d0+beta)/(1.0d0-beta) ))
1977c
1978 sum= sum + w * WHIT5G(y,mc**2/r)* L*XG(i)**(L-1.0d0)*XW(i)
1979c
1980 10 continue
1981c
1982 cs = 0.5d0/x * (x1-rx) * sum
1983c
1984 return
1985 end
1986c
1987ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1988 DOUBLE PRECISION function WHIT6G(x,Q2)
1989c input: x,Q2
1990c output: clg
1991c (gluon dist.)
1992ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1993 implicit none
1994c arg
1995 DOUBLE PRECISION Q2,x
1996c const
1997 DOUBLE PRECISION q42it,q52it,lam42,lam52
1998 DOUBLE PRECISION alinv
1999c local
2000 DOUBLE PRECISION A0g,B0g,C0g,A1g,AA1g,B1g,C1g
2001 DOUBLE PRECISION s,s2,s3,s4,prsccf
2002 DOUBLE PRECISION x1
2003c parameters
2004 parameter(lam42=0.16d0, lam52=0.091411319d0)
2005 parameter(Q42IT=4.0d0, Q52IT=100.0d0)
2006 parameter(alinv=137.036d0)
2007 common /scale/ s,s2,s3,s4,prsccf
2008c
2009c begin
2010 x1=1.0d0-x
2011c
2012 if(Q2.lt.100.0d0) then
2013c under 100 GeV^2
2014c
2015cccccc WHIT6 gluon (U100)
2016c
2017 A0g = 1.600000d+01+s*(-6.100000d+01)+s2*( 1.278000d+02)
2018 $ +s3*(-1.399000d+02)+s4*( 5.990000d+01)
2019 B0g = s*(-1.109000d+00)+s2*( 8.450000d-01)
2020 $ +s3*(-3.510000d-01)
2021 C0g = 1.500000d+01+s*( 1.596000d-01)+s2*( 4.180000d+00)
2022 $ +s3*(-1.765000d-01)
2023 A1g = s*( 6.020000d-01)+s2*( 5.090000d-01)
2024 $ +s3*(-2.054000d+00)+s4*( 1.392000d+00)
2025 AA1g= s*(-9.220000d-02)+s2*(-1.899000d+00)
2026 $ +s3*( 4.180000d+00)+s4*(-2.494000d+00)
2027 B1g =-2.895000d-01+s*( 3.760000d-01)+s2*(-1.719000d+00)
2028 $ +s3*( 1.116000d+00)
2029 C1g = 1.439000d+00+s*(-5.570000d-01)+s2*( 3.660000d-01)
2030 $ +s3*( 7.330000d-01)+s4*(-7.620000d-01)
2031 else
2032c over 100 GeV^2
2033c
2034cccccc WHIT6 gluon (O100)
2035c
2036 A0g = 2.378000d+00+s*(-4.380000d+00)+s2*( 5.850000d-01)
2037 $ +s3*( 8.340000d+00)+s4*(-9.920000d+00)
2038 B0g =-4.790000d-01+s*(-6.070000d-01)+s2*( 1.458000d+00)
2039 $ +s3*(-6.030000d+00)+s4*( 9.330000d+00)
2040 C0g = 1.706000d+01+s*( 4.960000d+00)+s2*( 2.497000d+01)
2041 $ +s3*(-1.582000d+02)+s4*( 2.954000d+02)
2042 A1g = 2.992000d-01+s*( 1.179000d+00)+s2*(-1.915000d+00)
2043 $ +s3*( 7.260000d+00)+s4*(-1.839000d+01)
2044 AA1g=-1.600000d-01+s*(-1.114000d+00)+s2*( 2.939000d+00)
2045 $ +s3*(-6.660000d+00)+s4*( 1.923000d+01)
2046 B1g =-4.830000d-01+s*( 7.550000d-01)+s2*(-3.800000d+00)
2047 $ +s3*( 1.075000d+01)+s4*(-1.993000d+01)
2048 C1g = 1.297000d+00+s*(-1.669000d-01)+s2*( 1.906000d+00)
2049 $ +s3*(-2.057000d+00)
2050 endif
2051c
2052 WHIT6G = prsccf/alinv/x*
2053 $ ( A0g * x**B0g * x1**C0g
2054 $ +(A1g+AA1g*x) * x**B1g * x1**C1g )
2055c
2056 return
2057 end
2058c
2059cccccccccccccccccccccccccccccc
2060c QPM calculation
2061 subroutine WHIT6Q(x,r,cv,cs)
2062ccc INPUTS : x,r=mc^2/Q^2
2063ccc OUTPUTS: cv,cs (valence- and sea- charm quark dist)
2064ccc cv <-- cv / ( alpha / 2PI)
2065ccc cs <-- cs / ( alpha_s/2PI)
2066c
2067 implicit none
2068c arg
2069 DOUBLE PRECISION x,r
2070 DOUBLE PRECISION cv,cs
2071c CONST
2072 DOUBLE PRECISION ec,mc
2073 parameter(ec=2.0d0/3.0d0,mc=1.5d0)
2074c N=15 Gauss int. weights and points
2075 integer GN,i
2076 parameter(GN=15)
2077 DOUBLE PRECISION XG(GN), XW(GN)
2078 DATA (XG(i),i=1,GN)/6.003741d-03,
2079 $ 3.136330d-02, 7.589671d-02, 1.377911d-01, 2.145139d-01,
2080 $ 3.029243d-01, 3.994030d-01, 5.000000d-01, 6.005970d-01,
2081 $ 6.970757d-01, 7.854861d-01, 8.622089d-01, 9.241033d-01,
2082 $ 9.686367d-01, 9.939963d-01/
2083 DATA (XW(i),i=1,GN)/1.537662d-02,
2084 $ 3.518302d-02, 5.357961d-02, 6.978534d-02, 8.313460d-02,
2085 $ 9.308050d-02, 9.921574d-02, 1.012891d-01, 9.921574d-02,
2086 $ 9.308050d-02, 8.313460d-02, 6.978534d-02, 5.357961d-02,
2087 $ 3.518302d-02, 1.537662d-02/
2088c local
2089 DOUBLE PRECISION sum, y,z,beta,w,WHIT6G,L
2090 parameter(L=4.0d0)
2091 DOUBLE PRECISION x1,rx,z1,rz
2092c
2093c begin
2094 x1=1.0d0-x
2095 rx=4.0d0*r*x
2096c
2097c direct
2098 beta=dsqrt(1.0d0-rx/x1)
2099 w=x*( beta*(-1.0d0+8.0d0*x*x1-rx*x1)
2100 $ +(x**2+x1**2+rx*(1.0d0-3.0d0*x)-0.5d0*rx**2)
2101 $ * log( (1.0d0+beta)/(1.0d0-beta) ))
2102 cv = 3.0d0*ec**2 * w / x
2103c
2104c resolved
2105 sum=0.0d0
2106 do 10 i=1,GN
2107 y= x+rx + (x1-rx)*XG(i)**L
2108 z=x/y
2109 z1=1.0d0-z
2110 rz=4.0d0*r*z
2111 beta=dsqrt(1.0d0-rz/z1)
2112 w=z*( beta*(-1.0d0+8.0d0*z*z1-rz*z1)
2113 $ +(z**2+z1**2+rz*(1.0d0-3.0d0*z)-0.5d0*rz**2)
2114 $ * log( (1.0d0+beta)/(1.0d0-beta) ))
2115c
2116 sum= sum + w * WHIT6G(y,mc**2/r)* L*XG(i)**(L-1.0d0)*XW(i)
2117c
2118 10 continue
2119c
2120 cs = 0.5d0/x * (x1-rx) * sum
2121c
2122 return
2123 end