PROGRAM MAKECTEQ6 * This is an interface to the smoothened CTEQ6 parton distribution * routine that produces a table in the standard form. * D. E. Soper and P. Anandam. 17 May 1999. * D. E. Soper 26 November 2002. DOUBLE PRECISION X,LX,LXMIN,LXMAX,DELTALX DOUBLE PRECISION MU,LMU,LMUMIN,LMUMAX,DELTALMU INTEGER ISET DOUBLE PRECISION Ctq6Pdf CHARACTER*100 HEADER DOUBLE PRECISION LAMBDA(3) INTEGER LAMBDAFLAVORS(3) DOUBLE PRECISION XMIN,XMAX INTEGER NX DOUBLE PRECISION MUMIN,MUMAX INTEGER NMU CHARACTER*20 FILENAME(3) CHARACTER*20 NAME(3) DOUBLE PRECISION PTN(-2:6,32,32) DOUBLE PRECISION A, B, AP, BP, CP, TG, ANGL, PT * Loop variables INTEGER I, N,M DATA FILENAME / 'cteq6m.ptn', 'cteq6d.ptn', $ 'cteq6l.ptn' / DATA NAME /'CTEQ6M','CTEQ6D','CTEQ6L'/ DATA LAMBDA /0.226D0,0.226D0,0.226D0/ DATA LAMBDAFLAVORS /5,5,5/ XMIN = 0.1000000000D-05 XMAX = 0.9500000000D0 NX = 32 MUMIN = 1.3D0 MUMAX = 10000D0 - 1D-6 NMU = 32 DO ISET=1,3 CALL SetCtq6(ISET) WRITE(*,*) 'Making ',FILENAME(ISET) OPEN(UNIT = 99, FILE = FILENAME(ISET), STATUS = 'UNKNOWN') HEADER = '! ' // NAME(ISET)(1:7) // $ ' distribution generated by DES 26 Nov 2002' WRITE(99,*)HEADER WRITE(99,*)'!' WRITE(99,*)'PARTON DATA' WRITE(99,*)' VERSION = 110496' WRITE(99,103)LAMBDA(ISET) 103 FORMAT(' LAMBDA',' = ',G20.10) WRITE(99,104)LAMBDAFLAVORS(ISET) 104 FORMAT('NFL_LAMBDA',' = ',I10) WRITE(99,105)XMIN 105 FORMAT(' XMIN',' = ',G20.10) WRITE(99,106)XMAX 106 FORMAT(' XMAX',' = ',G20.10) WRITE(99,107)NX 107 FORMAT(' N_XPOINTS',' = ',I10) WRITE(99,108)MUMIN 108 FORMAT(' MUMIN',' = ',G20.10) WRITE(99,109)MUMAX 109 FORMAT(' MUMAX',' = ',G20.10) WRITE(99,110)NMU 110 FORMAT('N_MUPOINTS',' = ',I10) * Calculate lattice info. LXMIN = DLOG(XMIN/(1.0D0 - XMIN)) LXMAX = DLOG(XMAX/(1.0D0 - XMAX)) DELTALX = (LXMAX - LXMIN)/(NX-1) LMUMIN = DLOG(MUMIN) LMUMAX = DLOG(MUMAX) DELTALMU = (LMUMAX - LMUMIN)/(NMU-1) LX = LXMIN - DELTALX DO N=1,NX LX = LX + DELTALX X = DEXP(LX)/(1.0D0+DEXP(LX)) LMU = LMUMIN - DELTALMU DO M=1,NMU LMU = LMU + DELTALMU MU = DEXP(LMU) PTN(-2,N,M) = Ctq6Pdf(-2, X, MU) PTN(-1,N,M) = Ctq6Pdf(-1, X, MU) PTN(0,N,M) = Ctq6Pdf(0, X, MU) PTN(1,N,M) = Ctq6Pdf(1, X, MU) PTN(2,N,M) = Ctq6Pdf(2, X, MU) PTN(3,N,M) = Ctq6Pdf(3, X, MU) PTN(4,N,M) = Ctq6Pdf(4, X, MU) PTN(5,N,M) = Ctq6Pdf(5, X, MU) C PTN(6,N,M) = Ctq6Pdf(6, X, MU) * We check for negative values and set those values to zero. IF (PTN(-2,N,M).LT.0D0) PTN(-2,N,M) = 0D0 IF (PTN(-1,N,M).LT.0D0) PTN(-1,N,M) = 0D0 IF (PTN(0,N,M).LT.0D0) PTN(0,N,M) = 0D0 IF (PTN(1,N,M).LT.0D0) PTN(1,N,M) = 0D0 IF (PTN(2,N,M).LT.0D0) PTN(2,N,M) = 0D0 IF (PTN(3,N,M).LT.0D0) PTN(3,N,M) = 0D0 IF (PTN(4,N,M).LT.0D0) PTN(4,N,M) = 0D0 IF (PTN(5,N,M).LT.0D0) PTN(5,N,M) = 0D0 PTN(-2,N,M) = DLOG(1.0D-16 + PTN(-2,N,M)) PTN(-1,N,M) = DLOG(1.0D-16 + PTN(-1,N,M)) PTN(0,N,M) = DLOG(1.0D-16 + PTN(0,N,M)) PTN(1,N,M) = DLOG(1.0D-16 + PTN(1,N,M)) PTN(2,N,M) = DLOG(1.0D-16 + PTN(2,N,M)) PTN(3,N,M) = DLOG(1.0D-16 + PTN(3,N,M)) PTN(4,N,M) = DLOG(1.0D-16 + PTN(4,N,M)) PTN(5,N,M) = DLOG(1.0D-16 + PTN(5,N,M)) PTN(6,N,M) = DLOG(1.0D-16 + 0.0D0) ENDDO ENDDO * Smoothen out the distributions DO I = -2, 6 DO M = 1, NMU DO N = 2, NX * The distributions shouldn't fall to zero too sharply IF ( (PTN(I,N,M).LT.DLOG(1D-14)) . AND. $ (PTN(I,N-1,M).GT.DLOG(1D-10)) ) $ PTN(I,N,M) = (PTN(I,N-1,M)+PTN(I,N,M)) / 2D0 ENDDO ENDDO ENDDO * Check for smoothness by looking at angles DO I = -2, 6 DO M = 1, NMU DO N = 2, NX-1 IF (PTN(I,N+1,M).GT.DLOG(1D-14)) THEN A = PTN(I,N,M) - PTN(I,N-1,M) B = PTN(I,N+1,M) - PTN(I,N,M) ANGL = DATAN((B-A)/(1D0+B*A))/3.1415D0*180D0 IF (ANGL.LT.-20D0) THEN TG = DTAN( DSIGN(20D0 / 180D0 * 3.14159D0,ANGL)) AP = - TG BP = TG*PTN(I,N+1,M) + TG*PTN(I,N-1,M) + 2D0 CP = TG - TG*PTN(I,N+1,M)*PTN(I,N-1,M) $ - PTN(I,N+1,M) - PTN(I,N-1,M) PT = (-BP + DSQRT(BP**2-4D0*AP*CP)) $ /(2*AP) PTN(I,N,M) = PT ENDIF IF (ANGL.GT.15D0) THEN PT = (PTN(I,N+1,M) + PTN(I,N-1,M))/2D0 PTN(I,N,M) = PT ENDIF ENDIF ENDDO ENDDO ENDDO * The smoothened parton distributions are written out DO N = 1, NX DO M = 1, NMU WRITE(99,'(9G13.5)') (PTN(I,N,M), I=-2,6) ENDDO ENDDO CLOSE(99) ENDDO STOP END C============================================================================ C============================================================================ C CTEQ Parton Distribution Functions: Version 6.0 C January 24, 2002 C C Ref: "New Generation of Parton Distributions with C Uncertainties from Global QCD Analysis" C By: J. Pumplin, D.R. Stump, J.Huston, H.L. Lai, P. Nadolsky, W.K. Tung C hep-ph/0201195 C C This package contains 3 standard sets of CTEQ6 PDF's and 40 up/down sets C with respect to CTEQ6M PDF's. Details are: C --------------------------------------------------------------------------- C Iset PDF Description Alpha_s(Mz)**Lam4 Lam5 Table_File C --------------------------------------------------------------------------- C 1 CTEQ6M Standard MSbar scheme 0.118 326 226 cteq6m.tbl C 2 CTEQ6D Standard DIS scheme 0.118 326 226 cteq6d.tbl C 3 CTEQ6L Leading Order 0.118** 326** 226 cteq6l.tbl C ------------------------------ C 1xx CTEQ6M1xx +/- w.r.t. CTEQ6M 0.118 326 226 cteq6m1xx.tbl C (where xx=01--40) C --------------------------------------------------------------------------- C ** ALL fits are obtained by using the same coupling strength \alpha_s(Mz)=0.118; C and the NLO running \alpha_s formula. For the LO fit, the evolution of the PDF C and the hard cross sections are calculated at LO. More detailed discussions are C given in hep-ph/0201195. C C The table grids are generated for 10^-6 < x < 1 and 1.3 < Q < 10,000 (GeV). C PDF values outside of the above range are returned using extrapolation. C Lam5 (Lam4) represents Lambda value (in MeV) for 5 (4) flavors. C The matching alpha_s between 4 and 5 flavors takes place at Q=4.5 GeV, C which is defined as the bottom quark mass, whenever it can be applied. C C The Table_Files are assumed to be in the working directory. C C Before using the PDF, it is necessary to do the initialization by C Call SetCtq6(Iset) C where Iset is the desired PDF specified in the above table. C C The function Ctq6Pdf (Iparton, X, Q) C returns the parton distribution inside the proton for parton [Iparton] C at [X] Bjorken_X and scale [Q] (GeV) in PDF set [Iset]. C Iparton is the parton label (5, 4, 3, 2, 1, 0, -1, ......, -5) C for (b, c, s, d, u, g, u_bar, ..., b_bar), C C For detailed information on the parameters used, e.q. quark masses, C QCD Lambda, ... etc., see info lines at the beginning of the C Table_Files. C C These programs, as provided, are in double precision. By removing the C "Implicit Double Precision" lines, they can also be run in single C precision. C C If you have detailed questions concerning these CTEQ6 distributions, C or if you find problems/bugs using this package, direct inquires to C Pumplin@pa.msu.edu or Tung@pa.msu.edu. C C=========================================================================== Function Ctq6Pdf (Iparton, X, Q) Implicit Double Precision (A-H,O-Z) Logical Warn Common > / CtqPar2 / Nx, Nt, NfMx > / QCDtable / Alambda, Nfl, Iorder Data Warn /.true./ save Warn If (X .lt. 0D0 .or. X .gt. 1D0) Then Print *, 'X out of range in Ctq6Pdf: ', X Stop Endif If (Q .lt. Alambda) Then Print *, 'Q out of range in Ctq6Pdf: ', Q Stop Endif If ((Iparton .lt. -NfMx .or. Iparton .gt. NfMx)) Then If (Warn) Then C put a warning for calling extra flavor. Warn = .false. Print *, 'Warning: Iparton out of range in Ctq6Pdf: ' > , Iparton Endif Ctq6Pdf = 0D0 Return Endif Ctq6Pdf = PartonX6 (Iparton, X, Q) if(Ctq6Pdf.lt.0.D0) Ctq6Pdf = 0.D0 Return C ******************** End Subroutine SetCtq6 (Iset) Implicit Double Precision (A-H,O-Z) Parameter (Isetmax0=3) Character Flnm(Isetmax0)*6, nn*3, Tablefile*40 Data (Flnm(I), I=1,Isetmax0) > / 'cteq6m', 'cteq6d', 'cteq6l' / Data Isetold, Isetmin0, Isetmin1, Isetmax1 /-987,1,101,140/ save C If data file not initialized, do so. If(Iset.ne.Isetold) then IU= NextUn() If (Iset.ge.Isetmin0 .and. Iset.le.Isetmax0) Then Tablefile=Flnm(Iset)//'.tbl' Elseif (Iset.ge.Isetmin1 .and. Iset.le.Isetmax1) Then write(nn,'(I3)') Iset Tablefile=Flnm(1)//nn//'.tbl' Else Print *, 'Invalid Iset number in SetCtq6 :', Iset Stop Endif Open(IU, File=Tablefile, Status='OLD', Err=100) 21 Call ReadTbl (IU) Close (IU) Isetold=Iset Endif Return 100 Print *, ' Data file ', Tablefile, ' cannot be opened ' >//'in SetCtq6!!' Stop C ******************** End Subroutine ReadTbl (Nu) Implicit Double Precision (A-H,O-Z) Character Line*80 PARAMETER (MXX = 96, MXQ = 20, MXF = 5) PARAMETER (MXPQX = (MXF + 3) * MXQ * MXX) Common > / CtqPar1 / Al, XV(0:MXX), TV(0:MXQ), UPD(MXPQX) > / CtqPar2 / Nx, Nt, NfMx > / XQrange / Qini, Qmax, Xmin > / QCDtable / Alambda, Nfl, Iorder > / Masstbl / Amass(6) Read (Nu, '(A)') Line Read (Nu, '(A)') Line Read (Nu, *) Dr, Fl, Al, (Amass(I),I=1,6) Iorder = Nint(Dr) Nfl = Nint(Fl) Alambda = Al Read (Nu, '(A)') Line Read (Nu, *) NX, NT, NfMx Read (Nu, '(A)') Line Read (Nu, *) QINI, QMAX, (TV(I), I =0, NT) Read (Nu, '(A)') Line Read (Nu, *) XMIN, (XV(I), I =0, NX) Do 11 Iq = 0, NT TV(Iq) = Log(Log (TV(Iq) /Al)) 11 Continue C C Since quark = anti-quark for nfl>2 at this stage, C we Read out only the non-redundent data points C No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence) Nblk = (NX+1) * (NT+1) Npts = Nblk * (NfMx+3) Read (Nu, '(A)') Line Read (Nu, *, IOSTAT=IRET) (UPD(I), I=1,Npts) Return C **************************** End Function NextUn() C Returns an unallocated FORTRAN i/o unit. Logical EX C Do 10 N = 10, 300 INQUIRE (UNIT=N, OPENED=EX) If (.NOT. EX) then NextUn = N Return Endif 10 Continue Stop ' There is no available I/O unit. ' C ************************* End C SUBROUTINE POLINT (XA,YA,N,X,Y,DY) IMPLICIT DOUBLE PRECISION (A-H, O-Z) C Adapted from "Numerical Recipes" PARAMETER (NMAX=10) DIMENSION XA(N),YA(N),C(NMAX),D(NMAX) NS=1 DIF=ABS(X-XA(1)) DO 11 I=1,N DIFT=ABS(X-XA(I)) IF (DIFT.LT.DIF) THEN NS=I DIF=DIFT ENDIF C(I)=YA(I) D(I)=YA(I) 11 CONTINUE Y=YA(NS) NS=NS-1 DO 13 M=1,N-1 DO 12 I=1,N-M HO=XA(I)-X HP=XA(I+M)-X W=C(I+1)-D(I) DEN=HO-HP IF(DEN.EQ.0.)PAUSE DEN=W/DEN D(I)=HP*DEN C(I)=HO*DEN 12 CONTINUE IF (2*NS.LT.N-M)THEN DY=C(NS+1) ELSE DY=D(NS) NS=NS-1 ENDIF Y=Y+DY 13 CONTINUE RETURN END Function PartonX6 (IPRTN, XX, QQ) c Given the parton distribution function in the array U in c COMMON / PEVLDT / , this routine interpolates to find c the parton distribution at an arbitray point in x and q. c Implicit Double Precision (A-H,O-Z) Parameter (MXX = 96, MXQ = 20, MXF = 5) Parameter (MXQX= MXQ * MXX, MXPQX = MXQX * (MXF+3)) Common > / CtqPar1 / Al, XV(0:MXX), TV(0:MXQ), UPD(MXPQX) > / CtqPar2 / Nx, Nt, NfMx > / XQrange / Qini, Qmax, Xmin Dimension fvec(4), fij(4) Dimension xvpow(0:mxx) Data OneP / 1.00001 / Data xpow / 0.3d0 / !**** choice of interpolation variable Data nqvec / 4 / Data ientry / 0 / Save ientry,xvpow c store the powers used for interpolation on first call... if(ientry .eq. 0) then ientry = 1 xvpow(0) = 0D0 do i = 1, nx xvpow(i) = xv(i)**xpow enddo endif X = XX Q = QQ tt = log(log(Q/Al)) c ------------- find lower end of interval containing x, i.e., c get jx such that xv(jx) .le. x .le. xv(jx+1)... JLx = -1 JU = Nx+1 11 If (JU-JLx .GT. 1) Then JM = (JU+JLx) / 2 If (X .Ge. XV(JM)) Then JLx = JM Else JU = JM Endif Goto 11 Endif C Ix 0 1 2 Jx JLx Nx-2 Nx C |---|---|---|...|---|-x-|---|...|---|---| C x 0 Xmin x 1 C If (JLx .LE. -1) Then Print '(A,1pE12.4)', 'Severe error: x <= 0 in PartonX6! x = ', x Stop ElseIf (JLx .Eq. 0) Then Jx = 0 Elseif (JLx .LE. Nx-2) Then C For interrior points, keep x in the middle, as shown above Jx = JLx - 1 Elseif (JLx.Eq.Nx-1 .or. x.LT.OneP) Then C We tolerate a slight over-shoot of one (OneP=1.00001), C perhaps due to roundoff or whatever, but not more than that. C Keep at least 4 points >= Jx Jx = JLx - 2 Else Print '(A,1pE12.4)', 'Severe error: x > 1 in PartonX6! x = ', x Stop Endif C ---------- Note: JLx uniquely identifies the x-bin; Jx does not. C This is the variable to be interpolated in ss = x**xpow If (JLx.Ge.2 .and. JLx.Le.Nx-2) Then c 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 c 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 c --------------Now find lower end of interval containing Q, i.e., c get jq such that qv(jq) .le. q .le. qv(jq+1)... JLq = -1 JU = NT+1 12 If (JU-JLq .GT. 1) Then JM = (JU+JLq) / 2 If (tt .GE. TV(JM)) Then JLq = JM Else JU = JM Endif Goto 12 Endif If (JLq .LE. 0) Then Jq = 0 Elseif (JLq .LE. Nt-2) Then C keep q in the middle, as shown above Jq = JLq - 1 Else C JLq .GE. Nt-1 case: Keep at least 4 points >= Jq. Jq = Nt - 3 Endif C This is the interpolation variable in Q If (JLq.GE.1 .and. JLq.LE.Nt-2) Then c store the lattice points in t... tvec1 = Tv(jq) tvec2 = Tv(jq+1) tvec3 = Tv(jq+2) tvec4 = Tv(jq+3) 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 c get the pdf function values at the lattice points... If (Iprtn .GE. 3) Then Ip = - Iprtn Else Ip = Iprtn EndIf jtmp = ((Ip + NfMx)*(NT+1)+(jq-1))*(NX+1)+jx+1 Do it = 1, nqvec J1 = jtmp + it*(NX+1) If (Jx .Eq. 0) Then C For the first 4 x points, interpolate x^2*f(x,Q) C This applies to the two lowest bins JLx = 0, 1 C We can not put the JLx.eq.1 bin into the "interrior" section C (as we do for q), since Upd(J1) is undefined. fij(1) = 0 fij(2) = Upd(J1+1) * XV(1)**2 fij(3) = Upd(J1+2) * XV(2)**2 fij(4) = Upd(J1+3) * XV(3)**2 C C Use Polint which allows x to be anywhere w.r.t. the grid Call Polint (XVpow(0), Fij(1), 4, ss, Fx, Dfx) If (x .GT. 0D0) Fvec(it) = Fx / x**2 C Pdf is undefined for x.eq.0 ElseIf (JLx .Eq. Nx-1) Then C This is the highest x bin: Call Polint (XVpow(Nx-3), Upd(J1), 4, ss, Fx, Dfx) Fvec(it) = Fx Else C for all interior points, use Jon's in-line function C This applied to (JLx.Ge.2 .and. JLx.Le.Nx-2) sf2 = Upd(J1+1) sf3 = Upd(J1+2) g1 = sf2*const1 - sf3*const2 g4 = -sf2*const3 + sf3*const4 Fvec(it) = (const5*(Upd(J1)-g1) & + const6*(Upd(J1+3)-g4) & + sf2*sy3 - sf3*sy2) / s23 Endif enddo C We now have the four values Fvec(1:4) c interpolate in t... If (JLq .LE. 0) Then C 1st Q-bin, as well as extrapolation to lower Q Call Polint (TV(0), Fvec(1), 4, tt, ff, Dfq) ElseIf (JLq .GE. Nt-1) Then C Last Q-bin, as well as extrapolation to higher Q Call Polint (TV(Nt-3), Fvec(1), 4, tt, ff, Dfq) Else C Interrior bins : (JLq.GE.1 .and. JLq.LE.Nt-2) C which include JLq.Eq.1 and JLq.Eq.Nt-2, since Upd is defined for C 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 PartonX6 = ff Return C ******************** End