end-of-line normalization
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.5.1 / src / eps09.f
CommitLineData
7fac8669 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
329200 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