Removing obsolete constants.
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.2.2 / wrapcteq6.f
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
21 c
22       call getnset(iset)
23       call getnmem(iset,imem)
24 c
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
35 c      
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
85 C
86 C                  Since quark = anti-quark for nfl>2 at this stage,
87 C                  we Read  out only the non-redundent data points
88 C     No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence)
89
90       Nblk = (NX+1) * (NT+1)
91       Npts =  Nblk  * (NfMx+3)
92
93
94 c **  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)
113 c        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)
123 C                                        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)
136 11    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
147 12      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
155 13    CONTINUE
156       RETURN
157       END
158
159 c===========================================================================
160       Function CtLhPartonX6 (iset,IPRTN, XX, QQ)
161 c  Given the parton distribution function in the array U in
162 c  COMMON / PEVLDT / , this routine interpolates to find
163 c  the parton distribution at an arbitray point in x and q.
164 c
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
189 c 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
201 c skip the initialization in x if same as in the previous call.
202         if(x .eq. xlast) goto 100
203         xlast = x
204
205 c      -------------    find lower end of interval containing x, i.e.,
206 c                       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
218 C                     Ix    0   1   2      Jx  JLx         Nx-2     Nx
219 C                           |---|---|---|...|---|-x-|---|...|---|---|
220 C                     x     0  Xmin               x                 1
221 C
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
229 C                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
233 C                  We tolerate a slight over-shoot of one (OneP=1.00001),
234 C              perhaps due to roundoff or whatever, but not more than that.
235 C                                      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
241 C          ---------- Note: JLx uniquely identifies the x-bin; Jx does not.
242
243 C                       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
248 c     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
263 c 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
277 100     continue
278
279 c 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
285 c         --------------Now find lower end of interval containing Q, i.e.,
286 c                          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
302 C                                  keep q in the middle, as shown above
303          Jq = JLq - 1
304       Else
305 C                         JLq .GE. Nt-1 case:  Keep at least 4 points >= Jq.
306         Jq = Nt - 3
307
308       Endif
309 C                                   This is the interpolation variable in Q
310
311       If (JLq.GE.1 .and. JLq.LE.Nt-2) Then
312 c                                        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
334 110     continue
335
336 c 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
350 C                          For the first 4 x points, interpolate x^2*f(x,Q)
351 C                           This applies to the two lowest bins JLx = 0, 1
352 C            We cannot put the JLx.eq.1 bin into the "interior" section
353 C                           (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
358 C
359 C                 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
364 C                                              Pdf is undefined for x.eq.0
365        ElseIf  (JLx .Eq. Nx-1) Then
366 C                                                This is the highest x bin:
367
368 c** 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
379 C for all interior points, use Jon's in-line function
380 C This applied to (JLx.Ge.2 .and. JLx.Le.Nx-2)
381 c (This is cubic spline interpolation, as used by cteq; it was 
382 c 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
395 C                                   We now have the four values Fvec(1:4)
396 c     interpolate in t...
397
398       If (JLq .LE. 0) Then
399 C                         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
403 C                         Last Q-bin, as well as extrapolation to higher Q
404         Call CtLhPolint4(TV(Nt-3), Fvec(1), 4, tt, ff, Dfq)
405       Else
406 C                         Interrior bins : (JLq.GE.1 .and. JLq.LE.Nt-2)
407 C       which include JLq.Eq.1 and JLq.Eq.Nt-2, since Upd is defined for
408 C                         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
424 C                                       ********************
425       End
426 c===========================================================================
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
446 c 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
451 c
452       If ((Iparton .lt. -NfMx .or. Iparton .gt. NfMx)) Then
453          If (Warn) Then
454 C        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
468 C                             ********************
469       End
470       SUBROUTINE CtLhPOLINT4 (XA,YA,N,X,Y,DY)
471 c fast version of polint, valid only for N=4
472 c 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)
595 c fast version of polint, valid only for N=3
596 c 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