! -*- F90 -*- subroutine CT12evolve(x,Q,pdf) implicit real*8(a-h,o-z) include 'parmsetup.inc' character*16 name(nmxset) integer nmem(nmxset),ndef(nmxset),mmem common/NAME/name,nmem,ndef,mmem double precision gridx(nmxgridx),gridq(nmxgridq) integer ngridx,ngridq,jx,jq real*8 pdf(-6:6) integer nset Character Line*80 PARAMETER (MXX = 204, MXQ = 40, MXF = 6, MaxVal=4, nhess = 300) !JG PARAMETER (MXPQX = (MXF+1+MaxVal) * MXQ * MXX) Common & & / CTPar1nhess12 / & & Al(nmxset), QV(0:MXX,nmxset), XV(0:MXX,nmxset),TV(0:MXQ,nmxset), & & UPD(0:nhess,MXPQX,nmxset), AlsCTEQ(0:MXQ,nmxset) & & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset), Nvg(nmxset) & & / XQrange / Qini(nmxset), Qmax(nmxset), Xmin(nmxset) & & / QCDtable / Alambda(nmxset), Nfl(nmxset), Iorder(nmxset) & & / Masstbl / Amass(6,nmxset) common/masses_LHA/cMass(nmxset),bMass(nmxset),tMass(nmxset) data pi / 3.141592653589793d0 / save ! call getnset(iset) call getnmem(iset,imem) U = X * CtLhCT12Pdf(imem,1,X,Q) D = X * CtLhCT12Pdf(imem,2,X,Q) USEA = X * CtLhCT12Pdf(imem,-1,X,Q) DSEA = X * CtLhCT12Pdf(imem,-2,X,Q) ASTR = X * CtLhCT12Pdf(imem,-3,X,Q) STR = X * CtLhCT12Pdf(imem,3,X,Q) !JG for meta CHM = X * CtLhCT12Pdf(imem,-4,X,Q) BOT = X * CtLhCT12Pdf(imem,-5,X,Q) GLU = X * CtLhCT12Pdf(imem,0,X,Q) ! pdf(0) = glu pdf(1) = d pdf(-1) = dsea pdf(2) = u pdf(-2) = usea pdf(3) = str !JG for meta pdf(-3) = astr pdf(4) = chm pdf(-4) = chm pdf(5) = bot pdf(-5) = bot pdf(6) = 0.0d0 pdf(-6) = 0.0d0 return ! entry CT12getgrid(nset,ngridx,ngridq,gridx,gridq) do jx=0,nx(nset) gridx(jx+1)=xv(jx,nset) enddo do jq=0,nt(nset) gridq(jq+1)=qv(jq,nset)*qv(jq,nset) enddo ngridx=nx(nset) ngridq=nt(nset) return ! entry CT12read(nset) call CtLhbldat1 call CtLhbldat2 call LHCT12set !*** nmem+1=number of member read(1,*)nmem(nset),ndef(nset) if(nmem(nset) .gt. nhess) then print *,'fatal error: nmem=',nmem(nset),' > nhess=',nhess stop endif MxVal = 2 Read (1, '(A)') Line Read (1, '(A)') Line if(Line(1:11) .eq. ' ipk, Ordr') then ipdsformat = 10 Read(1, *) ipk, Dr, QQalfa, alfaQ, (Amass(I,nset),I=1,6) Iorder(nset) = Nint(Dr) Read(1, '(A)') Line if(Line(1:7) .eq. ' IMASS' ) then ipdsformat = 11 read(1, *) aimass, fswitch, N0, N0, N0, Nfmx(nset), MxVal else Read (1, *) N0, N0, N0, NfMx(nset), MxVal endif else ipdsformat = 6 Read (1, *) Dr, Fl, Al(nset), (Amass(I,nset),I=1,6) Iorder(nset) = Nint(Dr) Nfmx(nset) = Nint(Fl) Alambda(nset) = Al(nset) Read (1, '(A)') Line Read (1, *) N0, N0, N0, NfMx(nset), N0, N0 endif cMass(nset) = Amass(4,nset) bMass(nset) = Amass(5,nset) tMass(nset) = Amass(6,nset) Read (1, '(A)') Line Read (1, *) NX(nset), NT(nset), N0, NG, N0 if(NG.gt.0) Read(1,'(A)') (line, i=1,NG+1) Read (1, '(A)') Line if(ipdsformat.eq.11) then Read (1, *) & & QINI(nset), QMAX(nset), (QV(I,nset),TV(I,nset),AlsCTEQ(I,nset), I =0, NT(nset)) qbase1=QV(1,nset)/Exp(Exp(Tv(1,nset))) qbase2=QV(NT(nset),nset)/Exp(Exp(Tv(NT(nset),nset))) if(abs(qbase1-qbase2).gt.1e-5) then print *,'something wrong with qbase' print *,'qbase1, qbase2=',qbase1,qbase2 stop else Al(nset)=(qbase1+qbase2)/2.0d0 Alambda(nset) = Al(nset) endif else Read (1, *) & & QINI(nset), QMAX(nset), (aa,TV(I,nset), I =0, NT(nset)) endif Read (1, '(A)') Line Read (1, *) XMIN(nset), aa, (XV(I,nset), I =1, NX(nset)) XV(0,nset)=0D0 Nvg(nset)=Mxval !JG for meta ! ! Since quark = anti-quark for nfl>2 at this stage, ! we Read out only the non-redundent data points ! No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence) Nblk = (NX(nset)+1) * (NT(nset)+1) Npts = Nblk * (NfMx(nset)+1+MxVal) !*** new version: allows nm do ihess = 0,nmem(nset) Read (1, '(A)') Line Read (1, '(A)') Line Read (1, *, IOSTAT=IRET) (UPD(ihess,I,nset), I=1,Npts) enddo return ! ! entry CT12alfa(alfas,Qalfa) if(ipdsformat.eq.11) then alfas = CtLhCT12Alphas(Qalfa) else alfas = pi*CtLhALPI(Qalfa) endif return ! entry CT12init(Eorder,Q2fit) return ! entry CT12pdf(mem) ! imem = mem call getnset(iset) call setnmem(iset,mem) return ! END subroutine LHCT12set Implicit Double Precision (A-H,O-Z) common /CT12co/ xlast, qlast, nxsave nxsave = -1000 xlast = -2. qlast = -2. return END !======================================================================= Function CtLhPartonX12 (imem,IPRTN, XX, QQ) ! Given the parton distribution function in the array U in ! COMMON / PEVLDT / , this routine interpolates to find ! the parton distribution at an arbitray point in x and q. ! Implicit Double Precision (A-H,O-Z) include 'parmsetup.inc' PARAMETER (MXX = 204, MXQ = 40, MXF = 6, MaxVal=4, nhess = 300) PARAMETER (MXPQX = (MXF+1+MaxVal) * MXQ * MXX) Common & & / CTPar1nhess12 / & & Al(nmxset), QV(0:MXX,nmxset), XV(0:MXX,nmxset), TV(0:MXQ,nmxset), & & UPD(0:nhess,MXPQX,nmxset), AlsCTEQ(0:MXQ,nmxset) & & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset), Nvg(nmxset) & & / XQrange / Qini(nmxset), Qmax(nmxset), Xmin(nmxset) common /CT12co/ xlast,qlast, nxsave parameter(nqvec = 4) Dimension fvec(4), fij(4) Dimension xvpow(0:mxx) Data OneP / 1.00001 / !**** choice of interpolation variable Data xpow / 0.3d0 / Save xvpow Data ixprint,iqprint/0,0/ save ixprint,iqprint save jq, jx, JLx, JLq, ss, sy2, sy3, s23, ty2, ty3 save const1 , const2, const3, const4, const5, const6 save tt, t13, t12, t23, t34 , t24, tmp1, tmp2, tdet call getnset(iset) call getnmem(iset,imem) ! store the powers used for interpolation on first call... if(nx(iset) .ne. nxsave) then nxsave = nx(iset) xvpow(0) = 0.D0 do i = 1, nx(iset) xvpow(i) = xv(i,iset)**xpow enddo endif X = XX Q = QQ if((x.lt.xmin(iset)).or.(x.gt.1.d0)) then ixprint=ixprint+1 if(ixprint.lt.11) print 98,x if(ixprint.eq.10) print *, & & 'more warning messages like the last suppressed.' endif 98 format(' WARNING: X=',e12.5,' OUT OF RANGE') if((q.lt.qini(iset)).or.(q.gt.qmax(iset))) then iqprint=iqprint+1 if(iqprint.lt.11) print 99,q if(iqprint.eq.10) print *, & & 'more warning messages like the last suppressed.' endif 99 format(' WARNING: Q=',e12.5,' OUT OF RANGE') ! skip the initialization in x if same as in the previous call. if(x .eq. xlast) goto 100 xlast = x ! ------------- find lower end of interval containing x, i.e., ! get jx such that xv(jx) .le. x .le. xv(jx+1)... JLx = -1 JU = Nx(iset)+1 11 If (JU-JLx .GT. 1) Then JM = (JU+JLx) / 2 If (X .Ge. XV(JM,iset)) Then JLx = JM Else JU = JM Endif Goto 11 Endif ! Ix 0 1 2 Jx JLx Nx-2 Nx ! |---|---|---|...|---|-x-|---|...|---|---| ! x 0 Xmin x 1 ! If (JLx .LE. -1) Then Print '(A,1pE12.4)','Severe error: x <= 0 in CtLhPartonX12 x=',x Stop ElseIf (JLx .Eq. 0) Then Jx = 0 Elseif (JLx .LE. Nx(iset)-2) Then ! For interior points, keep x in the middle, as shown abo Jx = JLx - 1 Elseif (JLx.Eq.Nx(iset)-1 .or. x.LT.OneP) Then ! We tolerate a slight over-shoot of one (OneP=1.00001) ! perhaps due to roundoff or whatever, but not more than th ! Keep at least 4 points >= Jx Jx = JLx - 2 Else Print '(A,1pE12.4)','Severe error: x > 1 in CtLhPartonX12 x=',x Stop Endif ! ---------- Note: JLx uniquely identifies the x-bin; Jx does n ! This is the variable to be interpolated in ss = x**xpow If (JLx.Ge.2 .and. JLx.Le.Nx(iset)-2) Then ! initiation work for "interior bins": store the lattice points in s svec1 = xvpow(jx) svec2 = xvpow(jx+1) svec3 = xvpow(jx+2) svec4 = xvpow(jx+3) s12 = svec1 - svec2 s13 = svec1 - svec3 s23 = svec2 - svec3 s24 = svec2 - svec4 s34 = svec3 - svec4 sy2 = ss - svec2 sy3 = ss - svec3 ! constants needed for interpolating in s at fixed t lattice points... const1 = s13/s23 const2 = s12/s23 const3 = s34/s23 const4 = s24/s23 s1213 = s12 + s13 s2434 = s24 + s34 sdet = s12*s34 - s1213*s2434 tmp = sy2*sy3/sdet const5 = (s34*sy2-s2434*sy3)*tmp/s12 const6 = (s1213*sy2-s12*sy3)*tmp/s34 EndIf 100 continue ! skip the initialization in q if same as in the previous call. if(q .eq. qlast) goto 110 qlast = q tt = log(log(Q/Al(iset))) ! --------------Now find lower end of interval containing Q, i.e ! get jq such that qv(jq) .le. q .le. qv(jq+1). JLq = -1 JU = NT(iset)+1 12 If (JU-JLq .GT. 1) Then JM = (JU+JLq) / 2 If (tt .GE. TV(JM,iset)) Then JLq = JM Else JU = JM Endif Goto 12 Endif If (JLq .LE. 0) Then Jq = 0 Elseif (JLq .LE. Nt(iset)-2) Then ! keep q in the middle, as shown above Jq = JLq - 1 Else ! JLq .GE. Nt-1 case: Keep at least 4 points >= Jq = Nt(iset) - 3 Endif ! This is the interpolation variable i If (JLq.GE.1 .and. JLq.LE.Nt(iset)-2) Then ! store the lattice points in t.. tvec1 = Tv(jq,iset) tvec2 = Tv(jq+1,iset) tvec3 = Tv(jq+2,iset) tvec4 = Tv(jq+3,iset) t12 = tvec1 - tvec2 t13 = tvec1 - tvec3 t23 = tvec2 - tvec3 t24 = tvec2 - tvec4 t34 = tvec3 - tvec4 ty2 = tt - tvec2 ty3 = tt - tvec3 tmp1 = t12 + t13 tmp2 = t24 + t34 tdet = t12*t34 - tmp1*tmp2 EndIf 110 continue ! get the pdf function values at the lattice points... ! In this code, we store 10 flavors: u,ubar,d,dbar,s,sbar,c,cbar,b=bbar, ! hence Iprtn=5 (b) is obtained from -5 (bbar) If (Iprtn .GE. 5) Then Ip = - Iprtn Else Ip = Iprtn EndIf jtmp = ((Ip + NfMx(iset))*(NT(iset)+1)+(jq-1))*(NX(iset)+1)+jx+1 Do it = 1, nqvec J1 = jtmp + it*(NX(iset)+1) If (Jx .Eq. 0) Then ! For the first 4 x points, interpolate x^2*f(x ! This applies to the two lowest bins JLx = 0, ! We cannot put the JLx.eq.1 bin into the "interior" section ! (as we do for q), since Upd(J1) is undefined fij(1) = 0 fij(2) = Upd(imem,J1+1,iset) * XV(1,iset)**2 fij(3) = Upd(imem,J1+2,iset) * XV(2,iset)**2 fij(4) = Upd(imem,J1+3,iset) * XV(3,iset)**2 ! ! Use CtLhPolint which allows x to be anywhere w.r.t. th Call CtLhPolint4(XVpow(0), Fij(1), 4, ss, Fx, Dfx) If (x .GT. 0D0) fvec(it) = Fx / x**2 ! Pdf is undefined for x.eq ElseIf (JLx .Eq. Nx(iset)-1) Then ! This is the highest x b !** fix allow 4 consecutive elements with iset... mrw 19.9.2005 fij(1) = Upd(imem,j1,iset) fij(2) = Upd(imem,j1+1,iset) fij(3) = Upd(imem,j1+2,iset) fij(4) = Upd(imem,j1+3,iset) Call CtLhPolint4 (XVpow(Nx(iset)-3), Fij(1), 4, ss, Fx, Dfx) fvec(it) = Fx Else ! for all interior points, use Jon's in-line function ! This applied to (JLx.Ge.2 .and. JLx.Le.Nx-2) ! (This is cubic spline interpolation, as used by cteq; it was ! changed to polint in previous Durham releases (jcp).) sf2 = Upd(imem,J1+1,iset) sf3 = Upd(imem,J1+2,iset) Fvec(it) = (const5*(Upd(imem,J1,iset) & & - sf2*const1 + sf3*const2) & & + const6*(Upd(imem,J1+3,iset) & & + sf2*const3 - sf3*const4) & & + sf2*sy3 - sf3*sy2) / s23 Endif enddo ! We now have the four values Fvec(1:4 ! interpolate in t... If (JLq .LE. 0) Then ! 1st Q-bin, as well as extrapolation to lower Q Call CtLhPolint4(TV(0,iset), Fvec(1), 4, tt, ff, Dfq) ElseIf (JLq .GE. Nt(iset)-1) Then ! Last Q-bin, as well as extrapolation to higher Call CtLhPolint4(TV(Nt(iset)-3,iset), Fvec(1), 4, tt, ff, Dfq) Else ! Interrior bins : (JLq.GE.1 .and. JLq.LE.Nt-2) ! which include JLq.Eq.1 and JLq.Eq.Nt-2, since Upd is defined for ! the full range QV(0:Nt) (in contrast to XV) tf2 = fvec(2) tf3 = fvec(3) g1 = ( tf2*t13 - tf3*t12) / t23 g4 = (-tf2*t34 + tf3*t24) / t23 h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12 & & + (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34) ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23 EndIf CtLhPartonX12 = ff Return END !======================================================================= Function CtLhCt12Pdf (imem,Iparton, X, Q) Implicit Double Precision (A-H,O-Z) include 'parmsetup.inc' Logical Warn Common & & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset), Nvg(nmxset) & & / QCDtable / Alambda(nmxset), Nfl(nmxset), Iorder(nmxset) Data Warn /.true./ save Warn call getnset(iset) If (X .lt. 0D0 .or. X .gt. 1D0) Then Print *, 'X out of range in CtLhCt12Pdf: ', X Stop Endif If (Q .lt. Alambda(iset)) Then Print *, 'Q out of range in CtLhCt12Pdf: ', Q Stop Endif ! added to force pdf = 0.0 at x=1.0 exactly - mrw if(x .eq. 1.0d0) then CtLhCt12Pdf = 0.0d0 return endif ! If ((Iparton .lt. -NfMx(iset) .or. Iparton .gt. NfMx(iset))) Then If (Warn) Then ! put a warning for calling extra flavor. Warn = .false. Print *, 'Warning: Iparton out of range in CtLhCt12Pdf: ' & & , Iparton Endif CtLhCt12Pdf = 0D0 Return Endif If ((Iparton.eq.3).and.(Nvg(iset).eq.2)) then CtLhCt12Pdf = CtLhPartonX12 (imem, -3, X, Q) ! JG for meta else CtLhCt12Pdf = CtLhPartonX12 (imem,Iparton, X, Q) Endif if(CtLhCt12Pdf.lt.0.D0) CtLhCt12Pdf = 0.D0 Return ! ******************** !=========================================================================== END Function CtLhCT12Alphas (QQ) Implicit Double Precision (A-H,O-Z) include 'parmsetup.inc' PARAMETER (MXX = 204, MXQ = 40, MXF = 6, MaxVal=4, nhess = 300) PARAMETER (MXPQX = (MXF+1+MaxVal) * MXQ * MXX) double precision Alsout Common & & / CTPar1nhess12 / & & Al(nmxset), QV(0:MXX,nmxset), XV(0:MXX,nmxset),TV(0:MXQ,nmxset), & & UPD(0:nhess,MXPQX,nmxset), AlsCTEQ(0:MXQ,nmxset) & & / CtqPar2 / Nx(nmxset), Nt(nmxset), NfMx(nmxset), Nvg(nmxset) & & / XQrange / Qini(nmxset), Qmax(nmxset), Xmin(nmxset) Data Q, JQ /-1D0, 0/ save call getnset(iset) call getnmem(iset,imem) Q = QQ tt = log(log(Q/Al(iset))) ! -------------- Find lower end of interval containing Q, i.e., ! get jq such that qv(jq) .le. q .le. qv(jq+1)... JLq = -1 JU = NT(iset)+1 13 If (JU-JLq .GT. 1) Then JM = (JU+JLq) / 2 If (tt .GE. TV(JM,iset)) Then JLq = JM Else JU = JM Endif Goto 13 Endif If (JLq .LE. 0) Then Jq = 0 Elseif (JLq .LE. Nt(iset)-2) Then ! keep q in the middle, as shown above Jq = JLq - 1 Else ! JLq .GE. Nt-1 case: Keep at least 4 points >= Jq. Jq = Nt(iset) - 3 Endif ! This is the interpolation variable in Q Call CtLhPolint4 (TV(jq,iset), AlsCTEQ(jq,iset),4,tt,Alsout,Dfq) CtLhCT12Alphas = Alsout Return ! ******************** End !======================================================================