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