LHAPDF/lhapdf5.5.1/src/eps09.f

   1 !****************************************************************************
   2 !
   3 !           EPS09.f
   4 !
   5 ! An interface for the scale dependent nuclear modifications
   6 !       R_f^A(x,Q) = f_A(x,Q)/f_p(x,Q)
   7 ! where f_A is the distribution of the parton flavour f for a PROTON in a
   8 ! nucleus A, and f_p is the corresponding parton distribution in the
   9 ! free proton.
  10 !
  11 ! When using this interface, please refer to:
  12 !
  13 ! K.J. Eskola, H. Paukkunen and C.A. Salgado,
  14 ! "EPS09 - a New Generation of NLO and LO Nuclear Parton Distribution Functions,"
  15 ! Published as JHEP04(2009) 065.
  16 ! Eprint: arXiv:0902.4154 [hep-ph].
  17 !
  18 ! Questions & comments to:
  19 !   hannu.paukkunen@phys.jyu.fi
  20 !   kari.eskola@phys.jyu.fi
  21 !   carlos.salgado@usc.es
  22 !
  23 ! ***************************************************************************
  24 ! Instructions:
  25 !
  26 ! For given input values of
  27 !
  28 !     order: 1=LO, 2=NLO   ; integer
  29 !     pset : 1...31        ; integer
  30 !            1     = central fit
  31 !            2,3   = error sets S{+1}, S{-1}
  32 !            4,5   = error sets S{+2}, S{-2}
  33 !            ...   ...
  34 !            30,31 = error sets {S+15}, {S-15}
  35 !     A    : atomic number ; integer
  36 !     x    : Bjorken-x     ; double precision
  37 !     Q    : scale in GeV  ; double precision
  38 !
  39 ! the command
  40 !
  41 !   Call EPS09(order, pset, A, x, Q, ruv, rdv, ru, rd, rs, rc, rb, rg)
  42 !
  43 ! returns the bound proton nuclear corrections R_f^A(x,Q)
  44 ! (in double precision) for
  45 !
  46 !   ruv = up valence
  47 !   rdv = down valence
  48 !   ru  = up sea
  49 !   rd  = down sea
  50 !   rs  = strange
  51 !   rc  = charm
  52 !   rb  = bottom
  53 !   rg  = gluons
  54 !
  55 ! The nuclear corrections for bound neutrons can be obtained
  56 ! by the isospin symmetry, e.g. the total up quark distribution
  57 ! per nucleon in a nucleus A with Z protons is
  58 !
  59 !  u_A(x,Q) =    Z/A * [ruv*uV_p(x,Q) + ru*uSea_p(x,Q)] +
  60 !            (A-Z)/A * [rdv*dV_p(x,Q) + rd*dSea_p(x,Q)]
  61 !
  62 ! Note that the parametrization should only be applied at the
  63 ! kinematical domain
  64 !
  65 !             1e-6 <= x <= 1
  66 !              1.3 <= Q <= 1000 GeV.
  67 !
  68 ! No warning message is displayed if these limits are
  69 ! exceeded, and outside these boundaries the modifications
  70 ! are frozen to the boundary values, i.e
  71 !
  72 !   for Q > 1000, the modifications at Q=1000 are returned,
  73 !   for Q < 1.3,  the modifications at Q=1.3 are returned,
  74 !   for x > 1,    the modifications at x=1 are returned
  75 !   for x < 1e-6, the modifications at x=1e-6 are returned,
  76 !
  77 ! The data used by the program for required order
  78 ! and atomic number A, are stored in separate files
  79 !
  80 !   LO : EPS09LOR_A
  81 !   NLO: EPS09NLOR_A
  82 !
  83 ! which must be located in the working directory.
  84 !
  85 ! The error bands for absolute cross-sections and for
  86 ! their nuclear ratios should be computed as explained
  87 ! in Secs. 2.5 and 4 of arXiv:0902.4154 [hep-ph]. For
  88 ! the absolute cross sections, both the errors in the
  89 ! free-proton PDFs f_p(x,Q) and the errors in
  90 ! the modifications R_f^A(x,Q) should be accounted for.
  91 ! For the nuclear ratios, it is sufficient to account only
  92 ! for the errors in the modifications R_f^A(x,Q).
  93 !
  94 ! *********************************************************
  95 ! *********************************************************
  96
  97
  98
  99       Subroutine EPS09(order, pset, AAA, xxx, QQQ,ruv, rdv, ru, rd, rs, rc, rb, rg)
 100
 101       Implicit none
 102       Double precision :: ruv, rdv, ru, rd, rs, rc, rb, rg, QQQ, xxx
 103       Double precision :: LSTEP, x, Q, Q2, allvalues(1:31,1:8,0:50,0:50)
 104       Double precision :: x_i=0.000001, arg(4), fu(4), res, fg(3)
 105       Double precision :: result(9), dummy
 106       Double precision :: realQ, Q2min=1.69, Q2max=1000000.0, Qsteps=50.0
 107       Double precision :: n_x, zero=0.0
 108
 109       Character (Len=50) filenimi
 110
 111       Integer :: xlinsteps=25, xlogsteps=25, startline, lineno
 112       Integer :: k, p, t, Qpoint, xpoint, pset, iovar
 113       Integer :: setnumber,j, A, openchannel, order, AAA
 114       Integer :: psetlast = -10, Alast = -10, orderlast = -10
 115
 116       character*512 dirpath,setpath
 117
 118       save Alast
 119       save psetlast
 120       save orderlast
 121       save allvalues
 122
 123 ! *********************************************
 124 ! Stop if the set specifications are wrong ones
 125 ! *********************************************
 126
 127       If (order .NE. 1 .and. order .NE. 2) then
 128       Write(*,*) 'Wrong order!'
 129       Write(*,*) 'LO : order = 1'
 130       Write(*,*) 'NLO: order = 2'
 131       Stop
 132       End If
 133
 134       If (pset .LT. 1 .or. pset .GT. 31) then
 135       Write(*,*) 'Wrong set!'
 136       Write(*,*) 'Central set: pset = 1'
 137       Write(*,*) 'Error sets : pset = 2...31'
 138       Stop
 139       End If
 140
 141 ! ********************************
 142 ! Make sure not to change any
 143 ! specifications given by the user
 144 ! ********************************
 145
 146       A  = AAA
 147       x  = xxx
 148       Q  = QQQ
 149       Q2 = Q*Q
 150
 151 ! *******************************
 152 ! Freeze x if it's < 10E-6 or > 1
 153 ! *******************************
 154
 155       If (x .LT. x_i) Then
 156       x = x_i
 157       End If
 158       If (x .GT. 1) Then
 159       x = 1.0
 160       End If
 161
 162 ! ************************************
 163 ! Freeze Q^2 if it's < 1.69 or > 10E+6
 164 ! ************************************
 165
 166       If (Q2 .LT. Q2min) Then
 167       Q2 = Q2min
 168       End If
 169       If (Q2 .GT. Q2max) Then
 170       Q2 = Q2max
 171       End If
 172
 173 ! If the set specifications have been changed, read the tables again
 174
 175       If (A .NE. Alast .or. order .NE. orderlast) Then
 176
 177 !      Write(*,*) 'Set changed!'
 178
 179 ! Read the table1
 180
 181       If (order .EQ. 1) then
 182
 183         If (A < 10) Then
 184         Write(filenimi,'("EPS09LOR_", I1)'), A
 185         Else If (A < 100) Then
 186         Write(filenimi,'("EPS09LOR_", I2)'), A
 187         Else If (A < 1000) Then
 188         Write(filenimi,'("EPS09LOR_", I3)'), A
 189         End If
 190
 191       Else
 192
 193         If (A < 10) Then
 194         Write(filenimi,'("EPS09NLOR_", I1)'), A
 195         Else If (A < 100) Then
 196         Write(filenimi,'("EPS09NLOR_", I2)'), A
 197         Else If (A < 1000) Then
 198         Write(filenimi,'("EPS09NLOR_", I3)'), A
 199         End If
 200
 201       End If
 202
 203       call getdirpath(dirpath)
 204       setpath=dirpath(:LEN_TRIM(dirpath))//"/"//filenimi(:LEN_TRIM(filenimi))//".LHgrid"
 205       print *,setpath(:LEN_TRIM(setpath))
 206       OPEN (11, file = setpath(:LEN_TRIM(setpath)), status='OLD', IOSTAT=iovar)
 207
 208       If (iovar .NE. 0) Then
 209       Write(*,*) 'Missing file: ',filenimi
 210       stop
 211       End If
 212
 213       Do setnumber = 1, 31
 214
 215       Do k = 0,50
 216
 217       Read(11,*) dummy
 218
 219       Do t = 0,50
 220
 221       Read(11,*) (allvalues(setnumber,p,k,t), p=1,8)
 222
 223       End Do
 224       End Do
 225
 226       End Do
 227
 228       Close(11)
 229
 230       psetlast  = pset
 231       Alast     = A
 232       orderlast = order
 233
 234       End If
 235
 236 ! Find out the position in the loglog Q^2-grid
 237
 238       realQ  = Qsteps * (log(log(Q2)/log(Q2min)))/(log(log(Q2max)/log(Q2min)))
 239       Qpoint = Aint(realQ)
 240
 241       If (Qpoint .LE. 0) Then
 242          Qpoint = 1
 243       End If
 244       If (Qpoint .GE. Anint(Qsteps)-1) Then
 245          Qpoint = Anint(Qsteps)-1
 246       End If
 247
 248       LSTEP = (1.0/(xlogsteps)) * LOG(0.1/x_i)
 249
 250 ! *********************
 251 ! Interpolate the grids
 252 ! *********************
 253
 254       Do t=1,8
 255
 256 ! Find the position in the x-grid
 257
 258       If (x .LE. 0.1) then
 259          n_x  = ((1.0/LSTEP) * Log(x/x_i))
 260        xpoint = Aint(n_x)
 261       Else
 262        n_x    = ((x-0.1)*xlinsteps/(1.0-0.1) + xlogsteps)
 263        xpoint = Aint(n_x)
 264       End If
 265
 266       If (xpoint .LE. 0) Then
 267         xpoint = 1
 268       End If
 269
 270       If (t .EQ. 1 .or. t .EQ. 2) Then
 271          If (xpoint .GE. (xlinsteps+xlogsteps)-4) Then
 272          xpoint = (xlinsteps+xlogsteps)-4
 273          End If
 274       End If
 275
 276       If (t .EQ. 3 .or. t .EQ. 4 .or. t .EQ. 5 .or. t .EQ. 6 .or. t .EQ. 7) Then
 277         If (xpoint .GE. (xlinsteps+xlogsteps)-7) Then
 278         xpoint = (xlinsteps+xlogsteps)-7
 279         End If
 280       End If
 281
 282       If (t .EQ. 8) Then
 283         If (xpoint .GE. (xlinsteps+xlogsteps)-4) Then
 284         xpoint = (xlinsteps+xlogsteps)-4
 285         End If
 286       End If
 287
 288       Do k = 1, 4
 289       If (xpoint-2+k .LT. xlogsteps) Then
 290       arg(k) = (x_i) * exp(LSTEP * (xpoint-2+k))
 291       Else
 292       arg(k) = 0.1 + (xpoint-2+k-xlogsteps) * (1-0.1)/xlinsteps
 293       End If
 294       End Do
 295
 296       Do j=1,3
 297
 298       fu(1) = allvalues(pset,t,Qpoint-2+j,xpoint-1)
 299       fu(2) = allvalues(pset,t,Qpoint-2+j,xpoint)
 300       fu(3) = allvalues(pset,t,Qpoint-2+j,xpoint+1)
 301       fu(4) = allvalues(pset,t,Qpoint-2+j,xpoint+2)
 302       Call luovi(fu,arg,4,x,res)
 303       fg(j) = res
 304
 305 ! *****************************************
 306 ! *****************************************
 307
 308       End Do
 309
 310       arg(1) = Qpoint-1
 311       arg(2) = Qpoint
 312       arg(3) = Qpoint+1
 313
 314       Call luovi(fg,arg,3,realQ,res)
 315
 316       result(t) = res
 317
 318       End Do
 319
 320       ruv = max(result(1),zero)
 321       rdv = max(result(2),zero)
 322       ru  = max(result(3),zero)
 323       rd  = max(result(4),zero)
 324       rs  = max(result(5),zero)
 325       rc  = max(result(6),zero)
 326       rb  = max(result(7),zero)
 327       rg  = max(result(8),zero)
 328
 329 200   Continue
 330
 331       End Subroutine EPS09
 332
 333 ! ********************************
 334 ! Modified version of Cern Library
 335 ! interpolation routine E100
 336 ! ********************************
 337
 338       SUBROUTINE luovi(F,ARG,MMM,Z,SUM)
 339
 340       Implicit none
 341       INTEGER :: MMM
 342       Double precision  :: F(MMM), ARG(MMM), COF(MMM), SUM, Z
 343       INTEGER :: M, MM, I, J, JNDEX, INDEX
 344
 345       MM = MIN(MMM, 20)
 346       M = MM - 1
 347       DO 1780 I= 1, MM
 348       COF(I) = F(I)
 349  1780 Continue
 350       DO 1800 I= 1, M
 351       DO 1790 J= I, M
 352          JNDEX = MM - J
 353          INDEX = JNDEX + I
 354          COF(INDEX) = (COF(INDEX)-COF(INDEX-1))/(ARG(INDEX)-ARG(JNDEX))
 355  1790 CONTINUE
 356  1800 CONTINUE
 357       SUM = COF(MM)
 358       DO 1810 I= 1, M
 359          INDEX = MM - I
 360          SUM = (Z-ARG(INDEX))*SUM + COF(INDEX)
 361  1810 CONTINUE
 362
 363       End SUBROUTINE luovi
 364