eps09 added.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 29 Aug 2012 15:33:50 +0000 (15:33 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 29 Aug 2012 15:33:50 +0000 (15:33 +0000)
LHAPDF/lhapdf5.5.1/src/eps09.f [new file with mode: 0644]

diff --git a/LHAPDF/lhapdf5.5.1/src/eps09.f b/LHAPDF/lhapdf5.5.1/src/eps09.f
new file mode 100644 (file)
index 0000000..eadc639
--- /dev/null
@@ -0,0 +1,364 @@
+!****************************************************************************\r
+!\r
+!           EPS09.f\r
+!\r
+! An interface for the scale dependent nuclear modifications\r
+!       R_f^A(x,Q) = f_A(x,Q)/f_p(x,Q) \r
+! where f_A is the distribution of the parton flavour f for a PROTON in a\r
+! nucleus A, and f_p is the corresponding parton distribution in the \r
+! free proton.\r
+!  \r
+! When using this interface, please refer to:\r
+!  \r
+! K.J. Eskola, H. Paukkunen and C.A. Salgado,\r
+! "EPS09 - a New Generation of NLO and LO Nuclear Parton Distribution Functions,"\r
+! Published as JHEP04(2009) 065.\r
+! Eprint: arXiv:0902.4154 [hep-ph].\r
+!\r
+! Questions & comments to:\r
+!   hannu.paukkunen@phys.jyu.fi\r
+!   kari.eskola@phys.jyu.fi\r
+!   carlos.salgado@usc.es\r
+! \r
+! ***************************************************************************\r
+! Instructions:\r
+!\r
+! For given input values of\r
+!\r
+!     order: 1=LO, 2=NLO   ; integer\r
+!     pset : 1...31        ; integer\r
+!            1     = central fit\r
+!            2,3   = error sets S{+1}, S{-1}\r
+!            4,5   = error sets S{+2}, S{-2}\r
+!            ...   ...\r
+!            30,31 = error sets {S+15}, {S-15}\r
+!     A    : atomic number ; integer\r
+!     x    : Bjorken-x     ; double precision\r
+!     Q    : scale in GeV  ; double precision\r
+!\r
+! the command \r
+!\r
+!   Call EPS09(order, pset, A, x, Q, ruv, rdv, ru, rd, rs, rc, rb, rg)\r
+!\r
+! returns the bound proton nuclear corrections R_f^A(x,Q)\r
+! (in double precision) for\r
+!   \r
+!   ruv = up valence\r
+!   rdv = down valence\r
+!   ru  = up sea\r
+!   rd  = down sea\r
+!   rs  = strange\r
+!   rc  = charm\r
+!   rb  = bottom\r
+!   rg  = gluons\r
+!\r
+! The nuclear corrections for bound neutrons can be obtained\r
+! by the isospin symmetry, e.g. the total up quark distribution\r
+! per nucleon in a nucleus A with Z protons is\r
+!\r
+!  u_A(x,Q) =    Z/A * [ruv*uV_p(x,Q) + ru*uSea_p(x,Q)] +\r
+!            (A-Z)/A * [rdv*dV_p(x,Q) + rd*dSea_p(x,Q)]\r
+!\r
+! Note that the parametrization should only be applied at the\r
+! kinematical domain\r
+!\r
+!             1e-6 <= x <= 1\r
+!              1.3 <= Q <= 1000 GeV.\r
+!\r
+! No warning message is displayed if these limits are\r
+! exceeded, and outside these boundaries the modifications\r
+! are frozen to the boundary values, i.e\r
+!\r
+!   for Q > 1000, the modifications at Q=1000 are returned,\r
+!   for Q < 1.3,  the modifications at Q=1.3 are returned,\r
+!   for x > 1,    the modifications at x=1 are returned\r
+!   for x < 1e-6, the modifications at x=1e-6 are returned,\r
+!\r
+! The data used by the program for required order\r
+! and atomic number A, are stored in separate files\r
+!\r
+!   LO : EPS09LOR_A\r
+!   NLO: EPS09NLOR_A\r
+!\r
+! which must be located in the working directory.\r
+!\r
+! The error bands for absolute cross-sections and for\r
+! their nuclear ratios should be computed as explained\r
+! in Secs. 2.5 and 4 of arXiv:0902.4154 [hep-ph]. For\r
+! the absolute cross sections, both the errors in the\r
+! free-proton PDFs f_p(x,Q) and the errors in\r
+! the modifications R_f^A(x,Q) should be accounted for.\r
+! For the nuclear ratios, it is sufficient to account only\r
+! for the errors in the modifications R_f^A(x,Q).\r
+!\r
+! *********************************************************\r
+! *********************************************************\r
+\r
+\r
+\r
+      Subroutine EPS09(order, pset, AAA, xxx, QQQ,ruv, rdv, ru, rd, rs, rc, rb, rg)\r
+\r
+      Implicit none\r
+      Double precision :: ruv, rdv, ru, rd, rs, rc, rb, rg, QQQ, xxx\r
+      Double precision :: LSTEP, x, Q, Q2, allvalues(1:31,1:8,0:50,0:50)\r
+      Double precision :: x_i=0.000001, arg(4), fu(4), res, fg(3)\r
+      Double precision :: result(9), dummy\r
+      Double precision :: realQ, Q2min=1.69, Q2max=1000000.0, Qsteps=50.0\r
+      Double precision :: n_x, zero=0.0\r
+\r
+      Character (Len=50) filenimi\r
+\r
+      Integer :: xlinsteps=25, xlogsteps=25, startline, lineno\r
+      Integer :: k, p, t, Qpoint, xpoint, pset, iovar\r
+      Integer :: setnumber,j, A, openchannel, order, AAA\r
+      Integer :: psetlast = -10, Alast = -10, orderlast = -10\r
\r
+      character*512 dirpath,setpath\r
+\r
+      save Alast\r
+      save psetlast\r
+      save orderlast\r
+      save allvalues\r
+\r
+! *********************************************\r
+! Stop if the set specifications are wrong ones\r
+! *********************************************\r
+\r
+      If (order .NE. 1 .and. order .NE. 2) then\r
+      Write(*,*) 'Wrong order!'\r
+      Write(*,*) 'LO : order = 1'\r
+      Write(*,*) 'NLO: order = 2'\r
+      Stop\r
+      End If\r
+\r
+      If (pset .LT. 1 .or. pset .GT. 31) then\r
+      Write(*,*) 'Wrong set!'\r
+      Write(*,*) 'Central set: pset = 1'\r
+      Write(*,*) 'Error sets : pset = 2...31'\r
+      Stop\r
+      End If\r
+\r
+! ********************************\r
+! Make sure not to change any\r
+! specifications given by the user\r
+! ********************************\r
+\r
+      A  = AAA\r
+      x  = xxx\r
+      Q  = QQQ\r
+      Q2 = Q*Q \r
+\r
+! *******************************\r
+! Freeze x if it's < 10E-6 or > 1\r
+! *******************************\r
+\r
+      If (x .LT. x_i) Then\r
+      x = x_i\r
+      End If\r
+      If (x .GT. 1) Then\r
+      x = 1.0\r
+      End If\r
+\r
+! ************************************\r
+! Freeze Q^2 if it's < 1.69 or > 10E+6\r
+! ************************************\r
+\r
+      If (Q2 .LT. Q2min) Then\r
+      Q2 = Q2min\r
+      End If\r
+      If (Q2 .GT. Q2max) Then\r
+      Q2 = Q2max\r
+      End If\r
+\r
+! If the set specifications have been changed, read the tables again\r
+\r
+      If (A .NE. Alast .or. order .NE. orderlast) Then\r
+\r
+!      Write(*,*) 'Set changed!'\r
+\r
+! Read the table1\r
+\r
+      If (order .EQ. 1) then\r
+\r
+        If (A < 10) Then\r
+        Write(filenimi,'("EPS09LOR_", I1)'), A\r
+        Else If (A < 100) Then\r
+        Write(filenimi,'("EPS09LOR_", I2)'), A\r
+        Else If (A < 1000) Then\r
+        Write(filenimi,'("EPS09LOR_", I3)'), A\r
+        End If\r
+\r
+      Else\r
+\r
+        If (A < 10) Then\r
+        Write(filenimi,'("EPS09NLOR_", I1)'), A\r
+        Else If (A < 100) Then\r
+        Write(filenimi,'("EPS09NLOR_", I2)'), A\r
+        Else If (A < 1000) Then\r
+        Write(filenimi,'("EPS09NLOR_", I3)'), A\r
+        End If\r
+\r
+      End If\r
+\r
+      call getdirpath(dirpath) \r
+      setpath=dirpath(:LEN_TRIM(dirpath))//"/"//filenimi(:LEN_TRIM(filenimi))//".LHgrid"\r
+      print *,setpath(:LEN_TRIM(setpath))\r
+      OPEN (11, file = setpath(:LEN_TRIM(setpath)), status='OLD', IOSTAT=iovar)\r
+\r
+      If (iovar .NE. 0) Then\r
+      Write(*,*) 'Missing file: ',filenimi\r
+      stop\r
+      End If\r
+\r
+      Do setnumber = 1, 31\r
+\r
+      Do k = 0,50\r
+\r
+      Read(11,*) dummy\r
+\r
+      Do t = 0,50\r
+\r
+      Read(11,*) (allvalues(setnumber,p,k,t), p=1,8)\r
+\r
+      End Do\r
+      End Do\r
+\r
+      End Do\r
+\r
+      Close(11)\r
+\r
+      psetlast  = pset\r
+      Alast     = A\r
+      orderlast = order\r
+\r
+      End If\r
+\r
+! Find out the position in the loglog Q^2-grid\r
+\r
+      realQ  = Qsteps * (log(log(Q2)/log(Q2min)))/(log(log(Q2max)/log(Q2min)))\r
+      Qpoint = Aint(realQ)\r
+\r
+      If (Qpoint .LE. 0) Then\r
+         Qpoint = 1\r
+      End If\r
+      If (Qpoint .GE. Anint(Qsteps)-1) Then\r
+         Qpoint = Anint(Qsteps)-1\r
+      End If\r
+\r
+      LSTEP = (1.0/(xlogsteps)) * LOG(0.1/x_i)\r
+\r
+! *********************\r
+! Interpolate the grids \r
+! *********************\r
+\r
+      Do t=1,8\r
+\r
+! Find the position in the x-grid\r
+\r
+      If (x .LE. 0.1) then\r
+         n_x  = ((1.0/LSTEP) * Log(x/x_i))\r
+       xpoint = Aint(n_x)\r
+      Else\r
+       n_x    = ((x-0.1)*xlinsteps/(1.0-0.1) + xlogsteps)\r
+       xpoint = Aint(n_x)\r
+      End If\r
+\r
+      If (xpoint .LE. 0) Then\r
+        xpoint = 1\r
+      End If\r
+\r
+      If (t .EQ. 1 .or. t .EQ. 2) Then\r
+         If (xpoint .GE. (xlinsteps+xlogsteps)-4) Then\r
+         xpoint = (xlinsteps+xlogsteps)-4\r
+         End If\r
+      End If\r
+\r
+      If (t .EQ. 3 .or. t .EQ. 4 .or. t .EQ. 5 .or. t .EQ. 6 .or. t .EQ. 7) Then\r
+        If (xpoint .GE. (xlinsteps+xlogsteps)-7) Then\r
+        xpoint = (xlinsteps+xlogsteps)-7\r
+        End If\r
+      End If\r
+       \r
+      If (t .EQ. 8) Then\r
+        If (xpoint .GE. (xlinsteps+xlogsteps)-4) Then\r
+        xpoint = (xlinsteps+xlogsteps)-4\r
+        End If\r
+      End If\r
+\r
+      Do k = 1, 4\r
+      If (xpoint-2+k .LT. xlogsteps) Then\r
+      arg(k) = (x_i) * exp(LSTEP * (xpoint-2+k))\r
+      Else\r
+      arg(k) = 0.1 + (xpoint-2+k-xlogsteps) * (1-0.1)/xlinsteps\r
+      End If\r
+      End Do\r
\r
+      Do j=1,3\r
+\r
+      fu(1) = allvalues(pset,t,Qpoint-2+j,xpoint-1)\r
+      fu(2) = allvalues(pset,t,Qpoint-2+j,xpoint)\r
+      fu(3) = allvalues(pset,t,Qpoint-2+j,xpoint+1)\r
+      fu(4) = allvalues(pset,t,Qpoint-2+j,xpoint+2)\r
+      Call luovi(fu,arg,4,x,res)\r
+      fg(j) = res\r
+\r
+! *****************************************\r
+! *****************************************\r
+\r
+      End Do\r
+\r
+      arg(1) = Qpoint-1\r
+      arg(2) = Qpoint\r
+      arg(3) = Qpoint+1\r
+\r
+      Call luovi(fg,arg,3,realQ,res)\r
+  \r
+      result(t) = res\r
+\r
+      End Do\r
+\r
+      ruv = max(result(1),zero)\r
+      rdv = max(result(2),zero)\r
+      ru  = max(result(3),zero)\r
+      rd  = max(result(4),zero)\r
+      rs  = max(result(5),zero)\r
+      rc  = max(result(6),zero)\r
+      rb  = max(result(7),zero)\r
+      rg  = max(result(8),zero)\r
+\r
+200   Continue\r
+\r
+      End Subroutine EPS09\r
+\r
+! ********************************\r
+! Modified version of Cern Library\r
+! interpolation routine E100\r
+! ********************************\r
+\r
+      SUBROUTINE luovi(F,ARG,MMM,Z,SUM)\r
+\r
+      Implicit none\r
+      INTEGER :: MMM\r
+      Double precision  :: F(MMM), ARG(MMM), COF(MMM), SUM, Z\r
+      INTEGER :: M, MM, I, J, JNDEX, INDEX\r
+\r
+      MM = MIN(MMM, 20)\r
+      M = MM - 1\r
+      DO 1780 I= 1, MM\r
+      COF(I) = F(I)\r
+ 1780 Continue\r
+      DO 1800 I= 1, M\r
+      DO 1790 J= I, M\r
+         JNDEX = MM - J\r
+         INDEX = JNDEX + I\r
+         COF(INDEX) = (COF(INDEX)-COF(INDEX-1))/(ARG(INDEX)-ARG(JNDEX))\r
+ 1790 CONTINUE\r
+ 1800 CONTINUE\r
+      SUM = COF(MM)\r
+      DO 1810 I= 1, M\r
+         INDEX = MM - I\r
+         SUM = (Z-ARG(INDEX))*SUM + COF(INDEX)\r
+ 1810 CONTINUE\r
+\r
+      End SUBROUTINE luovi\r
+\r