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