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