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