]> git.uio.no Git - u/mrichter/AliRoot.git/blame - LHAPDF/lhapdf5.3.1/wrapcteq65.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
COVERITY reports many FORWARD_NULL, corrected; AliEMCALv2: initialization of fHits...
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapcteq65.f
CommitLineData
4e9e3152 1 subroutine CTEQ65evolve(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 real*8 pdf(-6:6)
8 integer nset
9 Character Line*80
10 PARAMETER (MXX = 204, MXQ = 25, MXF = 6, MaxVal=3, nhess = 40)
11 PARAMETER (MXPQX = (MXF+1+MaxVal) * MXQ * MXX)
12 Common
13 > / CtqPar1nhess / Al, XV(0:MXX), TV(0:MXQ), UPD(0:nhess,MXPQX)
14 > / CtqPar2 / Nx, Nt, NfMx
15 > / XQrange / Qini, Qmax, Xmin
16 > / QCDtable / Alambda, Nfl, Iorder
17 > / Masstbl / Amass(6)
18
19 common/masses_LHA/cMass,bMass,tMass
20
21 data pi / 3.141592653589793d0 /
22 save
23c
24 call getnset(iset)
25 call getnmem(iset,imem)
26
27 U = X * CtLhCtq65Pdf(imem,1,X,Q)
28 D = X * CtLhCtq65Pdf(imem,2,X,Q)
29 USEA = X * CtLhCtq65Pdf(imem,-1,X,Q)
30 DSEA = X * CtLhCtq65Pdf(imem,-2,X,Q)
31 STR = X * CtLhCtq65Pdf(imem,3,X,Q)
32 CHM = X * CtLhCtq65Pdf(imem,-4,X,Q)
33 BOT = X * CtLhCtq65Pdf(imem,-5,X,Q)
34 GLU = X * CtLhCtq65Pdf(imem,0,X,Q)
35c
36 pdf(0) = glu
37 pdf(1) = d
38 pdf(-1) = dsea
39 pdf(2) = u
40 pdf(-2) = usea
41 pdf(3) = str
42 pdf(-3) = str
43 pdf(4) = chm
44 pdf(-4) = chm
45 pdf(5) = bot
46 pdf(-5) = bot
47 pdf(6) = 0.0d0
48 pdf(-6) = 0.0d0
49 return
50*
51 entry CTEQ65cevolve(x,Q,pdf) !like cteq65evolve, but allows c.ne.cbar and s.ne.sbar
52c
53 call getnset(iset)
54 call getnmem(iset,imem)
55
56 U = X * CtLhCtq65Pdf(imem,1,X,Q)
57 D = X * CtLhCtq65Pdf(imem,2,X,Q)
58 USEA = X * CtLhCtq65Pdf(imem,-1,X,Q)
59 DSEA = X * CtLhCtq65Pdf(imem,-2,X,Q)
60 STR = X * CtLhCtq65Pdf(imem,3,X,Q)
61 SBAR = X * CtLhCtq65Pdf(imem,-3,X,Q)
62 CHM = X * CtLhCtq65Pdf(imem,4,X,Q)
63 CBAR = X * CtLhCtq65Pdf(imem,-4,X,Q)
64 BOT = X * CtLhCtq65Pdf(imem,5,X,Q)
65 GLU = X * CtLhCtq65Pdf(imem,0,X,Q)
66c
67 pdf(0) = glu
68 pdf(1) = d
69 pdf(-1) = dsea
70 pdf(2) = u
71 pdf(-2) = usea
72 pdf(3) = str
73 pdf(-3) = sbar
74 pdf(4) = chm
75 pdf(-4) = cbar
76 pdf(5) = bot
77 pdf(-5) = bot
78 pdf(6) = 0.0d0
79 pdf(-6) = 0.0d0
80
81 return
82*
83 entry CTEQ65read(nset)
84 call CtLhbldat1
85 call CtLhbldat2
86
87 call LHct65set
88
89 read(1,*)nmem(nset),ndef(nset) !*** nmem+1=number of members; ndef is not used for anything ***
90 if(nmem(nset) .gt. nhess) then
91 print *,'fatal error: nmem=',nmem(nset),' > nhess=',nhess
92 stop
93 endif
94
95 MxVal = 3
96 Read (1, '(A)') Line
97 Read (1, '(A)') Line
98 Read (1, *) Dr, Fl, Al, (Amass(I),I=1,6)
99 Iorder = Nint(Dr)
100 Nfl = Nint(Fl)
101 Alambda = Al
102
103 cMass = Amass(4)
104 bMass = Amass(5)
105 tMass = Amass(6)
106
107 Read (1, '(A)') Line
108C This is the .pds (WKT) format
109 Read (1, *) N0, N0, N0, NfMx, N0, N0
110 Read (1, '(A)') Line
111 Read (1, *) NX, NT, N0, N0, N0
112 Read (1, '(A)') (Line,I=1,4)
113 Read (1, *) QINI, QMAX, (aa,TV(I), I =0, NT)
114
115 Read (1, '(A)') Line
116 Read (1, *) XMIN, aa, (XV(I), I =1, NX)
117 XV(0)=0D0
118C
119C Since quark = anti-quark for nfl>2 at this stage,
120C we Read out only the non-redundent data points
121C No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence)
122
123
124 Nblk = (NX+1) * (NT+1)
125 Npts = Nblk * (NfMx+1+MxVal)
126
127 do ihess = 0,nmem(nset) !*** new version: allows nmem < nhess ***
128
129 Read (1, '(A)') Line
130 Read (1, '(A)') Line
131 Read (1, *, IOSTAT=IRET) (UPD(ihess,I), I=1,Npts)
132
133 enddo
134 return
135*
136 entry CTEQ65cread(nset) !like CTEQ65read, but c.ne.cbar allowed
137 call CtLhbldat1
138 call CtLhbldat2
139 call LHct65set
140
141 read(1,*)nmem(nset),ndef(nset) !*** nmem+1=number of members; ndef is not used for anything ***
142 if(nmem(nset) .gt. nhess) then
143 print *,'fatal error: nmem=',nmem(nset),' > nhess=',nhess
144 stop
145 endif
146
147 MxVal = 4 !one more species free
148 Read (1, '(A)') Line
149 Read (1, '(A)') Line
150 Read (1, *) Dr, Fl, Al, (Amass(I),I=1,6)
151 Iorder = Nint(Dr)
152 Nfl = Nint(Fl)
153 Alambda = Al
154
155 cMass = Amass(4)
156 bMass = Amass(5)
157 tMass = Amass(6)
158
159 Read (1, '(A)') Line
160C This is the .pds (WKT) format
161 Read (1, *) N0, N0, N0, NfMx, N0, N0
162 Read (1, '(A)') Line
163 Read (1, *) NX, NT, N0, N0, N0
164 Read (1, '(A)') (Line,I=1,4)
165 Read (1, *) QINI, QMAX, (aa,TV(I), I =0, NT)
166
167 Read (1, '(A)') Line
168 Read (1, *) XMIN, aa, (XV(I), I =1, NX)
169 XV(0)=0D0
170C
171C we Read only the non-redundent data points
172
173 Nblk = (NX+1) * (NT+1)
174 Npts = Nblk * (NfMx+1+MxVal)
175
176 do ihess = 0,nmem(nset) !*** new version: allows nmem < nhess ***
177
178 Read (1, '(A)') Line
179 Read (1, '(A)') Line
180 Read (1, *, IOSTAT=IRET) (UPD(ihess,I), I=1,Npts)
181
182 enddo
183 return
184*
185 entry CTEQ65alfa(alfas,Qalfa)
186 alfas = pi*CtLhALPI(Qalfa)
187 return
188*
189 entry CTEQ65init(Eorder,Q2fit)
190 return
191*
192 entry CTEQ65pdf(mem)
193c imem = mem
194 call getnset(iset)
195 call setnmem(iset,mem)
196 return
197*
198 end
199
200 subroutine LHct65set
201 Implicit Double Precision (A-H,O-Z)
202 common /ctq65co/ xlast, qlast, nxsave
203 nxsave = -1000
204 xlast = -2.
205 qlast = -2.
206 return
207 end
208
209c===========================================================================
210 Function CtLhPartonX65 (iset,IPRTN, XX, QQ)
211c Given the parton distribution function in the array U in
212c COMMON / PEVLDT / , this routine interpolates to find
213c the parton distribution at an arbitray point in x and q.
214c
215 Implicit Double Precision (A-H,O-Z)
216
217 PARAMETER (MXX = 204, MXQ = 25, MXF = 6, MaxVal=3, nhess = 40)
218 PARAMETER (MXPQX = (MXF+1+MaxVal) * MXQ * MXX)
219
220 Common
221 > / CtqPar1nhess / Al, XV(0:MXX), TV(0:MXQ), UPD(0:nhess,MXPQX)
222 > / CtqPar2 / Nx, Nt, NfMx
223 > / XQrange / Qini, Qmax, Xmin
224 common /ctq65co/ xlast,qlast, nxsave
225
226 parameter(nqvec = 4)
227
228 Dimension fvec(4), fij(4)
229 Dimension xvpow(0:mxx)
230 Data OneP / 1.00001 /
231 Data xpow / 0.3d0 / !**** choice of interpolation variable
232 Save xvpow
233
234 save jq, jx, JLx, JLq, ss, sy2, sy3, s23, ty2, ty3
235 save const1 , const2, const3, const4, const5, const6
236 save tt, t13, t12, t23, t34 , t24, tmp1, tmp2, tdet
237
238c store the powers used for interpolation on first call...
239 if(nx .ne. nxsave) then
240 nxsave = nx
241 xvpow(0) = 0.D0
242 do i = 1, nx
243 xvpow(i) = xv(i)**xpow
244 enddo
245 endif
246
247 X = XX
248 Q = QQ
249
250 if((x.lt.xmin).or.(x.gt.1.d0)) print 98,x
251 98 format(' WARNING: X=',e12.5,' OUT OF RANGE')
252 if((q.lt.qini).or.(q.gt.qmax)) print 99,q
253 99 format(' WARNING: Q=',e12.5,' OUT OF RANGE')
254
255c skip the initialization in x if same as in the previous call.
256 if(x .eq. xlast) goto 100
257 xlast = x
258
259c ------------- find lower end of interval containing x, i.e.,
260c get jx such that xv(jx) .le. x .le. xv(jx+1)...
261 JLx = -1
262 JU = Nx+1
263 11 If (JU-JLx .GT. 1) Then
264 JM = (JU+JLx) / 2
265 If (X .Ge. XV(JM)) Then
266 JLx = JM
267 Else
268 JU = JM
269 Endif
270 Goto 11
271 Endif
272C Ix 0 1 2 Jx JLx Nx-2 Nx
273C |---|---|---|...|---|-x-|---|...|---|---|
274C x 0 Xmin x 1
275C
276 If (JLx .LE. -1) Then
277 Print '(A,1pE12.4)','Severe error: x <= 0 in CtLhPartonX65 x=',x
278 Stop
279 ElseIf (JLx .Eq. 0) Then
280 Jx = 0
281 Elseif (JLx .LE. Nx-2) Then
282
283C For interior points, keep x in the middle, as shown above
284 Jx = JLx - 1
285 Elseif (JLx.Eq.Nx-1 .or. x.LT.OneP) Then
286
287C We tolerate a slight over-shoot of one (OneP=1.00001),
288C perhaps due to roundoff or whatever, but not more than that.
289C Keep at least 4 points >= Jx
290 Jx = JLx - 2
291 Else
292 Print '(A,1pE12.4)','Severe error: x > 1 in CtLhPartonX65 x=',x
293 Stop
294 Endif
295C ---------- Note: JLx uniquely identifies the x-bin; Jx does not.
296
297C This is the variable to be interpolated in
298 ss = x**xpow
299
300 If (JLx.Ge.2 .and. JLx.Le.Nx-2) Then
301
302c initiation work for "interior bins": store the lattice points in s...
303 svec1 = xvpow(jx)
304 svec2 = xvpow(jx+1)
305 svec3 = xvpow(jx+2)
306 svec4 = xvpow(jx+3)
307
308 s12 = svec1 - svec2
309 s13 = svec1 - svec3
310 s23 = svec2 - svec3
311 s24 = svec2 - svec4
312 s34 = svec3 - svec4
313
314 sy2 = ss - svec2
315 sy3 = ss - svec3
316
317c constants needed for interpolating in s at fixed t lattice points...
318 const1 = s13/s23
319 const2 = s12/s23
320 const3 = s34/s23
321 const4 = s24/s23
322 s1213 = s12 + s13
323 s2434 = s24 + s34
324 sdet = s12*s34 - s1213*s2434
325 tmp = sy2*sy3/sdet
326 const5 = (s34*sy2-s2434*sy3)*tmp/s12
327 const6 = (s1213*sy2-s12*sy3)*tmp/s34
328
329 EndIf
330
331100 continue
332
333c skip the initialization in q if same as in the previous call.
334 if(q .eq. qlast) goto 110
335 qlast = q
336
337 tt = log(log(Q/Al))
338
339c --------------Now find lower end of interval containing Q, i.e.,
340c get jq such that qv(jq) .le. q .le. qv(jq+1)...
341 JLq = -1
342 JU = NT+1
343 12 If (JU-JLq .GT. 1) Then
344 JM = (JU+JLq) / 2
345 If (tt .GE. TV(JM)) Then
346 JLq = JM
347 Else
348 JU = JM
349 Endif
350 Goto 12
351 Endif
352
353 If (JLq .LE. 0) Then
354 Jq = 0
355 Elseif (JLq .LE. Nt-2) Then
356C keep q in the middle, as shown above
357 Jq = JLq - 1
358 Else
359C JLq .GE. Nt-1 case: Keep at least 4 points >= Jq.
360 Jq = Nt - 3
361
362 Endif
363C This is the interpolation variable in Q
364
365 If (JLq.GE.1 .and. JLq.LE.Nt-2) Then
366c store the lattice points in t...
367 tvec1 = Tv(jq)
368 tvec2 = Tv(jq+1)
369 tvec3 = Tv(jq+2)
370 tvec4 = Tv(jq+3)
371
372 t12 = tvec1 - tvec2
373 t13 = tvec1 - tvec3
374 t23 = tvec2 - tvec3
375 t24 = tvec2 - tvec4
376 t34 = tvec3 - tvec4
377
378 ty2 = tt - tvec2
379 ty3 = tt - tvec3
380
381 tmp1 = t12 + t13
382 tmp2 = t24 + t34
383
384 tdet = t12*t34 - tmp1*tmp2
385
386 EndIf
387
388110 continue
389
390c get the pdf function values at the lattice points...
391c In this code, we store 10 flavors: u,ubar,d,dbar,s,sbar,c,cbar,b=bbar,g
392c hence Iprtn=5 (b) is obtained from -5 (bbar)
393
394 If (Iprtn .GE. 5) Then
395 Ip = - Iprtn
396 Else
397 Ip = Iprtn
398 EndIf
399 jtmp = ((Ip + NfMx)*(NT+1)+(jq-1))*(NX+1)+jx+1
400
401 Do it = 1, nqvec
402
403 J1 = jtmp + it*(NX+1)
404
405 If (Jx .Eq. 0) Then
406C For the first 4 x points, interpolate x^2*f(x,Q)
407C This applies to the two lowest bins JLx = 0, 1
408C We cannot put the JLx.eq.1 bin into the "interior" section
409C (as we do for q), since Upd(J1) is undefined.
410 fij(1) = 0
411 fij(2) = Upd(iset,J1+1) * XV(1)**2
412 fij(3) = Upd(iset,J1+2) * XV(2)**2
413 fij(4) = Upd(iset,J1+3) * XV(3)**2
414C
415C Use CtLhPolint which allows x to be anywhere w.r.t. the grid
416
417 Call CtLhPolint4(XVpow(0), Fij(1), 4, ss, Fx, Dfx)
418
419 If (x .GT. 0D0) fvec(it) = Fx / x**2
420C Pdf is undefined for x.eq.0
421 ElseIf (JLx .Eq. Nx-1) Then
422C This is the highest x bin:
423
424c** fix allow 4 consecutive elements with iset... mrw 19.9.2005
425 fij(1) = Upd(iset,j1)
426 fij(2) = Upd(iset,j1+1)
427 fij(3) = Upd(iset,j1+2)
428 fij(4) = Upd(iset,j1+3)
429 Call CtLhPolint4 (XVpow(Nx-3), Fij(1), 4, ss, Fx, Dfx)
430
431 fvec(it) = Fx
432
433 Else
434
435C for all interior points, use Jon's in-line function
436C This applied to (JLx.Ge.2 .and. JLx.Le.Nx-2)
437c (This is cubic spline interpolation, as used by cteq; it was
438c changed to polint in previous Durham releases (jcp).)
439 sf2 = Upd(iset,J1+1)
440 sf3 = Upd(iset,J1+2)
441
442 Fvec(it) = (const5*(Upd(iset,J1)
443 & - sf2*const1 + sf3*const2)
444 & + const6*(Upd(iset,J1+3)
445 & + sf2*const3 - sf3*const4)
446 & + sf2*sy3 - sf3*sy2) / s23
447
448 Endif
449
450 enddo
451C We now have the four values Fvec(1:4)
452c interpolate in t...
453
454 If (JLq .LE. 0) Then
455C 1st Q-bin, as well as extrapolation to lower Q
456 Call CtLhPolint4(TV(0), Fvec(1), 4, tt, ff, Dfq)
457
458 ElseIf (JLq .GE. Nt-1) Then
459C Last Q-bin, as well as extrapolation to higher Q
460 Call CtLhPolint4(TV(Nt-3), Fvec(1), 4, tt, ff, Dfq)
461 Else
462C Interrior bins : (JLq.GE.1 .and. JLq.LE.Nt-2)
463C which include JLq.Eq.1 and JLq.Eq.Nt-2, since Upd is defined for
464C the full range QV(0:Nt) (in contrast to XV)
465 tf2 = fvec(2)
466 tf3 = fvec(3)
467
468 g1 = ( tf2*t13 - tf3*t12) / t23
469 g4 = (-tf2*t34 + tf3*t24) / t23
470
471 h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12
472 & + (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34)
473
474 ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23
475 EndIf
476
477 CtLhPartonX65 = ff
478
479 Return
480 End
481c===========================================================================
482 Function CtLhCtq65Pdf (iset,Iparton, X, Q)
483 Implicit Double Precision (A-H,O-Z)
484 Logical Warn
485 Common
486 > / CtqPar2 / Nx, Nt, NfMx
487 > / QCDtable / Alambda, Nfl, Iorder
488
489 Data Warn /.true./
490 save Warn
491
492 If (X .lt. 0D0 .or. X .gt. 1D0) Then
493 Print *, 'X out of range in CtLhCtq65Pdf: ', X
494 Stop
495 Endif
496 If (Q .lt. Alambda) Then
497 Print *, 'Q out of range in CtLhCtq65Pdf: ', Q
498 Stop
499 Endif
500
501c added to force pdf = 0.0 at x=1.0 exactly - mrw
502 if(x .eq. 1.0d0) then
503 CtLhCtq65Pdf = 0.0d0
504 return
505 endif
506c
507 If ((Iparton .lt. -NfMx .or. Iparton .gt. NfMx)) Then
508 If (Warn) Then
509C put a warning for calling extra flavor.
510 Warn = .false.
511 Print *, 'Warning: Iparton out of range in CtLhCtq65Pdf: '
512 > , Iparton
513 Endif
514 CtLhCtq65Pdf = 0D0
515 Return
516 Endif
517
518 CtLhCtq65Pdf = CtLhPartonX65 (iset,Iparton, X, Q)
519 if(CtLhCtq65Pdf.lt.0.D0) CtLhCtq65Pdf = 0.D0
520
521 Return
522
523C ********************
524 End