PROGRAM MAKERNCTEQ5 * This is an interface to the smoothened CTEQ5 parton distribution * routine that produces a table in the standard form. * D. E. Soper and P. Anandam. 17 May 1999. DOUBLE PRECISION LAMBDA DOUBLE PRECISION X,LX,LXMIN,LXMAX,DELTALX DOUBLE PRECISION MU,LMU,LMUMIN,LMUMAX,DELTALMU INTEGER ISET EXTERNAL Ctq5Pdf, SetCtq5 DOUBLE PRECISION Ctq5Pdf CHARACTER*100 HEADER DOUBLE PRECISION LAMBDA(7) INTEGER LAMBDAFLAVORS(7) DOUBLE PRECISION XMIN,XMAX INTEGER NX DOUBLE PRECISION MUMIN,MUMAX INTEGER NMU CHARACTER*20 FILENAME(7) CHARACTER*20 NAME(7) DOUBLE PRECISION PTN(-2:6,64,64) DOUBLE PRECISION A, B, AP, BP, CP, TG, ANGL, PT * Loop variables INTEGER I, N,M DATA FILENAME / 'cteq5m.ptn', 'cteq5d.ptn', $ 'cteq5l.ptn', 'cteq5hj.ptn', 'cteq5hq.ptn', $ 'cteq5f3.ptn','cteq5f4.ptn' / DATA NAME /'CTEQ5M','CTEQ5D','CTEQ5L','CTEQ5HJ','CTEQ5HQ', $ 'CTEQ5F3','CTEQ5F4'/ DATA LAMBDA /0.226D0,0.226D0,0.146D0,0.226D0,0.226D0,0.395D0, $ 0.309D0/ DATA LAMBDAFLAVORS /5,5,5,5,5,3,4/ XMIN = 0.1000000000D-04 XMAX = 0.9500000000D0 NX = 32 MUMIN = 1D0 MUMAX = 10000D0 - 1D-6 NMU = 32 DO ISET=1,7 CALL SetCtq5(ISET) WRITE(*,*) 'Making ',FILENAME(ISET) OPEN(UNIT = 99, FILE = FILENAME(ISET), STATUS = 'UNKNOWN') HEADER = '! ' // NAME(ISET)(1:7) // $ ' distribution generated by DES and PA 17 May 1999' 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) = Ctq5Pdf(-2, X, MU) PTN(-1,N,M) = Ctq5Pdf(-1, X, MU) PTN(0,N,M) = Ctq5Pdf(0, X, MU) PTN(1,N,M) = Ctq5Pdf(1, X, MU) PTN(2,N,M) = Ctq5Pdf(2, X, MU) PTN(3,N,M) = Ctq5Pdf(3, X, MU) PTN(4,N,M) = Ctq5Pdf(4, X, MU) PTN(5,N,M) = Ctq5Pdf(5, X, MU) PTN(6,N,M) = Ctq5Pdf(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 CTEQ Parton Distribution Functions: Version 5.0 C March 1, 1999 C C Ref: "GLOBAL QCD ANALYSIS OF PARTON STRUCTURE OF THE NUCLEON: C CTEQ5 PPARTON DISTRIBUTIONS" C C hep-ph/9903282 C C These PDF's use quadratic interpolation of attached tables. A parametrized C version of the same PDF's without external tables is under construction. C They will become available later. C C This package contains 7 sets of CTEQ5 PDF's. Details are: C --------------------------------------------------------------------------- C Iset PDF Description Alpha_s(Mz) Lam4 Lam5 Table_File C --------------------------------------------------------------------------- C 1 CTEQ5M Standard MSbar scheme 0.118 326 226 cteq5m.tbl C 2 CTEQ5D Standard DIS scheme 0.118 326 226 cteq5d.tbl C 3 CTEQ5L Leading Order 0.127 192 146 cteq5l.tbl C 4 CTEQ5HJ Large-x gluon enhanced 0.118 326 226 cteq5hj.tbl C 5 CTEQ5HQ Heavy Quark 0.118 326 226 cteq5hq.tbl C 6 CTEQ5F3 Nf=3 FixedFlavorNumber 0.106 (Lam3=395) cteq5f3.tbl C 7 CTEQ5F4 Nf=4 FixedFlavorNumber 0.112 309 XXX cteq5f4.tbl C --------------------------------------------------------------------------- C C The available applied range is 10^-5 < x < 1 and 1.0 < Q < 10,000 (GeV). C Lam5 (Lam4, Lam3) represents Lambda value (in MeV) for 5 (4,3) 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 SetCtq5(Iset) C where Iset is the desired PDF specified in the above table. C C The function Ctq5Pdf (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 whereas CTEQ5F3 has, by definition, only 3 flavors and gluon; C CTEQ5F4 has only 4 flavors and gluon. 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 CTEQ5 distributions, C or if you find problems/bugs using this package, direct inquires to C Hung-Liang Lai(lai@phys.nthu.edu.tw) or Wu-Ki Tung(Tung@pa.msu.edu). C C=========================================================================== Function Ctq5Pdf (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 Ctq5Pdf: ', X Stop Endif If (Q .lt. Alambda) Then Print *, 'Q out of range in Ctq5Pdf: ', 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 Ctq5Pdf: ' > , Iparton Endif Ctq5Pdf = 0D0 Return Endif Ctq5Pdf = PartonX (Iparton, X, Q) if(Ctq5Pdf.lt.0.D0) Ctq5Pdf = 0.D0 Return C ******************** End FUNCTION PartonX (IPRTN, X, Q) C C Given the parton distribution function in the array Upd in C COMMON / CtqPar1 / , this routine fetches u(fl, x, q) at any value of C x and q using Mth-order polynomial interpolation for x and Ln(Q/Lambda). C IMPLICIT DOUBLE PRECISION (A-H, O-Z) C PARAMETER (MXX = 105, MXQ = 25, MXF = 6) PARAMETER (MXPQX = (MXF *2 +2) * MXQ * MXX) PARAMETER (M= 2, M1 = M + 1) C Logical First Common > / CtqPar1 / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX) > / CtqPar2 / Nx, Nt, NfMx > / XQrange / Qini, Qmax, Xmin C Dimension Fq(M1), Df(M1) Data First /.true./ save First C Work with Log (Q) QG = LOG (Q/AL) C Find lower end of interval containing X JL = -1 JU = Nx+1 11 If (JU-JL .GT. 1) Then JM = (JU+JL) / 2 If (X .GT. XV(JM)) Then JL = JM Else JU = JM Endif Goto 11 Endif Jx = JL - (M-1)/2 If (X .lt. Xmin .and. First ) Then First = .false. Print '(A, 2(1pE12.4))', > ' WARNING: X < Xmin, extrapolation used; X, Xmin =', X, Xmin If (Jx .LT. 0) Jx = 0 Elseif (Jx .GT. Nx-M) Then Jx = Nx - M Endif C Find the interval where Q lies JL = -1 JU = NT+1 12 If (JU-JL .GT. 1) Then JM = (JU+JL) / 2 If (QG .GT. QL(JM)) Then JL = JM Else JU = JM Endif Goto 12 Endif Jq = JL - (M-1)/2 If (Jq .LT. 0) Then Jq = 0 If (Q .lt. Qini) Print '(A, 2(1pE12.4))', > ' WARNING: Q < Qini, extrapolation used; Q, Qini =', Q, Qini Elseif (Jq .GT. Nt-M) Then Jq = Nt - M If (Q .gt. Qmax) Print '(A, 2(1pE12.4))', > ' WARNING: Q > Qmax, extrapolation used; Q, Qmax =', Q, Qmax Endif If (Iprtn .GE. 3) Then Ip = - Iprtn Else Ip = Iprtn EndIf C Find the off-set in the linear array Upd JFL = Ip + NfMx J0 = (JFL * (NT+1) + Jq) * (NX+1) + Jx C C Now interpolate in x for M1 Q's Do 21 Iq = 1, M1 J1 = J0 + (Nx+1)*(Iq-1) + 1 Call Polint (XV(Jx), Upd(J1), M1, X, Fq(Iq), Df(Iq)) 21 Continue C Finish off by interpolating in Q Call Polint (QL(Jq), Fq(1), M1, QG, Ftmp, Ddf) PartonX = Ftmp C RETURN C **************************** END Subroutine SetCtq5 (Iset) Implicit Double Precision (A-H,O-Z) Parameter (Isetmax=7) Character Flnm(Isetmax)*12, Tablefile*40 Data (Flnm(I), I=1,Isetmax) > / 'cteq5m.tbl', 'cteq5d.tbl', 'cteq5l.tbl', 'cteq5hj.tbl' > , 'cteq5hq.tbl', 'cteq5f3.tbl', 'cteq5f4.tbl' / Data Tablefile / 'test.tbl' / Data Isetold, Isetmin, Isettest / -987, 1, 911 / save C If data file not initialized, do so. If(Iset.ne.Isetold) then IU= NextUn() If (Iset .eq. Isettest) then Print* ,'Opening ', Tablefile 21 Open(IU, File=Tablefile, Status='OLD', Err=101) ElseIf (Iset.lt.Isetmin .or. Iset.gt.Isetmax) Then Print *, 'Invalid Iset number in SetCtq5 :', Iset Stop Else Tablefile=Flnm(Iset) Open(IU, File=Tablefile, Status='OLD', Err=100) Endif Call ReadTbl (IU) Close (IU) Isetold=Iset Endif Return 100 Print *, ' Data file ', Tablefile, ' cannot be opened ' >//'in SetCtq5!!' Stop 101 Print*, Tablefile, ' cannot be opened ' Print*, 'Please input the .tbl file:' Read (*,'(A)') Tablefile Goto 21 C ******************** End Subroutine ReadTbl (Nu) Implicit Double Precision (A-H,O-Z) Character Line*80 PARAMETER (MXX = 105, MXQ = 25, MXF = 6) PARAMETER (MXPQX = (MXF *2 +2) * MXQ * MXX) Common > / CtqPar1 / Al, XV(0:MXX), QL(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, (QL(I), I =0, NT) Read (Nu, '(A)') Line Read (Nu, *) XMIN, (XV(I), I =0, NX) Do 11 Iq = 0, NT QL(Iq) = Log (QL(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