LHAPDF 5.2.2 source code.
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.2.2 / wrapcteq5.f
1       subroutine CTEQ5evolve(x,Q,pdf)
2       implicit real*8(a-h,o-z)
3       include 'parmsetup.inc'
4       character*16 name(nmxset)
5       integer nmem(nmxset),ndef(nmxset),mmem
6       common/NAME/name,nmem,ndef,mmem
7       integer nset
8       real*8 pdf(-6:6)
9       Character Line*80
10       PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
11       PARAMETER (MXPQX = (MXF *2 +2) * MXQ * MXX)
12       PARAMETER (M= 2, M1 = M + 1)
13       Common
14      > / CtqPar1 / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX)
15      > / CtqPar2 / Nx, Nt, NfMx
16      > / XQrange / Qini, Qmax, Xmin
17      > / QCDtable /  Alambda, Nfl, Iorder
18      > / Masstbl / Amass(6)
19       data pi / 3.141592653589793d0 /
20       save
21 c
22       U =         X * CtLhCtq5Pdf(1,X,Q)
23       D =         X * CtLhCtq5Pdf(2,X,Q)
24       USEA =      X * CtLhCtq5Pdf(-1,X,Q)
25       DSEA =      X * CtLhCtq5Pdf(-2,X,Q)
26       STR =       X * CtLhCtq5Pdf(3,X,Q)
27       CHM =       X * CtLhCtq5Pdf(4,X,Q)
28       BOT =       X * CtLhCtq5Pdf(5,X,Q)
29       GLU  =      X * CtLhCtq5Pdf(0,X,Q)
30       UPV=U-USEA
31       DNV=D-DSEA
32 c      
33       pdf(0)  = glu
34       pdf(1)  = dnv+dsea
35       pdf(-1) = dsea
36       pdf(2)  = upv+usea
37       pdf(-2) = usea
38       pdf(3)  = str
39       pdf(-3) = str
40       pdf(4)  = chm
41       pdf(-4) = chm
42       pdf(5)  = bot
43       pdf(-5) = bot
44       pdf(6)  = 0.0d0
45       pdf(-6) = 0.0d0
46
47       return
48 *
49       entry CTEQ5read(nset)
50
51       call CtLhbldat1           !this line was missing in previous releases (jcp)
52       call CtLhbldat2           !this line was missing in previous releases (jcp)
53  
54       read(1,*)nmem(nset),ndef(nset)
55       Read  (1, '(A)') Line     
56       Read  (1, '(A)') Line
57       Read  (1, *) Dr, Fl, Al, (Amass(I),I=1,6)
58       Iorder = Nint(Dr)
59       Nfl = Nint(Fl)
60       Alambda = Al
61
62       Read  (1, '(A)') Line 
63       Read  (1, *) NX,  NT, NfMx
64
65       Read  (1, '(A)') Line
66       Read  (1, *) QINI, QMAX, (QL(I), I =0, NT)
67
68       Read  (1, '(A)') Line
69       Read  (1, *) XMIN, (XV(I), I =0, NX)
70
71       Do 11 Iq = 0, NT
72          QL(Iq) = Log (QL(Iq) /Al)
73    11 Continue
74 C
75 C                  Since quark = anti-quark for nfl>2 at this stage, 
76 C                  we Read  out only the non-redundent data points
77 C     No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence) 
78
79       Nblk = (NX+1) * (NT+1)
80       Npts =  Nblk  * (NfMx+3)
81  
82       Read  (1, '(A)') Line
83       Read  (1, *, IOSTAT=IRET) (UPD(I), I=1,Npts)
84       Return
85 *    
86
87       entry CTEQ5alfa(alfas,Qalfa)
88
89        alfas = pi*CtLhALPI(Qalfa)
90
91       return
92 *
93       entry CTEQ5init(Eorder,Q2fit)
94
95       return
96 *
97       entry CTEQ5pdf(mem)
98         imem = mem
99       return
100
101       end
102
103       FUNCTION CtLhPartonX5 (IPRTN, X, Q)
104 C
105 C   Given the parton distribution function in the array UPD in
106 C   COMMON / CtqPar1 / , this routine fetches u(fl, x, q) at any value of
107 C   x and q using Mth-order polynomial interpolation for x and Ln(Q/Lambda).
108 C
109       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
110 C
111       PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
112       PARAMETER (MXPQX = (MXF *2 +2) * MXQ * MXX)
113       PARAMETER (M= 2, M1 = M + 1)
114 C
115       Logical First
116       Common 
117      > / CtqPar1 / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX)
118      > / CtqPar2 / Nx, Nt, NfMx
119      > / XQrange / Qini, Qmax, Xmin
120 C
121       Dimension Fq(M1), Df(M1)
122
123       Data First /.true./
124       save First
125 C                                                 Work with Log (Q)
126       QG  = LOG (Q/AL)
127
128 C                           Find lower end of interval containing X
129       JL = -1
130       JU = Nx+1
131  11   If (JU-JL .GT. 1) Then
132          JM = (JU+JL) / 2
133          If (X .GT. XV(JM)) Then
134             JL = JM
135          Else
136             JU = JM
137          Endif
138          Goto 11
139       Endif
140
141       Jx = JL - (M-1)/2
142       If (X .lt. Xmin .and. First ) Then
143          First = .false.
144 c         Print '(A, 2(1pE12.4))', 
145 c     >     ' WARNING: X << Xmin, extrapolation used; X,Xmin=',X,Xmin
146          If (Jx .LT. 0) Jx = 0
147       Elseif (Jx .GT. Nx-M) Then
148          Jx = Nx - M
149       Endif
150 C                                    Find the interval where Q lies
151       JL = -1
152       JU = NT+1
153  12   If (JU-JL .GT. 1) Then
154          JM = (JU+JL) / 2
155          If (QG .GT. QL(JM)) Then
156             JL = JM
157          Else
158             JU = JM
159          Endif
160          Goto 12
161       Endif
162
163       Jq = JL - (M-1)/2
164       If (Jq .LT. 0) Then
165          Jq = 0
166 c         If (Q .lt. Qini)  Print '(A, 2(1pE12.4))', 
167 c     >     ' WARNING: Q << Qini, extrapolation used; Q,Qini=',Q,Qini
168       Elseif (Jq .GT. Nt-M) Then
169          Jq = Nt - M
170 c         If (Q .gt. Qmax)  Print '(A, 2(1pE12.4))', 
171 c     >     ' WARNING: Q > Qmax, extrapolation used; Q, Qmax =', Q, Qmax
172       Endif
173
174       If (Iprtn .GE. 3) Then
175          Ip = - Iprtn
176       Else
177          Ip = Iprtn
178       EndIf
179 C                             Find the off-set in the linear array Upd
180       JFL = Ip + NfMx
181       J0  = (JFL * (NT+1) + Jq) * (NX+1) + Jx
182 C
183 C                                           Now interpolate in x for M1 Q's
184       Do 21 Iq = 1, M1
185          J1 = J0 + (Nx+1)*(Iq-1) + 1
186          Call CtLhPolint3 (XV(Jx), Upd(J1), M1, X, Fq(Iq), Df(Iq))
187  21   Continue
188 C                                          Finish off by interpolating in Q
189       Call CtLhPolint3 (QL(Jq), Fq(1), M1, QG, Ftmp, Ddf)
190
191       CtLhPartonX5 = Ftmp
192 C
193       RETURN
194 C                        ****************************
195       END
196
197       Function CtLhCtq5Pdf (Iparton, X, Q)
198       Implicit Double Precision (A-H,O-Z)
199       Logical Warn
200       Common
201      > / CtqPar2 / Nx, Nt, NfMx
202      > / QCDtable /  Alambda, Nfl, Iorder
203
204       Data Warn /.true./
205       save Warn
206
207       If (X .lt. 0D0 .or. X .gt. 1D0) Then
208         Print *, 'X out of range in CtLhCtq5Pdf: ', X
209         Stop
210       Endif
211       If (Q .lt. Alambda) Then
212         Print *, 'Q out of range in CtLhCtq5Pdf: ', Q
213         Stop
214       Endif
215       If ((Iparton .lt. -NfMx .or. Iparton .gt. NfMx)) Then
216          If (Warn) Then
217 C        put a warning for calling extra flavor.
218              Warn = .false.
219              Print *, 'Warning: Iparton out of range in CtLhCtq5Pdf: '
220      >              , Iparton
221          Endif
222          CtLhCtq5Pdf = 0D0
223          Return
224       Endif
225
226       CtLhCtq5Pdf = CtLhPartonX5 (Iparton, X, Q)
227       if(CtLhCtq5Pdf.lt.0.D0)  CtLhCtq5Pdf = 0.D0
228
229       Return
230
231 C                             ********************
232       End