COVERITY reports many FORWARD_NULL, corrected; AliEMCALv2: initialization of fHits...
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / 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
52       call CtLhbldat2
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)               !error corrected (jcp)
98 c        imem = mem
99         call getnset(iset)
100         call setnmem(iset,mem)
101       return
102
103       end
104
105       FUNCTION CtLhPartonX5 (IPRTN, X, Q)
106 C
107 C   Given the parton distribution function in the array UPD in
108 C   COMMON / CtqPar1 / , this routine fetches u(fl, x, q) at any value of
109 C   x and q using Mth-order polynomial interpolation for x and Ln(Q/Lambda).
110 C
111       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
112 C
113       PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
114       PARAMETER (MXPQX = (MXF *2 +2) * MXQ * MXX)
115       PARAMETER (M= 2, M1 = M + 1)
116 C
117       Common 
118      > / CtqPar1 / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX)
119      > / CtqPar2 / Nx, Nt, NfMx
120      > / XQrange / Qini, Qmax, Xmin
121 C
122       Dimension Fq(M1), Df(M1)
123
124 c note, this routine doesn't have the features of not recalculating 
125 c the x or q point values if they have not changed since the last call.
126 c that makes it slower than cteq6, but since cteq5 is already obsolete,
127 c speed is not so important for it. 
128 C                                                 Work with Log (Q)
129       QG  = LOG (Q/AL)
130
131 C                           Find lower end of interval containing X
132       JL = -1
133       JU = Nx+1
134  11   If (JU-JL .GT. 1) Then
135          JM = (JU+JL) / 2
136          If (X .GT. XV(JM)) Then
137             JL = JM
138          Else
139             JU = JM
140          Endif
141          Goto 11
142       Endif
143
144       Jx = JL - (M-1)/2
145 c crude treatment if outside the defined interval.
146       If (Jx .LT. 0) then
147          Jx = 0
148       Elseif (Jx .GT. Nx-M) Then
149          Jx = Nx - M
150       Endif
151 C                                    Find the interval where Q lies
152       JL = -1
153       JU = NT+1
154  12   If (JU-JL .GT. 1) Then
155          JM = (JU+JL) / 2
156          If (QG .GT. QL(JM)) Then
157             JL = JM
158          Else
159             JU = JM
160          Endif
161          Goto 12
162       Endif
163
164       Jq = JL - (M-1)/2
165       If (Jq .LT. 0) Then
166          Jq = 0
167 c         If (Q .lt. Qini)  Print '(A, 2(1pE12.4))', 
168 c     >     ' WARNING: Q << Qini, extrapolation used; Q,Qini=',Q,Qini
169       Elseif (Jq .GT. Nt-M) Then
170          Jq = Nt - M
171 c         If (Q .gt. Qmax)  Print '(A, 2(1pE12.4))', 
172 c     >     ' WARNING: Q > Qmax, extrapolation used; Q, Qmax =', Q, Qmax
173       Endif
174
175       If (Iprtn .GE. 3) Then
176          Ip = - Iprtn
177       Else
178          Ip = Iprtn
179       EndIf
180 C                             Find the off-set in the linear array Upd
181       JFL = Ip + NfMx
182       J0  = (JFL * (NT+1) + Jq) * (NX+1) + Jx
183 C
184 C                                           Now interpolate in x for M1 Q's
185       Do 21 Iq = 1, M1
186          J1 = J0 + (Nx+1)*(Iq-1) + 1
187          Call CtLhPolint3 (XV(Jx), Upd(J1), M1, X, Fq(Iq), Df(Iq))
188  21   Continue
189 C                                          Finish off by interpolating in Q
190       Call CtLhPolint3 (QL(Jq), Fq(1), M1, QG, Ftmp, Ddf)
191
192       CtLhPartonX5 = Ftmp
193 C
194       RETURN
195
196 C                        ****************************
197       END
198
199       Function CtLhCtq5Pdf (Iparton, X, Q)
200       Implicit Double Precision (A-H,O-Z)
201       Logical Warn
202       Common
203      > / CtqPar2 / Nx, Nt, NfMx
204      > / QCDtable /  Alambda, Nfl, Iorder
205
206       Data Warn /.true./
207       save Warn
208
209       If (X .lt. 0D0 .or. X .gt. 1D0) Then
210         Print *, 'X out of range in CtLhCtq5Pdf: ', X
211         Stop
212       Endif
213       If (Q .lt. Alambda) Then
214         Print *, 'Q out of range in CtLhCtq5Pdf: ', Q
215         Stop
216       Endif
217       If ((Iparton .lt. -NfMx .or. Iparton .gt. NfMx)) Then
218          If (Warn) Then
219 C        put a warning for calling extra flavor.
220              Warn = .false.
221              Print *, 'Warning: Iparton out of range in CtLhCtq5Pdf: '
222      >              , Iparton
223          Endif
224          CtLhCtq5Pdf = 0D0
225          Return
226       Endif
227
228       CtLhCtq5Pdf = CtLhPartonX5 (Iparton, X, Q)
229       if(CtLhCtq5Pdf.lt.0.D0)  CtLhCtq5Pdf = 0.D0
230
231       Return
232
233 C                             ********************
234       End