Cteq61Pdf.f

Go to the documentation of this file.

C============================================================================
C                             April 10, 2002, v6.01
C                             February 23, 2003, v6.1
C
C   Ref[1]: "New Generation of Parton Distributions with Uncertainties from Global QCD Analysis"
C       By: J. Pumplin, D.R. Stump, J.Huston, H.L. Lai, P. Nadolsky, W.K. Tung
C       JHEP 0207:012(2002), hep-ph/0201195
C
C   Ref[2]: "Inclusive Jet Production, Parton Distributions, and the Search for New Physics"
C       By : D. Stump, J. Huston, J. Pumplin, W.K. Tung, H.L. Lai, S. Kuhlmann, J. Owens
C       hep-ph/0303013
C
C   This package contains
C   (1) 4 standard sets of CTEQ6 PDF's (CTEQ6M, CTEQ6D, CTEQ6L, CTEQ6L1) ;
C   (2) 40 up/down sets (with respect to CTEQ6M) for uncertainty studies from Ref[1];
C   (3) updated version of the above: CTEQ6.1M and its 40 up/down eigenvector sets from Ref[2].
C
C  The CTEQ6.1M set provides a global fit that is almost equivalent in every respect
C  to the published CTEQ6M, Ref[1], although some parton distributions (e.g., the gluon)
C  may deviate from CTEQ6M in some kinematic ranges by amounts that are well within the
C  specified uncertainties.
C  The more significant improvements of the new version are associated with some of the
C  40 eigenvector sets, which are made more symmetrical and reliable in (3), compared to (2).
 
C  Details about calling convention are:
C ---------------------------------------------------------------------------
C  Iset   PDF-set     Description       Alpha_s(Mz)**Lam4  Lam5   Table_File
C ===========================================================================
C Standard, "best-fit", sets:
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   4    CTEQ6L1  Leading Order           0.130**   215** 165    cteq6l1.tbl
C ============================================================================
C For uncertainty calculations using eigenvectors of the Hessian:
C ---------------------------------------------------------------
C     central + 40 up/down sets along 20 eigenvector directions
C                             -----------------------------
C                Original version, Ref[1]:  central fit: CTEQ6M (=CTEQ6M.00)
C                             -----------------------
C  1xx  CTEQ6M.xx  +/- sets               0.118     326   226    cteq6m1xx.tbl
C        where xx = 01-40: 01/02 corresponds to +/- for the 1st eigenvector, ... etc.
C        e.g. 100      is CTEQ6M.00 (=CTEQ6M),
C             101/102 are CTEQ6M.01/02, +/- sets of 1st eigenvector, ... etc.
C                              -----------------------
C                Updated version, Ref[2]:  central fit: CTEQ6.1M (=CTEQ61.00)
C                              -----------------------
C  2xx  CTEQ61.xx  +/- sets               0.118     326   226    ctq61.xx.tbl
C        where xx = 01-40: 01/02 corresponds to +/- for the 1st eigenvector, ... etc.
C        e.g. 200      is CTEQ61.00 (=CTEQ6.1M),
C             201/202 are CTEQ61.01/02, +/- sets of 1st eigenvector, ... etc.
C ===========================================================================
C   ** ALL fits are obtained by using the same coupling strength
C   \alpha_s(Mz)=0.118 and the NLO running \alpha_s formula, except CTEQ6L1
C   which uses the LO running \alpha_s and its value determined from the fit.
C   For the LO fits, the evolution of the PDF and the hard cross sections are
C   calculated at LO.  More detailed discussions are given in the references.
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   [email protected] or [email protected].
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.
             print *, 'Warning: Iparton out of range in Ctq6Pdf! '
             print *, 'Iparton, MxFlvN0: ', iparton, nfmx
         Endif
         ctq6pdf = 0d0
         Return
      Endif
 
      ctq6pdf = partonx6(iparton, x, q)
      if(ctq6pdf.lt.0.d0)  ctq6pdf = 0.d0
      if(isnan(ctq6pdf)) then
         if (warn) then
            print *, 'Warning: Ctq6Pdf is NaN!'
            print *, 'Iparton, MxFlvN0, X, Q: ', iparton, nfmx, x, q
         endif
         ctq6pdf = 0.d0
      endif
 
      if (warn) warn = .false.
 
      Return
 
C                             ********************
      End
 
      Subroutine setctq6 (Iset)
      Implicit Double Precision (a-h,o-z)
      parameter(isetmax0=5)
      Character Flnm(Isetmax0)*6, nn*3, Tablefile*40
      Data (flnm(i), i=1,isetmax0)
     > / 'cteq6m', 'cteq6d', 'cteq6l', 'cteq6l','ctq61.'/
      Data isetold, isetmin0, isetmin1, isetmax1 /-987,1,100,140/
      Data isetmin2,isetmax2 /200,240/
      save
 
C             If data file not initialized, do so.
      If(iset.ne.isetold) then
         iu= nextun()
         If (iset.ge.isetmin0 .and. iset.le.3) Then
            tablefile=flnm(iset)//'.tbl'
         Elseif (iset.eq.4) Then
            tablefile=flnm(iset)//'1.tbl'
         Elseif (iset.ge.isetmin1 .and. iset.le.isetmax1) Then
            write(nn,'(I3)') iset
            tablefile=flnm(1)//nn//'.tbl'
         Elseif (iset.ge.isetmin2 .and. iset.le.isetmax2) Then
            write(nn,'(I3)') iset
            tablefile=flnm(5)//nn(2:3)//'.tbl'
         Else
            print *, 'Invalid Iset number in SetCtq6 :', iset
            stop
         Endif
         Open(iu, file='pdfs/'//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