1
C============================================================================
2
C CTEQ Parton Distribution Functions: Version 6
3
C January 24, 2002, v6.0
6
C Ref: "New Generation of Parton Distributions with
7
C Uncertainties from Global QCD Analysis"
8
C By: J. Pumplin, D.R. Stump, J.Huston, H.L. Lai, P. Nadolsky, W.K. Tung
11
C This package contains 3 standard sets of CTEQ6 PDF's and 40 up/down sets
12
C with respect to CTEQ6M PDF's. Details are:
13
C ---------------------------------------------------------------------------
14
C Iset PDF Description Alpha_s(Mz)**Lam4 Lam5 Table_File
15
C ---------------------------------------------------------------------------
16
C 1 CTEQ6M Standard MSbar scheme 0.118 326 226 cteq6m.tbl
17
C 2 CTEQ6D Standard DIS scheme 0.118 326 226 cteq6d.tbl
18
C 3 CTEQ6L Leading Order 0.118** 326** 226 cteq6l.tbl
19
C 4 CTEQ6L1 Leading Order 0.130 215 165 cteq6l1.tbl
20
C ------------------------------
21
C 1xx CTEQ6M1xx +/- w.r.t. CTEQ6M 0.118 326 226 cteq6m1xx.tbl
23
C ---------------------------------------------------------------------------
24
C ** ALL fits are obtained by using the same coupling strength
25
C \alpha_s(Mz)=0.118 and the NLO running \alpha_s formula, except CTEQ6L1
26
C which uses the LO running \alpha_s and its value determined from the fit.
27
C For the LO fits, the evolution of the PDF and the hard cross sections are
28
C calculated at LO. More detailed discussions are given in hep-ph/0201195.
30
C The table grids are generated for 10^-6 < x < 1 and 1.3 < Q < 10,000 (GeV).
31
C PDF values outside of the above range are returned using extrapolation.
32
C Lam5 (Lam4) represents Lambda value (in MeV) for 5 (4) flavors.
33
C The matching alpha_s between 4 and 5 flavors takes place at Q=4.5 GeV,
34
C which is defined as the bottom quark mass, whenever it can be applied.
36
C The Table_Files are assumed to be in the working directory.
38
C Before using the PDF, it is necessary to do the initialization by
40
C where Iset is the desired PDF specified in the above table.
42
C The function Ctq6Pdf (Iparton, X, Q)
43
C returns the parton distribution inside the proton for parton [Iparton]
44
C at [X] Bjorken_X and scale [Q] (GeV) in PDF set [Iset].
45
C Iparton is the parton label (5, 4, 3, 2, 1, 0, -1, ......, -5)
46
C for (b, c, s, d, u, g, u_bar, ..., b_bar),
48
C For detailed information on the parameters used, e.q. quark masses,
49
C QCD Lambda, ... etc., see info lines at the beginning of the
52
C These programs, as provided, are in double precision. By removing the
53
C "Implicit Double Precision" lines, they can also be run in single
56
C If you have detailed questions concerning these CTEQ6 distributions,
57
C or if you find problems/bugs using this package, direct inquires to
58
C Pumplin@pa.msu.edu or Tung@pa.msu.edu.
60
C===========================================================================
62
Function Ctq6Pdf (Iparton, X, Q)
63
Implicit Double Precision (A-H,O-Z)
66
> / CtqPar2 / Nx, Nt, NfMx
67
> / QCDtable / Alambda, Nfl, Iorder
72
If (X .lt. 0D0 .or. X .gt. 1D0) Then
73
Print *, 'X out of range in Ctq6Pdf: ', X
76
If (Q .lt. Alambda) Then
77
Print *, 'Q out of range in Ctq6Pdf: ', Q
80
If ((Iparton .lt. -NfMx .or. Iparton .gt. NfMx)) Then
82
C put a warning for calling extra flavor.
84
Print *, 'Warning: Iparton out of range in Ctq6Pdf: '
91
Ctq6Pdf = PartonX6 (Iparton, X, Q)
92
if(Ctq6Pdf.lt.0.D0) Ctq6Pdf = 0.D0
96
C ********************
99
Subroutine SetCtq6 (Iset)
100
Implicit Double Precision (A-H,O-Z)
101
Parameter (Isetmax0=4)
102
Character Flnm(Isetmax0)*6, nn*3, Tablefile*40
103
Data (Flnm(I), I=1,Isetmax0)
104
> / 'cteq6m', 'cteq6d', 'cteq6l', 'cteq6l'/
105
Data Isetold, Isetmin0, Isetmin1, Isetmax1 /-987,1,101,140/
108
C If data file not initialized, do so.
109
If(Iset.ne.Isetold) then
111
If (Iset.ge.Isetmin0 .and. Iset.le.3) Then
112
Tablefile=Flnm(Iset)//'.tbl'
113
Elseif (Iset.eq.Isetmax0) Then
114
Tablefile=Flnm(Iset)//'1.tbl'
115
Elseif (Iset.ge.Isetmin1 .and. Iset.le.Isetmax1) Then
116
write(nn,'(I3)') Iset
117
Tablefile=Flnm(1)//nn//'.tbl'
119
Print *, 'Invalid Iset number in SetCtq6 :', Iset
122
c Open(IU, File='Pdfdata/'//Tablefile, Status='OLD', Err=100)
123
call OpenData(TableFile)
124
21 Call ReadTbl6 (IU)
130
100 Print *, ' Data file ', Tablefile, ' cannot be opened '
133
C ********************
136
Subroutine ReadTbl6 (Nu)
137
Implicit Double Precision (A-H,O-Z)
139
PARAMETER (MXX = 96, MXQ = 20, MXF = 5)
140
PARAMETER (MXPQX = (MXF + 3) * MXQ * MXX)
142
> / CtqPar1 / Al, XV(0:MXX), TV(0:MXQ), UPD(MXPQX)
143
> / CtqPar2 / Nx, Nt, NfMx
144
> / XQrange / Qini, Qmax, Xmin
145
> / QCDtable / Alambda, Nfl, Iorder
146
> / Masstbl / Amass(6)
148
Read (Nu, '(A)') Line
149
Read (Nu, '(A)') Line
150
Read (Nu, *) Dr, Fl, Al, (Amass(I),I=1,6)
155
Read (Nu, '(A)') Line
156
Read (Nu, *) NX, NT, NfMx
158
Read (Nu, '(A)') Line
159
Read (Nu, *) QINI, QMAX, (TV(I), I =0, NT)
161
Read (Nu, '(A)') Line
162
Read (Nu, *) XMIN, (XV(I), I =0, NX)
165
TV(Iq) = Log(Log (TV(Iq) /Al))
168
C Since quark = anti-quark for nfl>2 at this stage,
169
C we Read out only the non-redundent data points
170
C No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence)
172
Nblk = (NX+1) * (NT+1)
173
Npts = Nblk * (NfMx+3)
174
Read (Nu, '(A)') Line
175
Read (Nu, *, IOSTAT=IRET) (UPD(I), I=1,Npts)
178
C ****************************
182
C Returns an unallocated FORTRAN i/o unit.
186
INQUIRE (UNIT=N, OPENED=EX)
192
Stop ' There is no available I/O unit. '
193
C *************************
197
SUBROUTINE POLINT6 (XA,YA,N,X,Y,DY)
199
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
200
C Adapted from "Numerical Recipes"
202
DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
207
IF (DIFT.LT.DIF) THEN
238
Function PartonX6 (IPRTN, XX, QQ)
240
c Given the parton distribution function in the array U in
241
c COMMON / PEVLDT / , this routine interpolates to find
242
c the parton distribution at an arbitray point in x and q.
244
Implicit Double Precision (A-H,O-Z)
246
Parameter (MXX = 96, MXQ = 20, MXF = 5)
247
Parameter (MXQX= MXQ * MXX, MXPQX = MXQX * (MXF+3))
250
> / CtqPar1 / Al, XV(0:MXX), TV(0:MXQ), UPD(MXPQX)
251
> / CtqPar2 / Nx, Nt, NfMx
252
> / XQrange / Qini, Qmax, Xmin
254
Dimension fvec(4), fij(4)
255
Dimension xvpow(0:mxx)
256
Data OneP / 1.00001 /
257
Data xpow / 0.3d0 / !**** choice of interpolation variable
262
c store the powers used for interpolation on first call...
263
if(ientry .eq. 0) then
268
xvpow(i) = xv(i)**xpow
276
c ------------- find lower end of interval containing x, i.e.,
277
c get jx such that xv(jx) .le. x .le. xv(jx+1)...
280
11 If (JU-JLx .GT. 1) Then
282
If (X .Ge. XV(JM)) Then
289
C Ix 0 1 2 Jx JLx Nx-2 Nx
290
C |---|---|---|...|---|-x-|---|...|---|---|
293
If (JLx .LE. -1) Then
294
Print '(A,1pE12.4)', 'Severe error: x <= 0 in PartonX6! x = ', x
296
ElseIf (JLx .Eq. 0) Then
298
Elseif (JLx .LE. Nx-2) Then
300
C For interrior points, keep x in the middle, as shown above
302
Elseif (JLx.Eq.Nx-1 .or. x.LT.OneP) Then
304
C We tolerate a slight over-shoot of one (OneP=1.00001),
305
C perhaps due to roundoff or whatever, but not more than that.
306
C Keep at least 4 points >= Jx
309
Print '(A,1pE12.4)', 'Severe error: x > 1 in PartonX6! x = ', x
312
C ---------- Note: JLx uniquely identifies the x-bin; Jx does not.
314
C This is the variable to be interpolated in
317
If (JLx.Ge.2 .and. JLx.Le.Nx-2) Then
319
c initiation work for "interior bins": store the lattice points in s...
334
c constants needed for interpolating in s at fixed t lattice points...
341
sdet = s12*s34 - s1213*s2434
343
const5 = (s34*sy2-s2434*sy3)*tmp/s12
344
const6 = (s1213*sy2-s12*sy3)*tmp/s34
348
c --------------Now find lower end of interval containing Q, i.e.,
349
c get jq such that qv(jq) .le. q .le. qv(jq+1)...
352
12 If (JU-JLq .GT. 1) Then
354
If (tt .GE. TV(JM)) Then
364
Elseif (JLq .LE. Nt-2) Then
365
C keep q in the middle, as shown above
368
C JLq .GE. Nt-1 case: Keep at least 4 points >= Jq.
372
C This is the interpolation variable in Q
374
If (JLq.GE.1 .and. JLq.LE.Nt-2) Then
375
c store the lattice points in t...
393
tdet = t12*t34 - tmp1*tmp2
398
c get the pdf function values at the lattice points...
400
If (Iprtn .GE. 3) Then
405
jtmp = ((Ip + NfMx)*(NT+1)+(jq-1))*(NX+1)+jx+1
409
J1 = jtmp + it*(NX+1)
412
C For the first 4 x points, interpolate x^2*f(x,Q)
413
C This applies to the two lowest bins JLx = 0, 1
414
C We can not put the JLx.eq.1 bin into the "interrior" section
415
C (as we do for q), since Upd(J1) is undefined.
417
fij(2) = Upd(J1+1) * XV(1)**2
418
fij(3) = Upd(J1+2) * XV(2)**2
419
fij(4) = Upd(J1+3) * XV(3)**2
421
C Use Polint6 which allows x to be anywhere w.r.t. the grid
423
Call Polint6 (XVpow(0), Fij(1), 4, ss, Fx, Dfx)
425
If (x .GT. 0D0) Fvec(it) = Fx / x**2
426
C Pdf is undefined for x.eq.0
427
ElseIf (JLx .Eq. Nx-1) Then
428
C This is the highest x bin:
430
Call Polint6 (XVpow(Nx-3), Upd(J1), 4, ss, Fx, Dfx)
435
C for all interior points, use Jon's in-line function
436
C This applied to (JLx.Ge.2 .and. JLx.Le.Nx-2)
440
g1 = sf2*const1 - sf3*const2
441
g4 = -sf2*const3 + sf3*const4
443
Fvec(it) = (const5*(Upd(J1)-g1)
444
& + const6*(Upd(J1+3)-g4)
445
& + sf2*sy3 - sf3*sy2) / s23
450
C We now have the four values Fvec(1:4)
451
c interpolate in t...
454
C 1st Q-bin, as well as extrapolation to lower Q
455
Call Polint6 (TV(0), Fvec(1), 4, tt, ff, Dfq)
457
ElseIf (JLq .GE. Nt-1) Then
458
C Last Q-bin, as well as extrapolation to higher Q
459
Call Polint6 (TV(Nt-3), Fvec(1), 4, tt, ff, Dfq)
461
C Interrior bins : (JLq.GE.1 .and. JLq.LE.Nt-2)
462
C which include JLq.Eq.1 and JLq.Eq.Nt-2, since Upd is defined for
463
C the full range QV(0:Nt) (in contrast to XV)
467
g1 = ( tf2*t13 - tf3*t12) / t23
468
g4 = (-tf2*t34 + tf3*t24) / t23
470
h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12
471
& + (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34)
473
ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23
479
C ********************
added
removed
Template/Source/PDF/Ctq6Pdf.f
