4 subroutine CTEQ65evolve(x,Q,pdf)
5 implicit real*8(a-h,o-z)
6 include 'parmsetup.inc'
7 character*16 name(nmxset)
9 integer nmem(nmxset),ndef(nmxset),mmem
10 common/NAME/name,nmem,ndef,mmem
11 double precision gridx(nmxgridx),gridq(nmxgridq)
12 integer ngridx,ngridq,jx,jq
16 PARAMETER (MXX = 204, MXQ = 25, MXF = 6, MaxVal=3, nhess = 0)
17 PARAMETER (MXPQX = (MXF+1+MaxVal) * MXQ * MXX)
19 & / CtqPar1nhess65 / &
20 & Al(nmxset), XV(0:MXX,nmxset), TV(0:MXQ,nmxset), &
21 & UPD(0:nhess,MXPQX,nmxset) &
22 & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset) &
23 & / XQrange / Qini(nmxset), Qmax(nmxset), Xmin(nmxset) &
24 & / QCDtable / Alambda(nmxset), Nfl(nmxset), Iorder(nmxset) &
25 & / Masstbl / Amass(6,nmxset)
27 common/masses_LHA/cMass(nmxset),bMass(nmxset),tMass(nmxset)
29 data pi / 3.141592653589793d0 /
33 call getnmem(iset,imem)
35 U = X * CtLhCtq65Pdf(0,1,X,Q)
36 D = X * CtLhCtq65Pdf(0,2,X,Q)
37 USEA = X * CtLhCtq65Pdf(0,-1,X,Q)
38 DSEA = X * CtLhCtq65Pdf(0,-2,X,Q)
39 STR = X * CtLhCtq65Pdf(0,-3,X,Q)
40 CHM = X * CtLhCtq65Pdf(0,-4,X,Q)
41 BOT = X * CtLhCtq65Pdf(0,-5,X,Q)
42 GLU = X * CtLhCtq65Pdf(0,0,X,Q)
59 !like cteq65evolve, but allows
60 entry CTEQ65cevolve(x,Q,pdf)
63 call getnmem(iset,imem)
65 U = X * CtLhCtq65Pdf(0,1,X,Q)
66 D = X * CtLhCtq65Pdf(0,2,X,Q)
67 USEA = X * CtLhCtq65Pdf(0,-1,X,Q)
68 DSEA = X * CtLhCtq65Pdf(0,-2,X,Q)
69 STR = X * CtLhCtq65Pdf(0,3,X,Q)
70 SBAR = X * CtLhCtq65Pdf(0,-3,X,Q)
71 CHM = X * CtLhCtq65Pdf(0,4,X,Q)
72 CBAR = X * CtLhCtq65Pdf(0,-4,X,Q)
73 BOT = X * CtLhCtq65Pdf(0,5,X,Q)
74 GLU = X * CtLhCtq65Pdf(0,0,X,Q)
92 !like cteq65evolve, but allows
93 entry CTEQ65sevolve(x,Q,pdf)
96 call getnmem(iset,imem)
98 U = X * CtLhCtq65Pdf(0,1,X,Q)
99 D = X * CtLhCtq65Pdf(0,2,X,Q)
100 USEA = X * CtLhCtq65Pdf(0,-1,X,Q)
101 DSEA = X * CtLhCtq65Pdf(0,-2,X,Q)
102 STR = X * CtLhCtq65Pdf(0,3,X,Q)
103 SBAR = X * CtLhCtq65Pdf(0,-3,X,Q)
104 CHM = X * CtLhCtq65Pdf(0,-4,X,Q)
105 BOT = X * CtLhCtq65Pdf(0,5,X,Q)
106 GLU = X * CtLhCtq65Pdf(0,0,X,Q)
124 entry CTEQ65getgrid(nset,ngridx,ngridq,gridx,gridq)
126 gridx(jx+1)=xv(jx,nset)
129 qtemp = alambda(nset)*dexp(dexp(tv(jq,nset)))
130 gridq(jq+1)=qtemp*qtemp
136 entry CTEQ65read(nset)
142 !*** nmem+1=number of member
143 read(1,*)nmem(nset),ndef(nset)
144 ! if(nmem(nset) .gt. nhess) then
145 ! print *,'fatal error: nmem=',nmem(nset),' > nhess=',nhess
152 Read (1, *) Dr, Fl, Al(nset), (Amass(I,nset),I=1,6)
153 Iorder(nset) = Nint(Dr)
155 Alambda(nset) = Al(nset)
157 cMass(nset) = Amass(4,nset)
158 bMass(nset) = Amass(5,nset)
159 tMass(nset) = Amass(6,nset)
162 ! This is the .pds (WKT) f
163 Read (1, *) N0, N0, N0, NfMx(nset), N0, N0
165 Read (1, *) NX(nset), NT(nset), N0, N0, N0
166 Read (1, '(A)') (Line,I=1,4)
168 & QINI(nset), QMAX(nset), (aa,TV(I,nset), I =0, NT(nset))
171 Read (1, *) XMIN(nset), aa, (XV(I,nset), I =1, NX(nset))
174 ! Since quark = anti-quark for nfl>2 at this stage,
175 ! we Read out only the non-redundent data points
176 ! No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence)
179 Nblk = (NX(nset)+1) * (NT(nset)+1)
180 Npts = Nblk * (NfMx(nset)+1+MxVal)
182 !*** new version: allows nm
183 do ihess = 0,nmem(nset)
187 Read (1, *, IOSTAT=IRET) (UPD(0,I,nset), I=1,Npts)
192 entry CTEQ66read(nset)
198 !*** nmem+1=number of member
199 read(1,*)nmem(nset),ndef(nset)
200 ! if(nmem(nset) .gt. nhess) then
201 ! print *,'fatal error: nmem=',nmem(nset),' > nhess=',nhess
208 Read (1, *) Dr, Fl, Al(nset), (Amass(I,nset),I=1,6)
209 Iorder(nset) = Nint(Dr)
211 Alambda(nset) = Al(nset)
213 cMass(nset) = Amass(4,nset)
214 bMass(nset) = Amass(5,nset)
215 tMass(nset) = Amass(6,nset)
218 ! This is the .pds (WKT) f
219 Read (1, *) N0, N0, N0, NfMx(nset), N0, N0
221 Read (1, *) NX(nset), NT(nset), N0, N0, N0
222 Read (1, '(A)') (Line,I=1,4)
224 & QINI(nset), QMAX(nset), (aa,TV(I,nset), I =0, NT(nset))
227 Read (1, *) XMIN(nset), aa, (XV(I,nset), I =1, NX(nset))
230 ! Since quark = anti-quark for nfl>2 at this stage,
231 ! we Read out only the non-redundent data points
232 ! No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence)
235 Nblk = (NX(nset)+1) * (NT(nset)+1)
236 Npts = Nblk * (NfMx(nset)+1+MxVal)
238 !*** new version: allows nm
239 do ihess = 0,nmem(nset)
243 Read (1, *, IOSTAT=IRET) (UPD(0,I,nset), I=1,Npts)
248 !like CTEQ65read, but c.ne.
249 entry CTEQ65cread(nset)
254 !*** nmem+1=number of member
255 read(1,*)nmem(nset),ndef(nset)
256 ! if(nmem(nset) .gt. nhess) then
257 ! print *,'fatal error: nmem=',nmem(nset),' > nhess=',nhess
261 !one more species free
265 Read (1, *) Dr, Fl, Al(nset), (Amass(I,nset),I=1,6)
266 Iorder(nset) = Nint(Dr)
268 Alambda(nset) = Al(nset)
270 cMass(nset) = Amass(4,nset)
271 bMass(nset) = Amass(5,nset)
272 tMass(nset) = Amass(6,nset)
275 ! This is the .pds (WKT) f
276 Read (1, *) N0, N0, N0, NfMx(nset), N0, N0
278 Read (1, *) NX(nset), NT(nset), N0, N0, N0
279 Read (1, '(A)') (Line,I=1,4)
281 & QINI(nset), QMAX(nset), (aa,TV(I,nset), I =0, NT(nset))
284 Read (1, *) XMIN(nset), aa, (XV(I,nset), I =1, NX(nset))
287 ! we Read only the non-redundent data points
289 Nblk = (NX(nset)+1) * (NT(nset)+1)
290 Npts = Nblk * (NfMx(nset)+1+MxVal)
292 !*** new version: allows nm
293 do ihess = 0,nmem(nset)
297 Read (1, *, IOSTAT=IRET) (UPD(0,I,nset), I=1,Npts)
302 entry CTEQ65alfa(alfas,Qalfa)
303 alfas = pi*CtLhALPI(Qalfa)
306 entry CTEQ65init(Eorder,Q2fit)
312 call setnmem(iset,mem)
314 ! have to reopen stream 1
315 call getsetpath(setpath)
316 open(1,file=setpath(1:len_trim(setpath)),action='READ')
319 do while (line(1:3).ne.'==>')
322 ! - backspace by one record
324 ! - dummy read up to the member requested
328 Read (1, *, IOSTAT=IRET) (UPD(0,I,iset), I=1,Npts)
330 !- read in the data of the requested member
333 Read (1, *, IOSTAT=IRET) (UPD(0,I,iset), I=1,Npts)
340 Implicit Double Precision (A-H,O-Z)
341 common /ctq65co/ xlast, qlast, nxsave
348 !=======================================================================
349 Function CtLhPartonX65 (imem,IPRTN, XX, QQ)
350 ! Given the parton distribution function in the array U in
351 ! COMMON / PEVLDT / , this routine interpolates to find
352 ! the parton distribution at an arbitray point in x and q.
354 Implicit Double Precision (A-H,O-Z)
356 include 'parmsetup.inc'
357 PARAMETER (MXX = 204, MXQ = 25, MXF = 6, MaxVal=3, nhess = 0)
358 PARAMETER (MXPQX = (MXF+1+MaxVal) * MXQ * MXX)
361 & / CtqPar1nhess65 / &
362 & Al(nmxset), XV(0:MXX,nmxset), TV(0:MXQ,nmxset), &
363 & UPD(0:nhess,MXPQX,nmxset) &
364 & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset) &
365 & / XQrange / Qini(nmxset), Qmax(nmxset), Xmin(nmxset)
366 common /ctq65co/ xlast,qlast, nxsave
370 Dimension fvec(4), fij(4)
371 Dimension xvpow(0:mxx)
372 Data OneP / 1.00001 /
373 !**** choice of interpolation variable
376 Data ixprint,iqprint/0,0/
379 save jq, jx, JLx, JLq, ss, sy2, sy3, s23, ty2, ty3
380 save const1 , const2, const3, const4, const5, const6
381 save tt, t13, t12, t23, t34 , t24, tmp1, tmp2, tdet
385 ! call getnmem(iset,imem)
387 ! store the powers used for interpolation on first call...
388 if(nx(iset) .ne. nxsave) then
392 xvpow(i) = xv(i,iset)**xpow
399 if((x.lt.xmin(iset)).or.(x.gt.1.d0)) then
401 if(ixprint.lt.11) print 98,x
402 if(ixprint.eq.10) print *, &
403 & 'more warning messages like the last suppressed.'
405 98 format(' WARNING: X=',e12.5,' OUT OF RANGE')
406 if((q.lt.qini(iset)).or.(q.gt.qmax(iset))) then
408 if(iqprint.lt.11) print 99,q
409 if(iqprint.eq.10) print *, &
410 & 'more warning messages like the last suppressed.'
412 99 format(' WARNING: Q=',e12.5,' OUT OF RANGE')
414 ! skip the initialization in x if same as in the previous call.
415 if(x .eq. xlast) goto 100
418 ! ------------- find lower end of interval containing x, i.e.,
419 ! get jx such that xv(jx) .le. x .le. xv(jx+1)...
422 11 If (JU-JLx .GT. 1) Then
424 If (X .Ge. XV(JM,iset)) Then
431 ! Ix 0 1 2 Jx JLx Nx-2 Nx
432 ! |---|---|---|...|---|-x-|---|...|---|---|
435 If (JLx .LE. -1) Then
436 Print '(A,1pE12.4)','Severe error: x <= 0 in CtLhPartonX65 x=',x
438 ElseIf (JLx .Eq. 0) Then
440 Elseif (JLx .LE. Nx(iset)-2) Then
442 ! For interior points, keep x in the middle, as shown abo
444 Elseif (JLx.Eq.Nx(iset)-1 .or. x.LT.OneP) Then
446 ! We tolerate a slight over-shoot of one (OneP=1.00001)
447 ! perhaps due to roundoff or whatever, but not more than th
448 ! Keep at least 4 points >= Jx
451 Print '(A,1pE12.4)','Severe error: x > 1 in CtLhPartonX65 x=',x
454 ! ---------- Note: JLx uniquely identifies the x-bin; Jx does n
456 ! This is the variable to be interpolated in
459 If (JLx.Ge.2 .and. JLx.Le.Nx(iset)-2) Then
461 ! initiation work for "interior bins": store the lattice points in s
476 ! constants needed for interpolating in s at fixed t lattice points...
483 sdet = s12*s34 - s1213*s2434
485 const5 = (s34*sy2-s2434*sy3)*tmp/s12
486 const6 = (s1213*sy2-s12*sy3)*tmp/s34
492 ! skip the initialization in q if same as in the previous call.
493 if(q .eq. qlast) goto 110
496 tt = log(log(Q/Al(iset)))
498 ! --------------Now find lower end of interval containing Q, i.e
499 ! get jq such that qv(jq) .le. q .le. qv(jq+1).
502 12 If (JU-JLq .GT. 1) Then
504 If (tt .GE. TV(JM,iset)) Then
514 Elseif (JLq .LE. Nt(iset)-2) Then
515 ! keep q in the middle, as shown above
518 ! JLq .GE. Nt-1 case: Keep at least 4 points >=
522 ! This is the interpolation variable i
524 If (JLq.GE.1 .and. JLq.LE.Nt(iset)-2) Then
525 ! store the lattice points in t..
527 tvec2 = Tv(jq+1,iset)
528 tvec3 = Tv(jq+2,iset)
529 tvec4 = Tv(jq+3,iset)
543 tdet = t12*t34 - tmp1*tmp2
549 ! get the pdf function values at the lattice points...
550 ! In this code, we store 10 flavors: u,ubar,d,dbar,s,sbar,c,cbar,b=bbar,
551 ! hence Iprtn=5 (b) is obtained from -5 (bbar)
553 If (Iprtn .GE. 5) Then
558 jtmp = ((Ip + NfMx(iset))*(NT(iset)+1)+(jq-1))*(NX(iset)+1)+jx+1
562 J1 = jtmp + it*(NX(iset)+1)
565 ! For the first 4 x points, interpolate x^2*f(x
566 ! This applies to the two lowest bins JLx = 0,
567 ! We cannot put the JLx.eq.1 bin into the "interior" section
568 ! (as we do for q), since Upd(J1) is undefined
570 fij(2) = Upd(imem,J1+1,iset) * XV(1,iset)**2
571 fij(3) = Upd(imem,J1+2,iset) * XV(2,iset)**2
572 fij(4) = Upd(imem,J1+3,iset) * XV(3,iset)**2
574 ! Use CtLhPolint which allows x to be anywhere w.r.t. th
576 Call CtLhPolint4(XVpow(0), Fij(1), 4, ss, Fx, Dfx)
578 If (x .GT. 0D0) fvec(it) = Fx / x**2
579 ! Pdf is undefined for x.eq
580 ElseIf (JLx .Eq. Nx(iset)-1) Then
581 ! This is the highest x b
583 !** fix allow 4 consecutive elements with iset... mrw 19.9.2005
584 fij(1) = Upd(imem,j1,iset)
585 fij(2) = Upd(imem,j1+1,iset)
586 fij(3) = Upd(imem,j1+2,iset)
587 fij(4) = Upd(imem,j1+3,iset)
588 Call CtLhPolint4 (XVpow(Nx(iset)-3), Fij(1), 4, ss, Fx, Dfx)
594 ! for all interior points, use Jon's in-line function
595 ! This applied to (JLx.Ge.2 .and. JLx.Le.Nx-2)
596 ! (This is cubic spline interpolation, as used by cteq; it was
597 ! changed to polint in previous Durham releases (jcp).)
598 sf2 = Upd(imem,J1+1,iset)
599 sf3 = Upd(imem,J1+2,iset)
601 Fvec(it) = (const5*(Upd(imem,J1,iset) &
602 & - sf2*const1 + sf3*const2) &
603 & + const6*(Upd(imem,J1+3,iset) &
604 & + sf2*const3 - sf3*const4) &
605 & + sf2*sy3 - sf3*sy2) / s23
610 ! We now have the four values Fvec(1:4
611 ! interpolate in t...
614 ! 1st Q-bin, as well as extrapolation to lower Q
615 Call CtLhPolint4(TV(0,iset), Fvec(1), 4, tt, ff, Dfq)
617 ElseIf (JLq .GE. Nt(iset)-1) Then
618 ! Last Q-bin, as well as extrapolation to higher
619 Call CtLhPolint4(TV(Nt(iset)-3,iset), Fvec(1), 4, tt, ff, Dfq)
621 ! Interrior bins : (JLq.GE.1 .and. JLq.LE.Nt-2)
622 ! which include JLq.Eq.1 and JLq.Eq.Nt-2, since Upd is defined for
623 ! the full range QV(0:Nt) (in contrast to XV)
627 g1 = ( tf2*t13 - tf3*t12) / t23
628 g4 = (-tf2*t34 + tf3*t24) / t23
630 h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12 &
631 & + (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34)
633 ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23
640 !=======================================================================
641 Function CtLhCtq65Pdf (imem,Iparton, X, Q)
642 Implicit Double Precision (A-H,O-Z)
643 include 'parmsetup.inc'
646 & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset) &
647 & / QCDtable / Alambda(nmxset), Nfl(nmxset), Iorder(nmxset)
653 If (X .lt. 0D0 .or. X .gt. 1D0) Then
654 Print *, 'X out of range in CtLhCtq65Pdf: ', X
657 If (Q .lt. Alambda(iset)) Then
658 Print *, 'Q out of range in CtLhCtq65Pdf: ', Q
662 ! added to force pdf = 0.0 at x=1.0 exactly - mrw
663 if(x .eq. 1.0d0) then
668 If ((Iparton .lt. -NfMx(iset) .or. Iparton .gt. NfMx(iset))) Then
670 ! put a warning for calling extra flavor.
672 Print *, 'Warning: Iparton out of range in CtLhCtq65Pdf: ' &
679 CtLhCtq65Pdf = CtLhPartonX65 (imem,Iparton, X, Q)
680 if(CtLhCtq65Pdf.lt.0.D0) CtLhCtq65Pdf = 0.D0
684 ! ********************