partonx6.F

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file partonx6.F

      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 Polint6 which allows x to be anywhere w.r.t. the grid

         Call polint6(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 polint6(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 polint6(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 polint6(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