Line data Source code
1 : // PartonDistributions.cc is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 Torbjorn Sjostrand.
3 : // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4 : // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 :
6 : // Function definitions (not found in the header) for the PDF, LHAPDF,
7 : // GRV94L, CTEQ5L, MSTWpdf, CTEQ6pdf, GRVpiL, PomFix, PomH1FitAB,
8 : // PomH1Jets, Lepton and NNPDF classes.
9 :
10 : #include "Pythia8/PartonDistributions.h"
11 :
12 : namespace Pythia8 {
13 :
14 : //==========================================================================
15 :
16 : // Base class for parton distribution functions.
17 :
18 : //--------------------------------------------------------------------------
19 :
20 : // Resolve valence content for assumed meson. Possibly modified later.
21 :
22 : void PDF::setValenceContent() {
23 :
24 : // Subdivide meson by flavour content.
25 0 : if (idBeamAbs < 100 || idBeamAbs > 1000) return;
26 0 : int idTmp1 = idBeamAbs/100;
27 0 : int idTmp2 = (idBeamAbs/10)%10;
28 :
29 : // Find which is quark and which antiquark.
30 0 : if (idTmp1%2 == 0) {
31 0 : idVal1 = idTmp1;
32 0 : idVal2 = -idTmp2;
33 0 : } else {
34 0 : idVal1 = idTmp2;
35 0 : idVal2 = -idTmp1;
36 : }
37 0 : if (idBeam < 0) {
38 0 : idVal1 = -idVal1;
39 0 : idVal2 = -idVal2;
40 0 : }
41 :
42 : // Special case for Pomeron, to start off.
43 0 : if (idBeamAbs == 990) {
44 0 : idVal1 = 1;
45 0 : idVal2 = -1;
46 0 : }
47 0 : }
48 :
49 : //--------------------------------------------------------------------------
50 :
51 : // Standard parton densities.
52 :
53 : double PDF::xf(int id, double x, double Q2) {
54 :
55 : // Need to update if flavour, x or Q2 changed.
56 : // Use idSav = 9 to indicate that ALL flavours are up-to-date.
57 : // Assume that flavour and antiflavour always updated simultaneously.
58 0 : if ( (abs(idSav) != abs(id) && idSav != 9) || x != xSav || Q2 != Q2Sav)
59 0 : {idSav = id; xfUpdate(id, x, Q2); xSav = x; Q2Sav = Q2;}
60 :
61 : // Baryon and nondiagonal meson beams: only p, pbar, pi+, pi- for now.
62 0 : if (idBeamAbs == 2212 || idBeamAbs == 211) {
63 0 : int idNow = (idBeam > 0) ? id : -id;
64 : int idAbs = abs(id);
65 0 : if (idNow == 0 || idAbs == 21) return max(0., xg);
66 0 : if (idNow == 1) return max(0., xd);
67 0 : if (idNow == -1) return max(0., xdbar);
68 0 : if (idNow == 2) return max(0., xu);
69 0 : if (idNow == -2) return max(0., xubar);
70 0 : if (idNow == 3) return max(0., xs);
71 0 : if (idNow == -3) return max(0., xsbar);
72 0 : if (idAbs == 4) return max(0., xc);
73 0 : if (idAbs == 5) return max(0., xb);
74 0 : if (idAbs == 22) return max(0., xgamma);
75 0 : return 0.;
76 :
77 : // Baryon beams: n and nbar by isospin conjugation of p and pbar.
78 0 : } else if (idBeamAbs == 2112) {
79 0 : int idNow = (idBeam > 0) ? id : -id;
80 : int idAbs = abs(id);
81 0 : if (idNow == 0 || idAbs == 21) return max(0., xg);
82 0 : if (idNow == 1) return max(0., xu);
83 0 : if (idNow == -1) return max(0., xubar);
84 0 : if (idNow == 2) return max(0., xd);
85 0 : if (idNow == -2) return max(0., xdbar);
86 0 : if (idNow == 3) return max(0., xs);
87 0 : if (idNow == -3) return max(0., xsbar);
88 0 : if (idAbs == 4) return max(0., xc);
89 0 : if (idAbs == 5) return max(0., xb);
90 0 : if (idAbs == 22) return max(0., xgamma);
91 0 : return 0.;
92 :
93 : // Diagonal meson beams: only pi0, Pomeron for now.
94 0 : } else if (idBeam == 111 || idBeam == 990) {
95 : int idAbs = abs(id);
96 0 : if (id == 0 || idAbs == 21) return max(0., xg);
97 0 : if (id == idVal1 || id == idVal2) return max(0., xu);
98 0 : if (idAbs <= 2) return max(0., xubar);
99 0 : if (idAbs == 3) return max(0., xs);
100 0 : if (idAbs == 4) return max(0., xc);
101 0 : if (idAbs == 5) return max(0., xb);
102 0 : if (idAbs == 22) return max(0., xgamma);
103 0 : return 0.;
104 :
105 :
106 : // Lepton beam.
107 : } else {
108 0 : if (id == idBeam ) return max(0., xlepton);
109 0 : if (abs(id) == 22) return max(0., xgamma);
110 0 : return 0.;
111 : }
112 :
113 0 : }
114 :
115 : //--------------------------------------------------------------------------
116 :
117 : // Only valence part of parton densities.
118 :
119 : double PDF::xfVal(int id, double x, double Q2) {
120 :
121 : // Need to update if flavour, x or Q2 changed.
122 : // Use idSav = 9 to indicate that ALL flavours are up-to-date.
123 : // Assume that flavour and antiflavour always updated simultaneously.
124 0 : if ( (abs(idSav) != abs(id) && idSav != 9) || x != xSav || Q2 != Q2Sav)
125 0 : {idSav = id; xfUpdate(id, x, Q2); xSav = x; Q2Sav = Q2;}
126 :
127 : // Baryon and nondiagonal meson beams: only p, pbar, n, nbar, pi+, pi-.
128 0 : if (idBeamAbs == 2212) {
129 0 : int idNow = (idBeam > 0) ? id : -id;
130 0 : if (idNow == 1) return max(0., xdVal);
131 0 : if (idNow == 2) return max(0., xuVal);
132 0 : return 0.;
133 0 : } else if (idBeamAbs == 2112) {
134 0 : int idNow = (idBeam > 0) ? id : -id;
135 0 : if (idNow == 1) return max(0., xuVal);
136 0 : if (idNow == 2) return max(0., xdVal);
137 0 : return 0.;
138 0 : } else if (idBeamAbs == 211) {
139 0 : int idNow = (idBeam > 0) ? id : -id;
140 0 : if (idNow == 2 || idNow == -1) return max(0., xuVal);
141 0 : return 0.;
142 :
143 : // Diagonal meson beams: only pi0, Pomeron for now.
144 0 : } else if (idBeam == 111 || idBeam == 990) {
145 0 : if (id == idVal1 || id == idVal2) return max(0., xuVal);
146 0 : return 0.;
147 :
148 : // Lepton beam.
149 : } else {
150 0 : if (id == idBeam) return max(0., xlepton);
151 0 : return 0.;
152 : }
153 :
154 0 : }
155 :
156 : //--------------------------------------------------------------------------
157 :
158 : // Only sea part of parton densities.
159 :
160 : double PDF::xfSea(int id, double x, double Q2) {
161 :
162 : // Need to update if flavour, x or Q2 changed.
163 : // Use idSav = 9 to indicate that ALL flavours are up-to-date.
164 : // Assume that flavour and antiflavour always updated simultaneously.
165 0 : if ( (abs(idSav) != abs(id) && idSav != 9) || x != xSav || Q2 != Q2Sav)
166 0 : {idSav = id; xfUpdate(id, x, Q2); xSav = x; Q2Sav = Q2;}
167 :
168 : // Hadron beams.
169 0 : if (idBeamAbs > 100) {
170 0 : int idNow = (idBeam > 0) ? id : -id;
171 : int idAbs = abs(id);
172 0 : if (idNow == 0 || idAbs == 21) return max(0., xg);
173 0 : if (idBeamAbs == 2212) {
174 0 : if (idNow == 1) return max(0., xdSea);
175 0 : if (idNow == -1) return max(0., xdbar);
176 0 : if (idNow == 2) return max(0., xuSea);
177 0 : if (idNow == -2) return max(0., xubar);
178 0 : } else if (idBeamAbs == 2112) {
179 0 : if (idNow == 1) return max(0., xuSea);
180 0 : if (idNow == -1) return max(0., xubar);
181 0 : if (idNow == 2) return max(0., xdSea);
182 0 : if (idNow == -2) return max(0., xdbar);
183 : } else {
184 0 : if (idAbs <= 2) return max(0., xuSea);
185 : }
186 0 : if (idNow == 3) return max(0., xs);
187 0 : if (idNow == -3) return max(0., xsbar);
188 0 : if (idAbs == 4) return max(0., xc);
189 0 : if (idAbs == 5) return max(0., xb);
190 0 : if (idAbs == 22) return max(0., xgamma);
191 0 : return 0.;
192 :
193 : // Lepton beam.
194 : } else {
195 0 : if (abs(id) == 22) return max(0., xgamma);
196 0 : return 0.;
197 : }
198 :
199 0 : }
200 :
201 : //==========================================================================
202 :
203 : // Gives the GRV 94 L (leading order) parton distribution function set
204 : // in parametrized form. Authors: M. Glueck, E. Reya and A. Vogt.
205 : // Ref: M. Glueck, E. Reya and A. Vogt, Z.Phys. C67 (1995) 433.
206 :
207 : void GRV94L::xfUpdate(int , double x, double Q2) {
208 :
209 : // Common expressions. Constrain Q2 for which parametrization is valid.
210 : double mu2 = 0.23;
211 : double lam2 = 0.2322 * 0.2322;
212 0 : double s = (Q2 > mu2) ? log( log(Q2/lam2) / log(mu2/lam2) ) : 0.;
213 0 : double ds = sqrt(s);
214 0 : double s2 = s * s;
215 0 : double s3 = s2 * s;
216 :
217 : // uv :
218 0 : double nu = 2.284 + 0.802 * s + 0.055 * s2;
219 0 : double aku = 0.590 - 0.024 * s;
220 0 : double bku = 0.131 + 0.063 * s;
221 0 : double au = -0.449 - 0.138 * s - 0.076 * s2;
222 0 : double bu = 0.213 + 2.669 * s - 0.728 * s2;
223 0 : double cu = 8.854 - 9.135 * s + 1.979 * s2;
224 0 : double du = 2.997 + 0.753 * s - 0.076 * s2;
225 0 : double uv = grvv (x, nu, aku, bku, au, bu, cu, du);
226 :
227 : // dv :
228 0 : double nd = 0.371 + 0.083 * s + 0.039 * s2;
229 : double akd = 0.376;
230 0 : double bkd = 0.486 + 0.062 * s;
231 0 : double ad = -0.509 + 3.310 * s - 1.248 * s2;
232 0 : double bd = 12.41 - 10.52 * s + 2.267 * s2;
233 0 : double cd = 6.373 - 6.208 * s + 1.418 * s2;
234 0 : double dd = 3.691 + 0.799 * s - 0.071 * s2;
235 0 : double dv = grvv (x, nd, akd, bkd, ad, bd, cd, dd);
236 :
237 : // udb :
238 : double alx = 1.451;
239 : double bex = 0.271;
240 0 : double akx = 0.410 - 0.232 * s;
241 0 : double bkx = 0.534 - 0.457 * s;
242 0 : double agx = 0.890 - 0.140 * s;
243 : double bgx = -0.981;
244 0 : double cx = 0.320 + 0.683 * s;
245 0 : double dx = 4.752 + 1.164 * s + 0.286 * s2;
246 0 : double ex = 4.119 + 1.713 * s;
247 0 : double esx = 0.682 + 2.978 * s;
248 0 : double udb = grvw (x, s, alx, bex, akx, bkx, agx, bgx, cx,
249 : dx, ex, esx);
250 :
251 : // del :
252 0 : double ne = 0.082 + 0.014 * s + 0.008 * s2;
253 0 : double ake = 0.409 - 0.005 * s;
254 0 : double bke = 0.799 + 0.071 * s;
255 0 : double ae = -38.07 + 36.13 * s - 0.656 * s2;
256 0 : double be = 90.31 - 74.15 * s + 7.645 * s2;
257 : double ce = 0.;
258 0 : double de = 7.486 + 1.217 * s - 0.159 * s2;
259 0 : double del = grvv (x, ne, ake, bke, ae, be, ce, de);
260 :
261 : // sb :
262 : double sts = 0.;
263 : double als = 0.914;
264 : double bes = 0.577;
265 0 : double aks = 1.798 - 0.596 * s;
266 0 : double as = -5.548 + 3.669 * ds - 0.616 * s;
267 0 : double bs = 18.92 - 16.73 * ds + 5.168 * s;
268 0 : double dst = 6.379 - 0.350 * s + 0.142 * s2;
269 0 : double est = 3.981 + 1.638 * s;
270 : double ess = 6.402;
271 0 : double sb = grvs (x, s, sts, als, bes, aks, as, bs, dst, est, ess);
272 :
273 : // cb :
274 : double stc = 0.888;
275 : double alc = 1.01;
276 : double bec = 0.37;
277 : double akc = 0.;
278 : double ac = 0.;
279 0 : double bc = 4.24 - 0.804 * s;
280 0 : double dct = 3.46 - 1.076 * s;
281 0 : double ect = 4.61 + 1.49 * s;
282 0 : double esc = 2.555 + 1.961 * s;
283 0 : double chm = grvs (x, s, stc, alc, bec, akc, ac, bc, dct, ect, esc);
284 :
285 : // bb :
286 : double stb = 1.351;
287 : double alb = 1.00;
288 : double beb = 0.51;
289 : double akb = 0.;
290 : double ab = 0.;
291 : double bb = 1.848;
292 0 : double dbt = 2.929 + 1.396 * s;
293 0 : double ebt = 4.71 + 1.514 * s;
294 0 : double esb = 4.02 + 1.239 * s;
295 0 : double bot = grvs (x, s, stb, alb, beb, akb, ab, bb, dbt, ebt, esb);
296 :
297 : // gl :
298 : double alg = 0.524;
299 : double beg = 1.088;
300 0 : double akg = 1.742 - 0.930 * s;
301 0 : double bkg = - 0.399 * s2;
302 0 : double ag = 7.486 - 2.185 * s;
303 0 : double bg = 16.69 - 22.74 * s + 5.779 * s2;
304 0 : double cg = -25.59 + 29.71 * s - 7.296 * s2;
305 0 : double dg = 2.792 + 2.215 * s + 0.422 * s2 - 0.104 * s3;
306 0 : double eg = 0.807 + 2.005 * s;
307 0 : double esg = 3.841 + 0.316 * s;
308 0 : double gl = grvw (x, s, alg, beg, akg, bkg, ag, bg, cg,
309 : dg, eg, esg);
310 :
311 : // Update values
312 0 : xg = gl;
313 0 : xu = uv + 0.5*(udb - del);
314 0 : xd = dv + 0.5*(udb + del);
315 0 : xubar = 0.5*(udb - del);
316 0 : xdbar = 0.5*(udb + del);
317 0 : xs = sb;
318 0 : xsbar = sb;
319 0 : xc = chm;
320 0 : xb = bot;
321 :
322 : // Subdivision of valence and sea.
323 0 : xuVal = uv;
324 0 : xuSea = xubar;
325 0 : xdVal = dv;
326 0 : xdSea = xdbar;
327 :
328 : // idSav = 9 to indicate that all flavours reset.
329 0 : idSav = 9;
330 :
331 0 : }
332 :
333 : //--------------------------------------------------------------------------
334 :
335 : double GRV94L::grvv (double x, double n, double ak, double bk, double a,
336 : double b, double c, double d) {
337 :
338 0 : double dx = sqrt(x);
339 0 : return n * pow(x, ak) * (1. + a * pow(x, bk) + x * (b + c * dx)) *
340 0 : pow(1. - x, d);
341 :
342 : }
343 :
344 : //--------------------------------------------------------------------------
345 :
346 : double GRV94L::grvw (double x, double s, double al, double be, double ak,
347 : double bk, double a, double b, double c, double d, double e, double es) {
348 :
349 0 : double lx = log(1./x);
350 0 : return (pow(x, ak) * (a + x * (b + x * c)) * pow(lx, bk) + pow(s, al)
351 0 : * exp(-e + sqrt(es * pow(s, be) * lx))) * pow(1. - x, d);
352 :
353 : }
354 :
355 : //--------------------------------------------------------------------------
356 :
357 : double GRV94L::grvs (double x, double s, double sth, double al, double be,
358 : double ak, double ag, double b, double d, double e, double es) {
359 :
360 0 : if(s <= sth) {
361 0 : return 0.;
362 : } else {
363 0 : double dx = sqrt(x);
364 0 : double lx = log(1./x);
365 0 : return pow(s - sth, al) / pow(lx, ak) * (1. + ag * dx + b * x) *
366 0 : pow(1. - x, d) * exp(-e + sqrt(es * pow(s, be) * lx));
367 : }
368 :
369 0 : }
370 :
371 : //==========================================================================
372 :
373 : // Gives the CTEQ 5 L (leading order) parton distribution function set
374 : // in parametrized form. Parametrization by J. Pumplin.
375 : // Ref: CTEQ Collaboration, H.L. Lai et al., Eur.Phys.J. C12 (2000) 375.
376 :
377 : // The range of (x, Q) covered by this parametrization of the QCD
378 : // evolved parton distributions is 1E-6 < x < 1, 1.1 GeV < Q < 10 TeV.
379 : // In the current implementation, densities are frozen at borders.
380 :
381 : void CTEQ5L::xfUpdate(int , double x, double Q2) {
382 :
383 : // Constrain x and Q2 to range for which parametrization is valid.
384 0 : double Q = sqrt( max( 1., min( 1e8, Q2) ) );
385 0 : x = max( 1e-6, min( 1.-1e-10, x) );
386 :
387 : // Derived kinematical quantities.
388 0 : double y = - log(x);
389 0 : double u = log( x / 0.00001);
390 0 : double x1 = 1. - x;
391 0 : double x1L = log(1. - x);
392 : double sumUbarDbar = 0.;
393 :
394 : // Parameters of parametrizations.
395 : const double Qmin[8] = { 0., 0., 0., 0., 0., 0., 1.3, 4.5};
396 : const double alpha[8] = { 0.2987216, 0.3407552, 0.4491863, 0.2457668,
397 : 0.5293999, 0.3713141, 0.03712017, 0.004952010 };
398 : const double ut1[8] = { 4.971265, 2.612618, -0.4656819, 3.862583,
399 : 0.1895615, 3.753257, 4.400772, 5.562568 };
400 : const double ut2[8] = { -1.105128, -1.258304e5, -274.2390, -1.265969,
401 : -3.069097, -1.113085, -1.356116, -1.801317 };
402 : const double am[8][9][3] = {
403 : // d.
404 : { { 0.5292616E+01, -0.2751910E+01, -0.2488990E+01 },
405 : { 0.9714424E+00, 0.1011827E-01, -0.1023660E-01 },
406 : { -0.1651006E+02, 0.7959721E+01, 0.8810563E+01 },
407 : { -0.1643394E+02, 0.5892854E+01, 0.9348874E+01 },
408 : { 0.3067422E+02, 0.4235796E+01, -0.5112136E+00 },
409 : { 0.2352526E+02, -0.5305168E+01, -0.1169174E+02 },
410 : { -0.1095451E+02, 0.3006577E+01, 0.5638136E+01 },
411 : { -0.1172251E+02, -0.2183624E+01, 0.4955794E+01 },
412 : { 0.1662533E-01, 0.7622870E-02, -0.4895887E-03 } },
413 : // u.
414 : { { 0.9905300E+00, -0.4502235E+00, 0.1624441E+00 },
415 : { 0.8867534E+00, 0.1630829E-01, -0.4049085E-01 },
416 : { 0.8547974E+00, 0.3336301E+00, 0.1371388E+00 },
417 : { 0.2941113E+00, -0.1527905E+01, 0.2331879E+00 },
418 : { 0.3384235E+02, 0.3715315E+01, 0.8276930E+00 },
419 : { 0.6230115E+01, 0.3134639E+01, -0.1729099E+01 },
420 : { -0.1186928E+01, -0.3282460E+00, 0.1052020E+00 },
421 : { -0.8545702E+01, -0.6247947E+01, 0.3692561E+01 },
422 : { 0.1724598E-01, 0.7120465E-02, 0.4003646E-04 } },
423 : // g.
424 : { { 0.1193572E+03, -0.3886845E+01, -0.1133965E+01 },
425 : { -0.9421449E+02, 0.3995885E+01, 0.1607363E+01 },
426 : { 0.4206383E+01, 0.2485954E+00, 0.2497468E+00 },
427 : { 0.1210557E+03, -0.3015765E+01, -0.1423651E+01 },
428 : { -0.1013897E+03, -0.7113478E+00, 0.2621865E+00 },
429 : { -0.1312404E+01, -0.9297691E+00, -0.1562531E+00 },
430 : { 0.1627137E+01, 0.4954111E+00, -0.6387009E+00 },
431 : { 0.1537698E+00, -0.2487878E+00, 0.8305947E+00 },
432 : { 0.2496448E-01, 0.2457823E-02, 0.8234276E-03 } },
433 : // ubar + dbar.
434 : { { 0.2647441E+02, 0.1059277E+02, -0.9176654E+00 },
435 : { 0.1990636E+01, 0.8558918E-01, 0.4248667E-01 },
436 : { -0.1476095E+02, -0.3276255E+02, 0.1558110E+01 },
437 : { -0.2966889E+01, -0.3649037E+02, 0.1195914E+01 },
438 : { -0.1000519E+03, -0.2464635E+01, 0.1964849E+00 },
439 : { 0.3718331E+02, 0.4700389E+02, -0.2772142E+01 },
440 : { -0.1872722E+02, -0.2291189E+02, 0.1089052E+01 },
441 : { -0.1628146E+02, -0.1823993E+02, 0.2537369E+01 },
442 : { -0.1156300E+01, -0.1280495E+00, 0.5153245E-01 } },
443 : // dbar/ubar.
444 : { { -0.6556775E+00, 0.2490190E+00, 0.3966485E-01 },
445 : { 0.1305102E+01, -0.1188925E+00, -0.4600870E-02 },
446 : { -0.2371436E+01, 0.3566814E+00, -0.2834683E+00 },
447 : { -0.6152826E+01, 0.8339877E+00, -0.7233230E+00 },
448 : { -0.8346558E+01, 0.2892168E+01, 0.2137099E+00 },
449 : { 0.1279530E+02, 0.1021114E+00, 0.5787439E+00 },
450 : { 0.5858816E+00, -0.1940375E+01, -0.4029269E+00 },
451 : { -0.2795725E+02, -0.5263392E+00, 0.1290229E+01 },
452 : { 0.0000000E+00, 0.0000000E+00, 0.0000000E+00 } },
453 : // sbar.
454 : { { 0.1580931E+01, -0.2273826E+01, -0.1822245E+01 },
455 : { 0.2702644E+01, 0.6763243E+00, 0.7231586E-02 },
456 : { -0.1857924E+02, 0.3907500E+01, 0.5850109E+01 },
457 : { -0.3044793E+02, 0.2639332E+01, 0.5566644E+01 },
458 : { -0.4258011E+01, -0.5429244E+01, 0.4418946E+00 },
459 : { 0.3465259E+02, -0.5532604E+01, -0.4904153E+01 },
460 : { -0.1658858E+02, 0.2923275E+01, 0.2266286E+01 },
461 : { -0.1149263E+02, 0.2877475E+01, -0.7999105E+00 },
462 : { 0.0000000E+00, 0.0000000E+00, 0.0000000E+00 } },
463 : // cbar.
464 : { { -0.8293661E+00, -0.3982375E+01, -0.6494283E-01 },
465 : { 0.2754618E+01, 0.8338636E+00, -0.6885160E-01 },
466 : { -0.1657987E+02, 0.1439143E+02, -0.6887240E+00 },
467 : { -0.2800703E+02, 0.1535966E+02, -0.7377693E+00 },
468 : { -0.6460216E+01, -0.4783019E+01, 0.4913297E+00 },
469 : { 0.3141830E+02, -0.3178031E+02, 0.7136013E+01 },
470 : { -0.1802509E+02, 0.1862163E+02, -0.4632843E+01 },
471 : { -0.1240412E+02, 0.2565386E+02, -0.1066570E+02 },
472 : { 0.0000000E+00, 0.0000000E+00, 0.0000000E+00 } },
473 : // bbar.
474 : { { -0.6031237E+01, 0.1992727E+01, -0.1076331E+01 },
475 : { 0.2933912E+01, 0.5839674E+00, 0.7509435E-01 },
476 : { -0.8284919E+01, 0.1488593E+01, -0.8251678E+00 },
477 : { -0.1925986E+02, 0.2805753E+01, -0.3015446E+01 },
478 : { -0.9480483E+01, -0.9767837E+00, -0.1165544E+01 },
479 : { 0.2193195E+02, -0.1788518E+02, 0.9460908E+01 },
480 : { -0.1327377E+02, 0.1201754E+02, -0.6277844E+01 },
481 : { 0.0000000E+00, 0.0000000E+00, 0.0000000E+00 },
482 : { 0.0000000E+00, 0.0000000E+00, 0.0000000E+00 } } };
483 :
484 : // Loop over 8 different parametrizations. Check if inside allowed region.
485 0 : for (int i = 0; i < 8; ++i) {
486 : double answer = 0.;
487 0 : if (Q > max(Qmin[i], alpha[i])) {
488 :
489 : // Evaluate answer.
490 0 : double tmp = log(Q / alpha[i]);
491 0 : double sb = log(tmp);
492 0 : double sb1 = sb - 1.2;
493 0 : double sb2 = sb1*sb1;
494 0 : double af[9];
495 0 : for (int j = 0; j < 9; ++j)
496 0 : af[j] = am[i][j][0] + sb1 * am[i][j][1] + sb2 * am[i][j][2];
497 0 : double part1 = af[1] * pow( y, 1. + 0.01 * af[4]) * (1. + af[8] * u);
498 0 : double part2 = af[0] * x1 + af[3] * x;
499 0 : double part3 = x * x1 * (af[5] + af[6] * x1 + af[7] * x * x1);
500 0 : double part4 = (ut2[i] < -100.) ? ut1[i] * x1L + af[2] * x1L
501 0 : : ut1[i] * x1L + af[2] * log(x1 + exp(ut2[i]));
502 0 : answer = x * exp( part1 + part2 + part3 + part4);
503 0 : answer *= 1. - Qmin[i] / Q;
504 0 : }
505 :
506 : // Store results.
507 0 : if (i == 0) xd = x * answer;
508 0 : else if (i == 1) xu = x * answer;
509 0 : else if (i == 2) xg = x * answer;
510 0 : else if (i == 3) sumUbarDbar = x * answer;
511 0 : else if (i == 4) { xubar = sumUbarDbar / (1. + answer);
512 0 : xdbar = sumUbarDbar * answer / (1. + answer); }
513 0 : else if (i == 5) {xs = x * answer; xsbar = xs;}
514 0 : else if (i == 6) xc = x * answer;
515 0 : else if (i == 7) xb = x * answer;
516 : }
517 :
518 : // Subdivision of valence and sea.
519 0 : xuVal = xu - xubar;
520 0 : xuSea = xubar;
521 0 : xdVal = xd - xdbar;
522 0 : xdSea = xdbar;
523 :
524 : // idSav = 9 to indicate that all flavours reset.
525 0 : idSav = 9;
526 :
527 0 : }
528 :
529 : //==========================================================================
530 :
531 : // The MSTWpdf class.
532 : // MSTW 2008 PDF's, specifically the LO one.
533 : // Original C++ version by Jeppe Andersen.
534 : // Modified by Graeme Watt <watt(at)hep.ucl.ac.uk>.
535 :
536 : //--------------------------------------------------------------------------
537 :
538 : // Constants: could be changed here if desired, but normally should not.
539 : // These are of technical nature, as described for each.
540 :
541 : // Number of parton flavours, x and Q2 grid points,
542 : // bins below c and b thresholds.
543 : const int MSTWpdf::np = 12;
544 : const int MSTWpdf::nx = 64;
545 : const int MSTWpdf::nq = 48;
546 : const int MSTWpdf::nqc0 = 4;
547 : const int MSTWpdf::nqb0 = 14;
548 :
549 : // Range of (x, Q2) grid.
550 : const double MSTWpdf::xmin = 1e-6;
551 : const double MSTWpdf::xmax = 1.0;
552 : const double MSTWpdf::qsqmin = 1.0;
553 : const double MSTWpdf::qsqmax = 1e9;
554 :
555 : // Array of x values.
556 : const double MSTWpdf::xxInit[65] = {0., 1e-6, 2e-6, 4e-6, 6e-6, 8e-6,
557 : 1e-5, 2e-5, 4e-5, 6e-5, 8e-5, 1e-4, 2e-4, 4e-4, 6e-4, 8e-4,
558 : 1e-3, 2e-3, 4e-3, 6e-3, 8e-3, 1e-2, 1.4e-2, 2e-2, 3e-2, 4e-2, 6e-2,
559 : 8e-2, 0.10, 0.125, 0.15, 0.175, 0.20, 0.225, 0.25, 0.275, 0.30,
560 : 0.325, 0.35, 0.375, 0.40, 0.425, 0.45, 0.475, 0.50, 0.525, 0.55,
561 : 0.575, 0.60, 0.625, 0.65, 0.675, 0.70, 0.725, 0.75, 0.775, 0.80,
562 : 0.825, 0.85, 0.875, 0.90, 0.925, 0.95, 0.975, 1.0 };
563 :
564 : // Array of Q values.
565 : const double MSTWpdf::qqInit[49] = {0., 1.0, 1.25, 1.5, 0., 0., 2.5, 3.2,
566 : 4.0, 5.0, 6.4, 8.0, 10., 12., 0., 0., 26.0, 40.0, 64.0, 1e2, 1.6e2,
567 : 2.4e2, 4e2, 6.4e2, 1e3, 1.8e3, 3.2e3, 5.6e3, 1e4, 1.8e4, 3.2e4, 5.6e4,
568 : 1e5, 1.8e5, 3.2e5, 5.6e5, 1e6, 1.8e6, 3.2e6, 5.6e6, 1e7, 1.8e7, 3.2e7,
569 : 5.6e7, 1e8, 1.8e8, 3.2e8, 5.6e8, 1e9 };
570 :
571 : //--------------------------------------------------------------------------
572 :
573 : // Initialize PDF: read in data grid from file and set up interpolation.
574 :
575 : void MSTWpdf::init(int iFitIn, string xmlPath, Info* infoPtr) {
576 :
577 : // Choice of fit among possibilities. Counters and temporary variables.
578 0 : iFit = iFitIn;
579 : int i,n,m,k,l,j;
580 0 : double dtemp;
581 :
582 : // Variables used for initialising c_ij array:
583 0 : double f[np+1][nx+1][nq+1];
584 0 : double f1[np+1][nx+1][nq+1]; // derivative w.r.t. x
585 0 : double f2[np+1][nx+1][nq+1]; // derivative w.r.t. q
586 0 : double f12[np+1][nx+1][nq+1];// cross derivative
587 0 : double f21[np+1][nx+1][nq+1];// cross derivative
588 0 : int wt[16][16]={{1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
589 : {0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0},
590 : {-3,0,0,3,0,0,0,0,-2,0,0,-1,0,0,0,0},
591 : {2,0,0,-2,0,0,0,0,1,0,0,1,0,0,0,0},
592 : {0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0},
593 : {0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0},
594 : {0,0,0,0,-3,0,0,3,0,0,0,0,-2,0,0,-1},
595 : {0,0,0,0,2,0,0,-2,0,0,0,0,1,0,0,1},
596 : {-3,3,0,0,-2,-1,0,0,0,0,0,0,0,0,0,0},
597 : {0,0,0,0,0,0,0,0,-3,3,0,0,-2,-1,0,0},
598 : {9,-9,9,-9,6,3,-3,-6,6,-6,-3,3,4,2,1,2},
599 : {-6,6,-6,6,-4,-2,2,4,-3,3,3,-3,-2,-1,-1,-2},
600 : {2,-2,0,0,1,1,0,0,0,0,0,0,0,0,0,0},
601 : {0,0,0,0,0,0,0,0,2,-2,0,0,1,1,0,0},
602 : {-6,6,-6,6,-3,-3,3,3,-4,4,2,-2,-2,-2,-1,-1},
603 : {4,-4,4,-4,2,2,-2,-2,2,-2,-2,2,1,1,1,1}};
604 0 : double xxd,d1d2,cl[16],x[16],d1,d2,y[5],y1[5],y2[5],y12[5];
605 : double mc2,mb2,eps=1e-6; // q^2 grid points at mc2+eps, mb2+eps
606 :
607 : // Select which data file to read for current fit.
608 0 : if (xmlPath[ xmlPath.length() - 1 ] != '/') xmlPath += "/";
609 0 : string fileName = " ";
610 0 : if (iFit == 1) fileName = "mrstlostar.00.dat";
611 0 : if (iFit == 2) fileName = "mrstlostarstar.00.dat";
612 0 : if (iFit == 3) fileName = "mstw2008lo.00.dat";
613 0 : if (iFit == 4) fileName = "mstw2008nlo.00.dat";
614 :
615 : // Open data file.
616 0 : ifstream data_file( (xmlPath + fileName).c_str() );
617 0 : if (!data_file.good()) {
618 0 : if (infoPtr != 0) infoPtr->errorMsg("Error from MSTWpdf::init: "
619 0 : "did not find parametrization file ", fileName);
620 0 : else cout << " Error from MSTWpdf::init: "
621 0 : << "did not find parametrization file " << fileName << endl;
622 0 : isSet = false;
623 0 : return;
624 : }
625 :
626 : // Read distance, tolerance, heavy quark masses
627 : // and alphaS values from file.
628 0 : char comma;
629 0 : int nExtraFlavours;
630 0 : data_file.ignore(256,'\n');
631 0 : data_file.ignore(256,'\n');
632 0 : data_file.ignore(256,'='); data_file >> distance >> tolerance;
633 0 : data_file.ignore(256,'='); data_file >> mCharm;
634 0 : data_file.ignore(256,'='); data_file >> mBottom;
635 0 : data_file.ignore(256,'='); data_file >> alphaSQ0;
636 0 : data_file.ignore(256,'='); data_file >> alphaSMZ;
637 0 : data_file.ignore(256,'='); data_file >> alphaSorder >> comma >> alphaSnfmax;
638 0 : data_file.ignore(256,'='); data_file >> nExtraFlavours;
639 0 : data_file.ignore(256,'\n');
640 0 : data_file.ignore(256,'\n');
641 0 : data_file.ignore(256,'\n');
642 :
643 : // Use c and b quark masses for outlay of qq array.
644 0 : for (int iqq = 0; iqq < 49; ++iqq) qq[iqq] = qqInit[iqq];
645 0 : mc2=mCharm*mCharm;
646 0 : mb2=mBottom*mBottom;
647 0 : qq[4]=mc2;
648 0 : qq[5]=mc2+eps;
649 0 : qq[14]=mb2;
650 0 : qq[15]=mb2+eps;
651 :
652 : // Check that the heavy quark masses are sensible.
653 0 : if (mc2 < qq[3] || mc2 > qq[6]) {
654 0 : if (infoPtr != 0) infoPtr->errorMsg("Error from MSTWpdf::init: "
655 : "invalid mCharm");
656 0 : else cout << " Error from MSTWpdf::init: invalid mCharm" << endl;
657 0 : isSet = false;
658 0 : return;
659 : }
660 0 : if (mb2 < qq[13] || mb2 > qq[16]) {
661 0 : if (infoPtr != 0) infoPtr->errorMsg("Error from MSTWpdf::init: "
662 : "invalid mBottom");
663 0 : else cout << " Error from MSTWpdf::init: invalid mBottom" << endl;
664 0 : isSet = false;
665 0 : return;
666 : }
667 :
668 : // The nExtraFlavours variable is provided to aid compatibility
669 : // with future grids where, for example, a photon distribution
670 : // might be provided (cf. the MRST2004QED PDFs).
671 0 : if (nExtraFlavours < 0 || nExtraFlavours > 1) {
672 0 : if (infoPtr != 0) infoPtr->errorMsg("Error from MSTWpdf::init: "
673 : "invalid nExtraFlavours");
674 0 : else cout << " Error from MSTWpdf::init: invalid nExtraFlavours" << endl;
675 0 : isSet = false;
676 0 : return;
677 : }
678 :
679 : // Now read in the grids from the grid file.
680 0 : for (n=1;n<=nx-1;n++)
681 0 : for (m=1;m<=nq;m++) {
682 0 : for (i=1;i<=9;i++)
683 0 : data_file >> f[i][n][m];
684 0 : if (alphaSorder==2) { // only at NNLO
685 0 : data_file >> f[10][n][m]; // = chm-cbar
686 0 : data_file >> f[11][n][m]; // = bot-bbar
687 : }
688 : else {
689 0 : f[10][n][m] = 0.; // = chm-cbar
690 0 : f[11][n][m] = 0.; // = bot-bbar
691 : }
692 0 : if (nExtraFlavours>0)
693 0 : data_file >> f[12][n][m]; // = photon
694 : else
695 0 : f[12][n][m] = 0.; // photon
696 0 : if (data_file.eof()) {
697 0 : if (infoPtr != 0) infoPtr->errorMsg("Error from MSTWpdf::init: "
698 : "failed to read in data file");
699 0 : else cout << " Error from MSTWpdf::init: failed to read in data file"
700 0 : << endl;
701 0 : isSet = false;
702 0 : return;
703 : }
704 : }
705 :
706 : // Check that ALL the file contents have been read in.
707 0 : data_file >> dtemp;
708 0 : if (!data_file.eof()) {
709 0 : if (infoPtr != 0) infoPtr->errorMsg("Error from MSTWpdf::init: "
710 : "failed to read in data file");
711 0 : else cout << " Error from MSTWpdf::init: failed to read in data file"
712 0 : << endl;
713 0 : isSet = false;
714 0 : return;
715 : }
716 :
717 : // Close the datafile.
718 0 : data_file.close();
719 :
720 : // PDFs are identically zero at x = 1.
721 0 : for (i=1;i<=np;i++)
722 0 : for (m=1;m<=nq;m++)
723 0 : f[i][nx][m]=0.0;
724 :
725 : // Set up the new array in log10(x) and log10(qsq).
726 0 : for (i=1;i<=nx;i++)
727 0 : xx[i]=log10(xxInit[i]);
728 0 : for (m=1;m<=nq;m++)
729 0 : qq[m]=log10(qq[m]);
730 :
731 : // Now calculate the derivatives used for bicubic interpolation.
732 0 : for (i=1;i<=np;i++) {
733 :
734 : // Start by calculating the first x derivatives
735 : // along the first x value:
736 0 : for (m=1;m<=nq;m++) {
737 0 : f1[i][1][m]=polderivative1(xx[1],xx[2],xx[3],f[i][1][m],f[i][2][m],
738 0 : f[i][3][m]);
739 : // Then along the rest (up to the last):
740 0 : for (k=2;k<nx;k++)
741 0 : f1[i][k][m]=polderivative2(xx[k-1],xx[k],xx[k+1],f[i][k-1][m],
742 0 : f[i][k][m],f[i][k+1][m]);
743 : // Then for the last column:
744 0 : f1[i][nx][m]=polderivative3(xx[nx-2],xx[nx-1],xx[nx],f[i][nx-2][m],
745 0 : f[i][nx-1][m],f[i][nx][m]);
746 : }
747 :
748 : // Then calculate the qq derivatives. At NNLO there are
749 : // discontinuities in the PDFs at mc2 and mb2, so calculate
750 : // the derivatives at q^2 = mc2, mc2+eps, mb2, mb2+eps in
751 : // the same way as at the endpoints qsqmin and qsqmax.
752 0 : for (m=1;m<=nq;m++) {
753 0 : if (m==1 || m==nqc0+1 || m==nqb0+1) {
754 0 : for (k=1;k<=nx;k++)
755 0 : f2[i][k][m]=polderivative1(qq[m],qq[m+1],qq[m+2],
756 0 : f[i][k][m],f[i][k][m+1],f[i][k][m+2]);
757 : }
758 0 : else if (m==nq || m==nqc0 || m==nqb0) {
759 0 : for (k=1;k<=nx;k++)
760 0 : f2[i][k][m]=polderivative3(qq[m-2],qq[m-1],qq[m],
761 0 : f[i][k][m-2],f[i][k][m-1],f[i][k][m]);
762 : }
763 : else {
764 : // The rest:
765 0 : for (k=1;k<=nx;k++)
766 0 : f2[i][k][m]=polderivative2(qq[m-1],qq[m],qq[m+1],
767 0 : f[i][k][m-1],f[i][k][m],f[i][k][m+1]);
768 : }
769 : }
770 :
771 : // Now, calculate the cross derivatives.
772 : // Calculate these as the average between (d/dx)(d/dy) and (d/dy)(d/dx).
773 :
774 : // First calculate (d/dx)(d/dy).
775 : // Start by calculating the first x derivatives
776 : // along the first x value:
777 0 : for (m=1;m<=nq;m++)
778 0 : f12[i][1][m]=polderivative1(xx[1],xx[2],xx[3],f2[i][1][m],
779 0 : f2[i][2][m],f2[i][3][m]);
780 : // Then along the rest (up to the last):
781 0 : for (k=2;k<nx;k++) {
782 0 : for (m=1;m<=nq;m++)
783 0 : f12[i][k][m]=polderivative2(xx[k-1],xx[k],xx[k+1],f2[i][k-1][m],
784 0 : f2[i][k][m],f2[i][k+1][m]);
785 : }
786 : // Then for the last column:
787 0 : for (m=1;m<=nq;m++)
788 0 : f12[i][nx][m]=polderivative3(xx[nx-2],xx[nx-1],xx[nx],
789 0 : f2[i][nx-2][m],f2[i][nx-1][m],f2[i][nx][m]);
790 :
791 : // Now calculate (d/dy)(d/dx).
792 0 : for (m=1;m<=nq;m++) {
793 0 : if (m==1 || m==nqc0+1 || m==nqb0+1) {
794 0 : for (k=1;k<=nx;k++)
795 0 : f21[i][k][m]=polderivative1(qq[m],qq[m+1],qq[m+2],
796 0 : f1[i][k][m],f1[i][k][m+1],f1[i][k][m+2]);
797 : }
798 0 : else if (m==nq || m==nqc0 || m==nqb0) {
799 0 : for (k=1;k<=nx;k++)
800 0 : f21[i][k][m]=polderivative3(qq[m-2],qq[m-1],qq[m],
801 0 : f1[i][k][m-2],f1[i][k][m-1],f1[i][k][m]);
802 : }
803 : else {
804 : // The rest:
805 0 : for (k=1;k<=nx;k++)
806 0 : f21[i][k][m]=polderivative2(qq[m-1],qq[m],qq[m+1],
807 0 : f1[i][k][m-1],f1[i][k][m],f1[i][k][m+1]);
808 : }
809 : }
810 :
811 : // Now take the average of (d/dx)(d/dy) and (d/dy)(d/dx).
812 0 : for (k=1;k<=nx;k++) {
813 0 : for (m=1;m<=nq;m++) {
814 0 : f12[i][k][m] = 0.5*(f12[i][k][m]+f21[i][k][m]);
815 : }
816 : }
817 :
818 : // Now calculate the coefficients c_ij.
819 0 : for (n=1;n<=nx-1;n++) {
820 0 : for (m=1;m<=nq-1;m++) {
821 0 : d1=xx[n+1]-xx[n];
822 0 : d2=qq[m+1]-qq[m];
823 0 : d1d2=d1*d2;
824 :
825 0 : y[1]=f[i][n][m];
826 0 : y[2]=f[i][n+1][m];
827 0 : y[3]=f[i][n+1][m+1];
828 0 : y[4]=f[i][n][m+1];
829 :
830 0 : y1[1]=f1[i][n][m];
831 0 : y1[2]=f1[i][n+1][m];
832 0 : y1[3]=f1[i][n+1][m+1];
833 0 : y1[4]=f1[i][n][m+1];
834 :
835 0 : y2[1]=f2[i][n][m];
836 0 : y2[2]=f2[i][n+1][m];
837 0 : y2[3]=f2[i][n+1][m+1];
838 0 : y2[4]=f2[i][n][m+1];
839 :
840 0 : y12[1]=f12[i][n][m];
841 0 : y12[2]=f12[i][n+1][m];
842 0 : y12[3]=f12[i][n+1][m+1];
843 0 : y12[4]=f12[i][n][m+1];
844 :
845 0 : for (k=1;k<=4;k++) {
846 0 : x[k-1]=y[k];
847 0 : x[k+3]=y1[k]*d1;
848 0 : x[k+7]=y2[k]*d2;
849 0 : x[k+11]=y12[k]*d1d2;
850 : }
851 :
852 0 : for (l=0;l<=15;l++) {
853 : xxd=0.0;
854 0 : for (k=0;k<=15;k++) xxd+= wt[l][k]*x[k];
855 0 : cl[l]=xxd;
856 : }
857 :
858 : l=0;
859 0 : for (k=1;k<=4;k++)
860 0 : for (j=1;j<=4;j++) c[i][n][m][k][j]=cl[l++];
861 : } //m
862 : } //n
863 : } // i
864 :
865 0 : }
866 :
867 : //--------------------------------------------------------------------------
868 :
869 : // Update PDF values.
870 :
871 : void MSTWpdf::xfUpdate(int , double x, double Q2) {
872 :
873 : // Update using MSTW routine.
874 0 : double q = sqrtpos(Q2);
875 : // Quarks:
876 0 : double dn = parton(1,x,q);
877 0 : double up = parton(2,x,q);
878 0 : double str = parton(3,x,q);
879 0 : double chm = parton(4,x,q);
880 0 : double bot = parton(5,x,q);
881 : // Valence quarks:
882 0 : double dnv = parton(7,x,q);
883 0 : double upv = parton(8,x,q);
884 0 : double sv = parton(9,x,q);
885 0 : double cv = parton(10,x,q);
886 0 : double bv = parton(11,x,q);
887 : // Antiquarks = quarks - valence quarks:
888 0 : double dsea = dn - dnv;
889 0 : double usea = up - upv;
890 0 : double sbar = str - sv;
891 0 : double cbar = chm - cv;
892 0 : double bbar = bot - bv;
893 : // Gluon:
894 0 : double glu = parton(0,x,q);
895 : // Photon (= zero unless considering QED contributions):
896 0 : double phot = parton(13,x,q);
897 :
898 : // Transfer to Pythia notation.
899 0 : xg = glu;
900 0 : xu = up;
901 0 : xd = dn;
902 0 : xubar = usea;
903 0 : xdbar = dsea;
904 0 : xs = str;
905 0 : xsbar = sbar;
906 0 : xc = 0.5 * (chm + cbar);
907 0 : xb = 0.5 * (bot + bbar);
908 0 : xgamma = phot;
909 :
910 : // Subdivision of valence and sea.
911 0 : xuVal = upv;
912 0 : xuSea = xubar;
913 0 : xdVal = dnv;
914 0 : xdSea = xdbar;
915 :
916 : // idSav = 9 to indicate that all flavours reset.
917 0 : idSav = 9;
918 :
919 0 : }
920 :
921 : //--------------------------------------------------------------------------
922 :
923 : // Returns the PDF value for parton of flavour 'f' at x,q.
924 :
925 : double MSTWpdf::parton(int f,double x,double q) {
926 :
927 : double qsq;
928 : int ip;
929 : int interpolate(1);
930 : double parton_pdf=0,parton_pdf1=0,anom;
931 : double xxx,qqq;
932 :
933 0 : qsq=q*q;
934 :
935 : // If mc2 < qsq < mc2+eps, then qsq = mc2+eps.
936 0 : if (qsq>pow(10.,qq[nqc0]) && qsq<pow(10.,qq[nqc0+1])) {
937 : qsq = pow(10.,qq[nqc0+1]);
938 0 : }
939 :
940 : // If mb2 < qsq < mb2+eps, then qsq = mb2+eps.
941 0 : if (qsq>pow(10.,qq[nqb0]) && qsq<pow(10.,qq[nqb0+1])) {
942 : qsq = pow(10.,qq[nqb0+1]);
943 0 : }
944 :
945 0 : if (x<xmin) {
946 : interpolate=0;
947 0 : if (x<=0.) return 0.;
948 : }
949 0 : else if (x>xmax) return 0.;
950 :
951 0 : if (qsq<qsqmin) {
952 : interpolate=-1;
953 0 : if (q<=0.) return 0.;
954 : }
955 0 : else if (qsq>qsqmax) {
956 : interpolate=0;
957 0 : }
958 :
959 0 : if (f==0) ip=1;
960 0 : else if (f>=1 && f<=5) ip=f+1;
961 0 : else if (f<=-1 && f>=-5) ip=-f+1;
962 0 : else if (f>=7 && f<=11) ip=f;
963 0 : else if (f==13) ip=12;
964 0 : else if (abs(f)==6 || f==12) return 0.;
965 : else return 0.;
966 :
967 : // Interpolation in log10(x), log10(qsq):
968 0 : xxx=log10(x);
969 0 : qqq=log10(qsq);
970 :
971 0 : if (interpolate==1) { // do usual interpolation
972 0 : parton_pdf=parton_interpolate(ip,xxx,qqq);
973 0 : if (f<=-1 && f>=-5) // antiquark = quark - valence
974 0 : parton_pdf -= parton_interpolate(ip+5,xxx,qqq);
975 : }
976 0 : else if (interpolate==-1) { // extrapolate to low Q^2
977 :
978 0 : if (x<xmin) { // extrapolate to low x
979 0 : parton_pdf = parton_extrapolate(ip,xxx,log10(qsqmin));
980 0 : parton_pdf1 = parton_extrapolate(ip,xxx,log10(1.01*qsqmin));
981 0 : if (f<=-1 && f>=-5) { // antiquark = quark - valence
982 0 : parton_pdf -= parton_extrapolate(ip+5,xxx,log10(qsqmin));
983 0 : parton_pdf1 -= parton_extrapolate(ip+5,xxx,log10(1.01*qsqmin));
984 0 : }
985 : }
986 : else { // do usual interpolation
987 0 : parton_pdf = parton_interpolate(ip,xxx,log10(qsqmin));
988 0 : parton_pdf1 = parton_interpolate(ip,xxx,log10(1.01*qsqmin));
989 0 : if (f<=-1 && f>=-5) { // antiquark = quark - valence
990 0 : parton_pdf -= parton_interpolate(ip+5,xxx,log10(qsqmin));
991 0 : parton_pdf1 -= parton_interpolate(ip+5,xxx,log10(1.01*qsqmin));
992 0 : }
993 : }
994 : // Calculate the anomalous dimension, dlog(xf)/dlog(qsq),
995 : // evaluated at qsqmin. Then extrapolate the PDFs to low
996 : // qsq < qsqmin by interpolating the anomalous dimenion between
997 : // the value at qsqmin and a value of 1 for qsq << qsqmin.
998 : // If value of PDF at qsqmin is very small, just set
999 : // anomalous dimension to 1 to prevent rounding errors.
1000 0 : if (fabs(parton_pdf) >= 1.e-5)
1001 0 : anom = max(-2.5, (parton_pdf1-parton_pdf)/parton_pdf/0.01);
1002 : else anom = 1.;
1003 0 : parton_pdf = parton_pdf*pow(qsq/qsqmin,anom*qsq/qsqmin+1.-qsq/qsqmin);
1004 :
1005 0 : }
1006 : else { // extrapolate outside PDF grid to low x or high Q^2
1007 0 : parton_pdf = parton_extrapolate(ip,xxx,qqq);
1008 0 : if (f<=-1 && f>=-5) // antiquark = quark - valence
1009 0 : parton_pdf -= parton_extrapolate(ip+5,xxx,qqq);
1010 : }
1011 :
1012 0 : return parton_pdf;
1013 0 : }
1014 :
1015 : //--------------------------------------------------------------------------
1016 :
1017 : // Interpolate PDF value inside data grid.
1018 :
1019 : double MSTWpdf::parton_interpolate(int ip, double xxx, double qqq) {
1020 :
1021 : double g, t, u;
1022 : int n, m, l;
1023 :
1024 0 : n=locate(xx,nx,xxx); // 0: below xmin, nx: above xmax
1025 0 : m=locate(qq,nq,qqq); // 0: below qsqmin, nq: above qsqmax
1026 :
1027 0 : t=(xxx-xx[n])/(xx[n+1]-xx[n]);
1028 0 : u=(qqq-qq[m])/(qq[m+1]-qq[m]);
1029 :
1030 : // Assume PDF proportional to (1-x)^p as x -> 1.
1031 0 : if (n==nx-1) {
1032 0 : double g0=((c[ip][n][m][1][4]*u+c[ip][n][m][1][3])*u
1033 0 : +c[ip][n][m][1][2])*u+c[ip][n][m][1][1]; // value at xx[n]
1034 0 : double g1=((c[ip][n-1][m][1][4]*u+c[ip][n-1][m][1][3])*u
1035 0 : +c[ip][n-1][m][1][2])*u+c[ip][n-1][m][1][1]; // value at xx[n-1]
1036 : double p = 1.0;
1037 0 : if (g0>0.0&&g1>0.0) p = log(g1/g0)/log((xx[n+1]-xx[n-1])/(xx[n+1]-xx[n]));
1038 0 : if (p<=1.0) p=1.0;
1039 0 : g=g0*pow((xx[n+1]-xxx)/(xx[n+1]-xx[n]),p);
1040 0 : }
1041 :
1042 : // Usual interpolation.
1043 : else {
1044 : g=0.0;
1045 0 : for (l=4;l>=1;l--) {
1046 0 : g=t*g+((c[ip][n][m][l][4]*u+c[ip][n][m][l][3])*u
1047 0 : +c[ip][n][m][l][2])*u+c[ip][n][m][l][1];
1048 : }
1049 : }
1050 :
1051 0 : return g;
1052 : }
1053 :
1054 : //--------------------------------------------------------------------------
1055 :
1056 : // Extrapolate PDF value outside data grid.
1057 :
1058 :
1059 : double MSTWpdf::parton_extrapolate(int ip, double xxx, double qqq) {
1060 :
1061 : double parton_pdf=0.;
1062 : int n,m;
1063 :
1064 0 : n=locate(xx,nx,xxx); // 0: below xmin, nx: above xmax
1065 0 : m=locate(qq,nq,qqq); // 0: below qsqmin, nq: above qsqmax
1066 :
1067 0 : if (n==0&&(m>0&&m<nq)) { // if extrapolation in small x only
1068 :
1069 : double f0,f1;
1070 0 : f0=parton_interpolate(ip,xx[1],qqq);
1071 0 : f1=parton_interpolate(ip,xx[2],qqq);
1072 0 : if ( f0>1e-3 && f1>1e-3 ) { // if values are positive, keep them so
1073 0 : f0=log(f0);
1074 0 : f1=log(f1);
1075 0 : parton_pdf=exp(f0+(f1-f0)/(xx[2]-xx[1])*(xxx-xx[1]));
1076 0 : } else // otherwise just extrapolate in the value
1077 0 : parton_pdf=f0+(f1-f0)/(xx[2]-xx[1])*(xxx-xx[1]);
1078 :
1079 0 : } if (n>0&&m==nq) { // if extrapolation into large q only
1080 :
1081 : double f0,f1;
1082 0 : f0=parton_interpolate(ip,xxx,qq[nq]);
1083 0 : f1=parton_interpolate(ip,xxx,qq[nq-1]);
1084 0 : if ( f0>1e-3 && f1>1e-3 ) { // if values are positive, keep them so
1085 0 : f0=log(f0);
1086 0 : f1=log(f1);
1087 0 : parton_pdf=exp(f0+(f0-f1)/(qq[nq]-qq[nq-1])*(qqq-qq[nq]));
1088 0 : } else // otherwise just extrapolate in the value
1089 0 : parton_pdf=f0+(f0-f1)/(qq[nq]-qq[nq-1])*(qqq-qq[nq]);
1090 :
1091 0 : } if (n==0&&m==nq) { // if extrapolation into large q AND small x
1092 :
1093 : double f0,f1;
1094 0 : f0=parton_extrapolate(ip,xx[1],qqq);
1095 0 : f1=parton_extrapolate(ip,xx[2],qqq);
1096 0 : if ( f0>1e-3 && f1>1e-3 ) { // if values are positive, keep them so
1097 0 : f0=log(f0);
1098 0 : f1=log(f1);
1099 0 : parton_pdf=exp(f0+(f1-f0)/(xx[2]-xx[1])*(xxx-xx[1]));
1100 0 : } else // otherwise just extrapolate in the value
1101 0 : parton_pdf=f0+(f1-f0)/(xx[2]-xx[1])*(xxx-xx[1]);
1102 :
1103 0 : }
1104 :
1105 0 : return parton_pdf;
1106 : }
1107 :
1108 : //--------------------------------------------------------------------------
1109 :
1110 : // Returns an integer j such that x lies inbetween xloc[j] and xloc[j+1].
1111 : // unit offset of increasing ordered array xloc assumed.
1112 : // n is the length of the array (xloc[n] highest element).
1113 :
1114 : int MSTWpdf::locate(double xloc[],int n,double x) {
1115 : int ju,jm,jl(0),j;
1116 0 : ju=n+1;
1117 :
1118 0 : while (ju-jl>1) {
1119 0 : jm=(ju+jl)/2; // compute a mid point.
1120 0 : if ( x>= xloc[jm])
1121 0 : jl=jm;
1122 : else ju=jm;
1123 : }
1124 0 : if (x==xloc[1]) j=1;
1125 0 : else if (x==xloc[n]) j=n-1;
1126 : else j=jl;
1127 :
1128 0 : return j;
1129 : }
1130 :
1131 : //--------------------------------------------------------------------------
1132 :
1133 : // Returns the estimate of the derivative at x1 obtained by a polynomial
1134 : // interpolation using the three points (x_i,y_i).
1135 :
1136 : double MSTWpdf::polderivative1(double x1, double x2, double x3, double y1,
1137 : double y2, double y3) {
1138 :
1139 0 : return (x3*x3*(y1-y2)+2.0*x1*(x3*(-y1+y2)+x2*(y1-y3))+x2*x2*(-y1+y3)
1140 0 : +x1*x1*(-y2+y3))/((x1-x2)*(x1-x3)*(x2-x3));
1141 :
1142 : }
1143 :
1144 : //--------------------------------------------------------------------------
1145 :
1146 : // Returns the estimate of the derivative at x2 obtained by a polynomial
1147 : // interpolation using the three points (x_i,y_i).
1148 :
1149 : double MSTWpdf::polderivative2(double x1, double x2, double x3, double y1,
1150 : double y2, double y3) {
1151 :
1152 0 : return (x3*x3*(y1-y2)-2.0*x2*(x3*(y1-y2)+x1*(y2-y3))+x2*x2*(y1-y3)
1153 0 : +x1*x1*(y2-y3))/((x1-x2)*(x1-x3)*(x2-x3));
1154 :
1155 : }
1156 :
1157 : //--------------------------------------------------------------------------
1158 :
1159 : // Returns the estimate of the derivative at x3 obtained by a polynomial
1160 : // interpolation using the three points (x_i,y_i).
1161 :
1162 : double MSTWpdf::polderivative3(double x1, double x2, double x3, double y1,
1163 : double y2, double y3) {
1164 :
1165 0 : return (x3*x3*(-y1+y2)+2.0*x2*x3*(y1-y3)+x1*x1*(y2-y3)+x2*x2*(-y1+y3)
1166 0 : +2.0*x1*x3*(-y2+y3))/((x1-x2)*(x1-x3)*(x2-x3));
1167 :
1168 : }
1169 :
1170 : //==========================================================================
1171 :
1172 : // The CTEQ6pdf class.
1173 : // Code for handling CTEQ6L, CTEQ6L1, CTEQ66.00, CT09MC1, CT09MC2, (CT09MCS?).
1174 :
1175 : // Constants: could be changed here if desired, but normally should not.
1176 : // These are of technical nature, as described for each.
1177 :
1178 : // Stay away from xMin, xMax, Qmin, Qmax limits.
1179 : const double CTEQ6pdf::EPSILON = 1e-6;
1180 :
1181 : // Assumed approximate power of small-x behaviour for interpolation.
1182 : const double CTEQ6pdf::XPOWER = 0.3;
1183 :
1184 : //--------------------------------------------------------------------------
1185 :
1186 : // Initialize PDF: read in data grid from file.
1187 :
1188 : void CTEQ6pdf::init(int iFitIn, string xmlPath, Info* infoPtr) {
1189 :
1190 : // Choice of fit among possibilities.
1191 0 : iFit = iFitIn;
1192 :
1193 : // Select which data file to read for current fit.
1194 0 : if (xmlPath[ xmlPath.length() - 1 ] != '/') xmlPath += "/";
1195 0 : string fileName = " ";
1196 0 : if (iFit == 1) fileName = "cteq6l.tbl";
1197 0 : if (iFit == 2) fileName = "cteq6l1.tbl";
1198 0 : if (iFit == 3) fileName = "ctq66.00.pds";
1199 0 : if (iFit == 4) fileName = "ct09mc1.pds";
1200 0 : if (iFit == 5) fileName = "ct09mc2.pds";
1201 0 : if (iFit == 6) fileName = "ct09mcs.pds";
1202 0 : bool isPdsGrid = (iFit > 2);
1203 :
1204 : // Open data file.
1205 0 : ifstream pdfgrid( (xmlPath + fileName).c_str() );
1206 0 : if (!pdfgrid.good()) {
1207 0 : if (infoPtr != 0) infoPtr->errorMsg("Error from CTEQ6pdf::init: "
1208 0 : "did not find parametrization file ", fileName);
1209 0 : else cout << " Error from CTEQ6pdf::init: "
1210 0 : << "did not find parametrization file " << fileName << endl;
1211 0 : isSet = false;
1212 0 : return;
1213 : }
1214 :
1215 : // Read in common information.
1216 0 : int iDum;
1217 0 : double orderTmp, nQTmp, qTmp, rDum;
1218 0 : string line;
1219 0 : getline( pdfgrid, line);
1220 0 : getline( pdfgrid, line);
1221 0 : getline( pdfgrid, line);
1222 0 : istringstream is1(line);
1223 0 : is1 >> orderTmp >> nQTmp >> lambda >> mQ[1] >> mQ[2] >> mQ[3]
1224 0 : >> mQ[4] >> mQ[5] >> mQ[6];
1225 0 : order = int(orderTmp + 0.5);
1226 0 : nQuark = int(nQTmp + 0.5);
1227 0 : getline( pdfgrid, line);
1228 :
1229 : // Read in information for the .pds grid format.
1230 0 : if (isPdsGrid) {
1231 0 : getline( pdfgrid, line);
1232 0 : istringstream is2(line);
1233 0 : is2 >> iDum >> iDum >> iDum >> nfMx >> mxVal >> iDum;
1234 0 : if (mxVal > 4) mxVal = 3;
1235 0 : getline( pdfgrid, line);
1236 0 : getline( pdfgrid, line);
1237 0 : istringstream is3(line);
1238 0 : is3 >> nX >> nT >> iDum >> nG >> iDum;
1239 0 : for (int i = 0; i < nG + 2; ++i) getline( pdfgrid, line);
1240 0 : getline( pdfgrid, line);
1241 0 : istringstream is4(line);
1242 0 : is4 >> qIni >> qMax;
1243 0 : for (int iT = 0; iT <= nT; ++iT) {
1244 0 : getline( pdfgrid, line);
1245 0 : istringstream is5(line);
1246 0 : is5 >> qTmp;
1247 0 : tv[iT] = log( log( qTmp/lambda));
1248 0 : }
1249 0 : getline( pdfgrid, line);
1250 0 : getline( pdfgrid, line);
1251 0 : istringstream is6(line);
1252 0 : is6 >> xMin >> rDum;
1253 : int nPackX = 6;
1254 0 : xv[0] = 0.;
1255 0 : for (int iXrng = 0; iXrng < int( (nX + nPackX - 1) / nPackX); ++iXrng) {
1256 0 : getline( pdfgrid, line);
1257 0 : istringstream is7(line);
1258 0 : for (int iX = nPackX * iXrng + 1; iX <= nPackX * (iXrng + 1); ++iX)
1259 0 : if (iX <= nX) is7 >> xv[iX];
1260 0 : }
1261 0 : }
1262 :
1263 : // Read in information for the .tbl grid format.
1264 : else {
1265 0 : mxVal = 2;
1266 0 : getline( pdfgrid, line);
1267 0 : istringstream is2(line);
1268 0 : is2 >> nX >> nT >> nfMx;
1269 0 : getline( pdfgrid, line);
1270 0 : getline( pdfgrid, line);
1271 0 : istringstream is3(line);
1272 0 : is3 >> qIni >> qMax;
1273 : int nPackT = 6;
1274 0 : for (int iTrng = 0; iTrng < int( (nT + nPackT) / nPackT); ++iTrng) {
1275 0 : getline( pdfgrid, line);
1276 0 : istringstream is4(line);
1277 0 : for (int iT = nPackT * iTrng; iT < nPackT * (iTrng + 1); ++iT)
1278 0 : if (iT <= nT) {
1279 0 : is4 >> qTmp;
1280 0 : tv[iT] = log( log( qTmp / lambda) );
1281 0 : }
1282 0 : }
1283 0 : getline( pdfgrid, line);
1284 0 : getline( pdfgrid, line);
1285 0 : istringstream is5(line);
1286 0 : is5 >> xMin;
1287 : int nPackX = 6;
1288 0 : for (int iXrng = 0; iXrng < int( (nX + nPackX) / nPackX); ++iXrng) {
1289 0 : getline( pdfgrid, line);
1290 0 : istringstream is6(line);
1291 0 : for (int iX = nPackX * iXrng; iX < nPackX * (iXrng + 1); ++iX)
1292 0 : if (iX <= nX) is6 >> xv[iX];
1293 0 : }
1294 0 : }
1295 :
1296 : // Read in the grid proper.
1297 0 : getline( pdfgrid, line);
1298 0 : int nBlk = (nX + 1) * (nT + 1);
1299 0 : int nPts = nBlk * (nfMx + 1 + mxVal);
1300 0 : int nPack = (isPdsGrid) ? 6 : 5;
1301 0 : for (int iRng = 0; iRng < int( (nPts + nPack - 1) / nPack); ++iRng) {
1302 0 : getline( pdfgrid, line);
1303 0 : istringstream is8(line);
1304 0 : for (int i = nPack * iRng + 1; i <= nPack * (iRng + 1); ++i)
1305 0 : if (i <= nPts) is8 >> upd[i];
1306 0 : }
1307 :
1308 : // Initialize x grid mapped to x^0.3.
1309 0 : xvpow[0] = 0.;
1310 0 : for (int iX = 1; iX <= nX; ++iX) xvpow[iX] = pow(xv[iX], XPOWER);
1311 :
1312 : // Set x and Q borders with some margin.
1313 0 : xMinEps = xMin * (1. + EPSILON);
1314 0 : xMaxEps = 1. - EPSILON;
1315 0 : qMinEps = qIni * (1. + EPSILON);
1316 0 : qMaxEps = qMax * (1. - EPSILON);
1317 :
1318 : // Initialize (x, Q) values of previous call.
1319 0 : xLast = 0.;
1320 0 : qLast = 0.;
1321 :
1322 0 : }
1323 :
1324 : //--------------------------------------------------------------------------
1325 :
1326 : // Update PDF values.
1327 :
1328 : void CTEQ6pdf::xfUpdate(int , double x, double Q2) {
1329 :
1330 : // Update using CTEQ6 routine, within allowed (x, q) range.
1331 0 : double xEps = max( xMinEps, x);
1332 0 : double qEps = max( qMinEps, min( qMaxEps, sqrtpos(Q2) ) );
1333 :
1334 : // Gluon:
1335 0 : double glu = xEps * parton6( 0, xEps, qEps);
1336 : // Sea quarks (note wrong order u, d):
1337 0 : double bot = xEps * parton6( 5, xEps, qEps);
1338 0 : double chm = xEps * parton6( 4, xEps, qEps);
1339 0 : double str = xEps * parton6( 3, xEps, qEps);
1340 0 : double usea = xEps * parton6(-1, xEps, qEps);
1341 0 : double dsea = xEps * parton6(-2, xEps, qEps);
1342 : // Valence quarks:
1343 0 : double upv = xEps * parton6( 1, xEps, qEps) - usea;
1344 0 : double dnv = xEps * parton6( 2, xEps, qEps) - dsea;
1345 :
1346 : // Transfer to Pythia notation.
1347 0 : xg = glu;
1348 0 : xu = upv + usea;
1349 0 : xd = dnv + dsea;
1350 0 : xubar = usea;
1351 0 : xdbar = dsea;
1352 0 : xs = str;
1353 0 : xsbar = str;
1354 0 : xc = chm;
1355 0 : xb = bot;
1356 0 : xgamma = 0.;
1357 :
1358 : // Subdivision of valence and sea.
1359 0 : xuVal = upv;
1360 0 : xuSea = usea;
1361 0 : xdVal = dnv;
1362 0 : xdSea = dsea;
1363 :
1364 : // idSav = 9 to indicate that all flavours reset.
1365 0 : idSav = 9;
1366 :
1367 0 : }
1368 :
1369 : //--------------------------------------------------------------------------
1370 :
1371 : // Returns the PDF value for parton of flavour iParton at x, q.
1372 :
1373 : double CTEQ6pdf::parton6(int iParton, double x, double q) {
1374 :
1375 : // Put zero for large x. Parton table and interpolation variables.
1376 0 : if (x > xMaxEps) return 0.;
1377 0 : int iP = (iParton > mxVal) ? -iParton : iParton;
1378 0 : double ss = pow( x, XPOWER);
1379 0 : double tt = log( log(q / lambda) );
1380 :
1381 : // Find location in grid.Skip if same as in latest call.
1382 0 : if (x != xLast || q != qLast) {
1383 :
1384 : // Binary search in x grid.
1385 0 : iGridX = 0;
1386 0 : iGridLX = -1;
1387 0 : int ju = nX + 1;
1388 : int jm = 0;
1389 0 : while (ju - iGridLX > 1 && jm >= 0) {
1390 0 : jm = (ju + iGridLX) / 2;
1391 0 : if (x >= xv[jm]) iGridLX = jm;
1392 : else ju = jm;
1393 : }
1394 :
1395 : // Separate acceptable from unacceptable grid points.
1396 0 : if (iGridLX <= -1) return 0.;
1397 0 : else if (iGridLX == 0) iGridX = 0;
1398 0 : else if (iGridLX <= nX - 2) iGridX = iGridLX - 1;
1399 0 : else if (iGridLX == nX - 1) iGridX = iGridLX - 2;
1400 0 : else return 0.;
1401 :
1402 : // Expressions for interpolation in x Grid.
1403 0 : if (iGridLX > 1 && iGridLX < nX - 1) {
1404 0 : double svec1 = xvpow[iGridX];
1405 0 : double svec2 = xvpow[iGridX+1];
1406 0 : double svec3 = xvpow[iGridX+2];
1407 0 : double svec4 = xvpow[iGridX+3];
1408 0 : double s12 = svec1 - svec2;
1409 0 : double s13 = svec1 - svec3;
1410 0 : xConst[8] = svec2 - svec3;
1411 0 : double s24 = svec2 - svec4;
1412 0 : double s34 = svec3 - svec4;
1413 0 : xConst[6] = ss - svec2;
1414 0 : xConst[7] = ss - svec3;
1415 0 : xConst[0] = s13 / xConst[8];
1416 0 : xConst[1] = s12 / xConst[8];
1417 0 : xConst[2] = s34 / xConst[8];
1418 0 : xConst[3] = s24 / xConst[8];
1419 0 : double s1213 = s12 + s13;
1420 0 : double s2434 = s24 + s34;
1421 0 : double sdet = s12 * s34 - s1213 * s2434;
1422 0 : double tmp = xConst[6] * xConst[7] / sdet;
1423 0 : xConst[4] = (s34 * xConst[6] - s2434 * xConst[7]) * tmp / s12;
1424 0 : xConst[5] = (s1213 * xConst[6] - s12 * xConst[7]) * tmp / s34;
1425 0 : }
1426 :
1427 : // Binary search in Q grid.
1428 0 : iGridQ = 0;
1429 0 : iGridLQ = -1;
1430 0 : ju = nT + 1;
1431 : jm = 0;
1432 0 : while (ju - iGridLQ > 1 && jm >= 0) {
1433 0 : jm = (ju + iGridLQ) / 2;
1434 0 : if (tt >= tv[jm]) iGridLQ = jm;
1435 : else ju = jm;
1436 : }
1437 0 : if (iGridLQ == 0) iGridQ = 0;
1438 0 : else if (iGridLQ <= nT - 2) iGridQ = iGridLQ - 1;
1439 0 : else iGridQ = nT - 3;
1440 :
1441 : // Expressions for interpolation in Q Grid.
1442 0 : if (iGridLQ > 0 && iGridLQ < nT - 1) {
1443 0 : double tvec1 = tv[iGridQ];
1444 0 : double tvec2 = tv[iGridQ+1];
1445 0 : double tvec3 = tv[iGridQ+2];
1446 0 : double tvec4 = tv[iGridQ+3];
1447 0 : double t12 = tvec1 - tvec2;
1448 0 : double t13 = tvec1 - tvec3;
1449 0 : tConst[8] = tvec2 - tvec3;
1450 0 : double t24 = tvec2 - tvec4;
1451 0 : double t34 = tvec3 - tvec4;
1452 0 : tConst[6] = tt - tvec2;
1453 0 : tConst[7] = tt - tvec3;
1454 0 : double tmp1 = t12 + t13;
1455 0 : double tmp2 = t24 + t34;
1456 0 : double tdet = t12 * t34 - tmp1 * tmp2;
1457 0 : tConst[0] = t13 / tConst[8];
1458 0 : tConst[1] = t12 / tConst[8];
1459 0 : tConst[2] = t34 / tConst[8];
1460 0 : tConst[3] = t24 / tConst[8];
1461 0 : tConst[4] = (t34 * tConst[6] - tmp2 * tConst[7]) / t12
1462 0 : * tConst[6] * tConst[7] / tdet;
1463 0 : tConst[5] = (tmp1 * tConst[6] - t12 * tConst[7]) / t34
1464 0 : * tConst[6] * tConst[7] / tdet;
1465 0 : }
1466 :
1467 : // Save x and q values so do not have to redo same again.
1468 0 : xLast = x;
1469 0 : qLast = q;
1470 0 : }
1471 :
1472 : // Jump to here if x and q are the same as for the last call.
1473 0 : int jtmp = ( (iP + nfMx) * (nT + 1) + (iGridQ - 1) ) * (nX + 1) + iGridX + 1;
1474 :
1475 : // Interpolate in x space for four different q values.
1476 0 : for(int it = 1; it <= 4; ++it) {
1477 0 : int j1 = jtmp + it * (nX + 1);
1478 0 : if (iGridX == 0) {
1479 0 : double fij[5];
1480 0 : fij[1] = 0.;
1481 0 : fij[2] = upd[j1+1] * pow2(xv[1]);
1482 0 : fij[3] = upd[j1+2] * pow2(xv[2]);
1483 0 : fij[4] = upd[j1+3] * pow2(xv[3]);
1484 0 : double fX = polint4F( &xvpow[0], &fij[1], ss);
1485 0 : fVec[it] = (x > 0.) ? fX / pow2(x) : 0.;
1486 0 : } else if (iGridLX==nX-1) {
1487 0 : fVec[it] = polint4F( &xvpow[nX-3], &upd[j1], ss);
1488 0 : } else {
1489 0 : double sf2 = upd[j1+1];
1490 0 : double sf3 = upd[j1+2];
1491 0 : double g1 = sf2 * xConst[0] - sf3 * xConst[1];
1492 0 : double g4 = -sf2 * xConst[2] + sf3 * xConst[3];
1493 0 : fVec[it] = (xConst[4] * (upd[j1] - g1) + xConst[5] * (upd[j1+3] - g4)
1494 0 : + sf2 * xConst[7] - sf3 * xConst[6]) / xConst[8];
1495 : }
1496 : }
1497 :
1498 : // Interpolate in q space for x-interpolated values found above.
1499 : double ff;
1500 0 : if( iGridLQ <= 0 ) {
1501 0 : ff = polint4F( &tv[0], &fVec[1], tt);
1502 0 : } else if (iGridLQ >= nT - 1) {
1503 0 : ff=polint4F( &tv[nT-3], &fVec[1], tt);
1504 0 : } else {
1505 0 : double tf2 = fVec[2];
1506 0 : double tf3 = fVec[3];
1507 0 : double g1 = tf2 * tConst[0] - tf3 * tConst[1];
1508 0 : double g4 = -tf2 * tConst[2] + tf3 * tConst[3];
1509 0 : ff = (tConst[4] * (fVec[1] - g1) + tConst[5] * (fVec[4] - g4)
1510 0 : + tf2 * tConst[7] - tf3 * tConst[6]) / tConst[8];
1511 : }
1512 :
1513 : // Done.
1514 : return ff;
1515 0 : }
1516 :
1517 : //--------------------------------------------------------------------------
1518 :
1519 : // The POLINT4 routine is based on the POLINT routine from "Numerical Recipes",
1520 : // but assuming N=4, and ignoring the error estimation.
1521 : // Suggested by Z. Sullivan.
1522 :
1523 : double CTEQ6pdf::polint4F(double xa[],double ya[],double x) {
1524 :
1525 : double y, h1, h2, h3, h4, w, den, d1, c1, d2, c2, d3, c3, cd1, cc1,
1526 : cd2, cc2, dd1, dc1;
1527 :
1528 0 : h1 = xa[0] - x;
1529 0 : h2 = xa[1] - x;
1530 0 : h3 = xa[2] - x;
1531 0 : h4 = xa[3] - x;
1532 :
1533 0 : w = ya[1] - ya[0];
1534 0 : den = w / (h1 - h2);
1535 0 : d1 = h2 * den;
1536 0 : c1 = h1 * den;
1537 :
1538 0 : w = ya[2] - ya[1];
1539 0 : den = w / (h2 - h3);
1540 0 : d2 = h3 * den;
1541 0 : c2 = h2 * den;
1542 :
1543 0 : w = ya[3] - ya[2];
1544 0 : den = w / (h3 - h4);
1545 0 : d3 = h4 * den;
1546 0 : c3 = h3 * den;
1547 :
1548 0 : w = c2 - d1;
1549 0 : den = w / (h1 - h3);
1550 0 : cd1 = h3 * den;
1551 0 : cc1 = h1 * den;
1552 :
1553 0 : w = c3 - d2;
1554 0 : den = w / (h2 - h4);
1555 0 : cd2 = h4 * den;
1556 0 : cc2 = h2 * den;
1557 :
1558 0 : w = cc2 - cd1;
1559 0 : den = w / (h1 - h4);
1560 0 : dd1 = h4 * den;
1561 0 : dc1 = h1 * den;
1562 :
1563 0 : if (h3 + h4 < 0.) y = ya[3] + d3 + cd2 + dd1;
1564 0 : else if (h2 + h3 < 0.) y = ya[2] + d2 + cd1 + dc1;
1565 0 : else if (h1 + h2 < 0.) y = ya[1] + c2 + cd1 + dc1;
1566 0 : else y = ya[0] + c1 + cc1 + dc1;
1567 :
1568 0 : return y;
1569 :
1570 : }
1571 :
1572 : //==========================================================================
1573 :
1574 : // SA Unresolved proton: equivalent photon spectrum from
1575 : // V.M. Budnev, I.F. Ginzburg, G.V. Meledin and V.G. Serbo,
1576 : // Phys. Rept. 15 (1974/1975) 181.
1577 :
1578 : // Constants:
1579 : const double ProtonPoint::ALPHAEM = 0.00729735;
1580 : const double ProtonPoint::Q2MAX = 2.0;
1581 : const double ProtonPoint::Q20 = 0.71;
1582 : const double ProtonPoint::A = 7.16;
1583 : const double ProtonPoint::B = -3.96;
1584 : const double ProtonPoint::C = 0.028;
1585 :
1586 : //--------------------------------------------------------------------------
1587 :
1588 : // Gives a generic Q2-independent equivalent photon spectrum.
1589 :
1590 : void ProtonPoint::xfUpdate(int , double x, double /*Q2*/ ) {
1591 :
1592 : // Photon spectrum
1593 0 : double tmpQ2Min = 0.88 * pow2(x);
1594 0 : double phiMax = phiFunc(x, Q2MAX / Q20);
1595 0 : double phiMin = phiFunc(x, tmpQ2Min / Q20);
1596 :
1597 : double fgm = 0;
1598 0 : if (phiMax < phiMin && m_infoPtr != 0) {
1599 0 : m_infoPtr->errorMsg("Error from ProtonPoint::xfUpdate: "
1600 : "phiMax - phiMin < 0!");
1601 0 : } else {
1602 : // Corresponds to: x*f(x)
1603 0 : fgm = (ALPHAEM / M_PI) * (1 - x) * (phiMax - phiMin);
1604 : }
1605 :
1606 : // Update values
1607 0 : xg = 0.;
1608 0 : xu = 0.;
1609 0 : xd = 0.;
1610 0 : xubar = 0.;
1611 0 : xdbar = 0.;
1612 0 : xs = 0.;
1613 0 : xsbar = 0.;
1614 0 : xc = 0.;
1615 0 : xb = 0.;
1616 0 : xgamma = fgm;
1617 :
1618 : // Subdivision of valence and sea.
1619 0 : xuVal = 0.;
1620 0 : xuSea = 0;
1621 0 : xdVal = 0.;
1622 0 : xdSea = 0;
1623 :
1624 : // idSav = 9 to indicate that all flavours reset.
1625 0 : idSav = 9;
1626 :
1627 0 : }
1628 :
1629 : //--------------------------------------------------------------------------
1630 :
1631 : // Function related to Q2 integration.
1632 :
1633 : double ProtonPoint::phiFunc(double x, double Q) {
1634 :
1635 0 : double tmpV = 1. + Q;
1636 : double tmpSum1 = 0;
1637 : double tmpSum2 = 0;
1638 0 : for (int k=1; k<4; ++k) {
1639 0 : tmpSum1 += 1. / (k * pow(tmpV, k));
1640 0 : tmpSum2 += pow(B, k) / (k * pow(tmpV, k));
1641 : }
1642 :
1643 0 : double tmpY = pow2(x) / (1 - x);
1644 0 : double funVal = (1 + A * tmpY) * (-1.*log(tmpV / Q) + tmpSum1)
1645 0 : + (1 - B) * tmpY / (4 * Q * pow(tmpV, 3))
1646 0 : + C * (1 + tmpY/4.)* (log((tmpV - B)/tmpV) + tmpSum2);
1647 :
1648 0 : return funVal;
1649 :
1650 : }
1651 :
1652 : //==========================================================================
1653 :
1654 : // Gives the GRV 1992 pi+ (leading order) parton distribution function set
1655 : // in parametrized form. Authors: Glueck, Reya and Vogt.
1656 : // Ref: M. Glueck, E. Reya and A. Vogt, Z. Phys. C53 (1992) 651.
1657 : // Allowed variable range: 0.25 GeV^2 < Q^2 < 10^8 GeV^2 and 10^-5 < x < 1.
1658 :
1659 : void GRVpiL::xfUpdate(int , double x, double Q2) {
1660 :
1661 : // Common expressions. Constrain Q2 for which parametrization is valid.
1662 : double mu2 = 0.25;
1663 : double lam2 = 0.232 * 0.232;
1664 0 : double s = (Q2 > mu2) ? log( log(Q2/lam2) / log(mu2/lam2) ) : 0.;
1665 0 : double s2 = s * s;
1666 0 : double x1 = 1. - x;
1667 0 : double xL = -log(x);
1668 0 : double xS = sqrt(x);
1669 :
1670 : // uv, dbarv.
1671 0 : double uv = (0.519 + 0.180 * s - 0.011 * s2) * pow(x, 0.499 - 0.027 * s)
1672 0 : * (1. + (0.381 - 0.419 * s) * xS) * pow(x1, 0.367 + 0.563 * s);
1673 :
1674 : // g.
1675 0 : double gl = ( pow(x, 0.482 + 0.341 * sqrt(s))
1676 0 : * ( (0.678 + 0.877 * s - 0.175 * s2) + (0.338 - 1.597 * s) * xS
1677 0 : + (-0.233 * s + 0.406 * s2) * x) + pow(s, 0.599)
1678 0 : * exp(-(0.618 + 2.070 * s) + sqrt(3.676 * pow(s, 1.263) * xL) ) )
1679 0 : * pow(x1, 0.390 + 1.053 * s);
1680 :
1681 : // sea: u, d, s.
1682 0 : double ub = pow(s, 0.55) * (1. - 0.748 * xS + (0.313 + 0.935 * s) * x)
1683 0 : * pow(x1, 3.359) * exp(-(4.433 + 1.301 * s) + sqrt((9.30 - 0.887 * s)
1684 0 : * pow(s, 0.56) * xL) ) / pow(xL, 2.538 - 0.763 * s);
1685 :
1686 : // c.
1687 0 : double chm = (s < 0.888) ? 0. : pow(s - 0.888, 1.02) * (1. + 1.008 * x)
1688 0 : * pow(x1, 1.208 + 0.771 * s) * exp(-(4.40 + 1.493 * s)
1689 0 : + sqrt( (2.032 + 1.901 * s) * pow(s, 0.39) * xL) );
1690 :
1691 : // b.
1692 0 : double bot = (s < 1.351) ? 0. : pow(s - 1.351, 1.03)
1693 0 : * pow(x1, 0.697 + 0.855 * s) * exp(-(4.51 + 1.490 * s)
1694 0 : + sqrt( (3.056 + 1.694 * s) * pow(s, 0.39) * xL) );
1695 :
1696 : // Update values.
1697 0 : xg = gl;
1698 0 : xu = uv + ub;
1699 0 : xd = ub;
1700 0 : xubar = ub;
1701 0 : xdbar = uv + ub;
1702 0 : xs = ub;
1703 0 : xsbar = ub;
1704 0 : xc = chm;
1705 0 : xb = bot;
1706 :
1707 : // Subdivision of valence and sea.
1708 0 : xuVal = uv;
1709 0 : xuSea = ub;
1710 0 : xdVal = uv;
1711 0 : xdSea = ub;
1712 :
1713 : // idSav = 9 to indicate that all flavours reset.
1714 0 : idSav = 9;
1715 :
1716 0 : }
1717 :
1718 : //==========================================================================
1719 :
1720 : // Pomeron PDF: simple Q2-independent parametrizations N x^a (1 - x)^b.
1721 :
1722 : //--------------------------------------------------------------------------
1723 :
1724 : // Calculate normalization factors once and for all.
1725 :
1726 : void PomFix::init() {
1727 :
1728 0 : normGluon = GammaReal(PomGluonA + PomGluonB + 2.)
1729 0 : / (GammaReal(PomGluonA + 1.) * GammaReal(PomGluonB + 1.));
1730 0 : normQuark = GammaReal(PomQuarkA + PomQuarkB + 2.)
1731 0 : / (GammaReal(PomQuarkA + 1.) * GammaReal(PomQuarkB + 1.));
1732 :
1733 0 : }
1734 :
1735 : //--------------------------------------------------------------------------
1736 :
1737 : // Gives a generic Q2-independent Pomeron PDF.
1738 :
1739 : void PomFix::xfUpdate(int , double x, double) {
1740 :
1741 : // Gluon and quark distributions.
1742 0 : double gl = normGluon * pow(x, PomGluonA) * pow( (1. - x), PomGluonB);
1743 0 : double qu = normQuark * pow(x, PomQuarkA) * pow( (1. - x), PomQuarkB);
1744 :
1745 : // Update values
1746 0 : xg = (1. - PomQuarkFrac) * gl;
1747 0 : xu = (PomQuarkFrac / (4. + 2. * PomStrangeSupp) ) * qu;
1748 0 : xd = xu;
1749 0 : xubar = xu;
1750 0 : xdbar = xu;
1751 0 : xs = PomStrangeSupp * xu;
1752 0 : xsbar = xs;
1753 0 : xc = 0.;
1754 0 : xb = 0.;
1755 :
1756 : // Subdivision of valence and sea.
1757 0 : xuVal = 0.;
1758 0 : xuSea = xu;
1759 0 : xdVal = 0.;
1760 0 : xdSea = xd;
1761 :
1762 : // idSav = 9 to indicate that all flavours reset.
1763 0 : idSav = 9;
1764 :
1765 0 : }
1766 :
1767 : //==========================================================================
1768 :
1769 : // Pomeron PDF: the H1 2006 Fit A and Fit B Q2-dependent parametrizations.
1770 :
1771 : //--------------------------------------------------------------------------
1772 :
1773 : void PomH1FitAB::init( int iFit, string xmlPath, Info* infoPtr) {
1774 :
1775 : // Open files from which grids should be read in.
1776 0 : if (xmlPath[ xmlPath.length() - 1 ] != '/') xmlPath += "/";
1777 0 : string dataFile = "pomH1FitBlo.data";
1778 0 : if (iFit == 1) dataFile = "pomH1FitA.data";
1779 0 : if (iFit == 2) dataFile = "pomH1FitB.data";
1780 0 : ifstream is( (xmlPath + dataFile).c_str() );
1781 0 : if (!is.good()) {
1782 0 : if (infoPtr != 0) infoPtr->errorMsg("Error from PomH1FitAB::init: "
1783 : "the H1 Pomeron parametrization file was not found");
1784 0 : else cout << " Error from PomH1FitAB::init: "
1785 0 : << "the H1 Pomeron parametrization file was not found" << endl;
1786 0 : isSet = false;
1787 0 : return;
1788 : }
1789 :
1790 : // Lower and upper bounds. Bin widths for logarithmic spacing.
1791 0 : nx = 100;
1792 0 : xlow = 0.001;
1793 0 : xupp = 0.99;
1794 0 : dx = log(xupp / xlow) / (nx - 1.);
1795 0 : nQ2 = 30;
1796 0 : Q2low = 1.0;
1797 0 : Q2upp = 30000.;
1798 0 : dQ2 = log(Q2upp / Q2low) / (nQ2 - 1.);
1799 :
1800 : // Read in quark data grid.
1801 0 : for (int i = 0; i < nx; ++i)
1802 0 : for (int j = 0; j < nQ2; ++j)
1803 0 : is >> quarkGrid[i][j];
1804 :
1805 : // Read in gluon data grid.
1806 0 : for (int i = 0; i < nx; ++i)
1807 0 : for (int j = 0; j < nQ2; ++j)
1808 0 : is >> gluonGrid[i][j];
1809 :
1810 : // Check for errors during read-in of file.
1811 0 : if (!is) {
1812 0 : if (infoPtr != 0) infoPtr->errorMsg("Error from PomH1FitAB::init: "
1813 : "the H1 Pomeron parametrization files could not be read");
1814 0 : else cout << " Error from PomH1FitAB::init: "
1815 0 : << "the H1 Pomeron parametrization files could not be read" << endl;
1816 0 : isSet = false;
1817 0 : return;
1818 : }
1819 :
1820 : // Done.
1821 0 : isSet = true;
1822 0 : return;
1823 0 : }
1824 :
1825 : //--------------------------------------------------------------------------
1826 :
1827 : void PomH1FitAB::xfUpdate(int , double x, double Q2) {
1828 :
1829 : // Retrict input to validity range.
1830 0 : double xt = min( xupp, max( xlow, x) );
1831 0 : double Q2t = min( Q2upp, max( Q2low, Q2) );
1832 :
1833 : // Lower grid point and distance above it.
1834 0 : double dlx = log( xt / xlow) / dx;
1835 0 : int i = min( nx - 2, int(dlx) );
1836 0 : dlx -= i;
1837 0 : double dlQ2 = log( Q2t / Q2low) / dQ2;
1838 0 : int j = min( nQ2 - 2, int(dlQ2) );
1839 0 : dlQ2 -= j;
1840 :
1841 : // Interpolate to derive quark PDF.
1842 0 : double qu = (1. - dlx) * (1. - dlQ2) * quarkGrid[i][j]
1843 0 : + dlx * (1. - dlQ2) * quarkGrid[i + 1][j]
1844 0 : + (1. - dlx) * dlQ2 * quarkGrid[i][j + 1]
1845 0 : + dlx * dlQ2 * quarkGrid[i + 1][j + 1];
1846 :
1847 : // Interpolate to derive gluon PDF.
1848 0 : double gl = (1. - dlx) * (1. - dlQ2) * gluonGrid[i][j]
1849 0 : + dlx * (1. - dlQ2) * gluonGrid[i + 1][j]
1850 0 : + (1. - dlx) * dlQ2 * gluonGrid[i][j + 1]
1851 0 : + dlx * dlQ2 * gluonGrid[i + 1][j + 1];
1852 :
1853 : // Update values.
1854 0 : xg = rescale * gl;
1855 0 : xu = rescale * qu;
1856 0 : xd = xu;
1857 0 : xubar = xu;
1858 0 : xdbar = xu;
1859 0 : xs = xu;
1860 0 : xsbar = xu;
1861 0 : xc = 0.;
1862 0 : xb = 0.;
1863 :
1864 : // Subdivision of valence and sea.
1865 0 : xuVal = 0.;
1866 0 : xuSea = xu;
1867 0 : xdVal = 0.;
1868 0 : xdSea = xu;
1869 :
1870 : // idSav = 9 to indicate that all flavours reset.
1871 0 : idSav = 9;
1872 :
1873 0 : }
1874 :
1875 : //==========================================================================
1876 :
1877 : // Pomeron PDF: the H1 2007 Jets Q2-dependent parametrization.
1878 :
1879 : //--------------------------------------------------------------------------
1880 :
1881 : void PomH1Jets::init( string xmlPath, Info* infoPtr) {
1882 :
1883 : // Open files from which grids should be read in.
1884 0 : if (xmlPath[ xmlPath.length() - 1 ] != '/') xmlPath += "/";
1885 0 : ifstream isg( (xmlPath + "pomH1JetsGluon.data").c_str() );
1886 0 : ifstream isq( (xmlPath + "pomH1JetsSinglet.data").c_str() );
1887 0 : ifstream isc( (xmlPath + "pomH1JetsCharm.data").c_str() );
1888 0 : if (!isg.good() || !isq.good() || !isc.good()) {
1889 0 : if (infoPtr != 0) infoPtr->errorMsg("Error from PomH1Jets::init: "
1890 : "the H1 Pomeron parametrization files were not found");
1891 0 : else cout << " Error from PomH1Jets::init: "
1892 0 : << "the H1 Pomeron parametrization files were not found" << endl;
1893 0 : isSet = false;
1894 0 : return;
1895 : }
1896 :
1897 : // Read in x and Q grids. Do interpolation logarithmically in Q2.
1898 0 : for (int i = 0; i < 100; ++i) {
1899 0 : isg >> setw(13) >> xGrid[i];
1900 : }
1901 0 : for (int j = 0; j < 88; ++j) {
1902 0 : isg >> setw(13) >> Q2Grid[j];
1903 0 : Q2Grid[j] = log( Q2Grid[j] );
1904 : }
1905 :
1906 : // Read in gluon data grid.
1907 0 : for (int j = 0; j < 88; ++j) {
1908 0 : for (int i = 0; i < 100; ++i) {
1909 0 : isg >> setw(13) >> gluonGrid[i][j];
1910 : }
1911 : }
1912 :
1913 : // Identical x and Q2 grid for singlet, so skip ahead.
1914 0 : double dummy;
1915 0 : for (int i = 0; i < 188; ++i) isq >> setw(13) >> dummy;
1916 :
1917 : // Read in singlet data grid.
1918 0 : for (int j = 0; j < 88; ++j) {
1919 0 : for (int i = 0; i < 100; ++i) {
1920 0 : isq >> setw(13) >> singletGrid[i][j];
1921 : }
1922 : }
1923 :
1924 : // Identical x and Q2 grid for charm, so skip ahead.
1925 0 : for (int i = 0; i < 188; ++i) isc >> setw(13) >> dummy;
1926 :
1927 : // Read in charm data grid.
1928 0 : for (int j = 0; j < 88; ++j) {
1929 0 : for (int i = 0; i < 100; ++i) {
1930 0 : isc >> setw(13) >> charmGrid[i][j];
1931 : }
1932 : }
1933 :
1934 : // Check for errors during read-in of files.
1935 0 : if (!isg || !isq || !isc) {
1936 0 : if (infoPtr != 0) infoPtr->errorMsg("Error from PomH1Jets::init: "
1937 : "the H1 Pomeron parametrization files could not be read");
1938 0 : else cout << " Error from PomH1Jets::init: "
1939 0 : << "the H1 Pomeron parametrization files could not be read" << endl;
1940 0 : isSet = false;
1941 0 : return;
1942 : }
1943 :
1944 : // Done.
1945 0 : isSet = true;
1946 0 : return;
1947 0 : }
1948 :
1949 : //--------------------------------------------------------------------------
1950 :
1951 : void PomH1Jets::xfUpdate(int , double x, double Q2) {
1952 :
1953 : // Find position in x array.
1954 0 : double xLog = log(x);
1955 : int i = 0;
1956 : double dx = 0.;
1957 0 : if (xLog <= xGrid[0]);
1958 0 : else if (xLog >= xGrid[99]) {
1959 : i = 98;
1960 : dx = 1.;
1961 0 : } else {
1962 0 : while (xLog > xGrid[i]) ++i;
1963 0 : --i;
1964 0 : dx = (xLog - xGrid[i]) / (xGrid[i + 1] - xGrid[i]);
1965 : }
1966 :
1967 : // Find position in y array.
1968 0 : double Q2Log = log(Q2);
1969 : int j = 0;
1970 : double dQ2 = 0.;
1971 0 : if (Q2Log <= Q2Grid[0]);
1972 0 : else if (Q2Log >= Q2Grid[87]) {
1973 : j = 86;
1974 : dQ2 = 1.;
1975 0 : } else {
1976 0 : while (Q2Log > Q2Grid[j]) ++j;
1977 0 : --j;
1978 0 : dQ2 = (Q2Log - Q2Grid[j]) / (Q2Grid[j + 1] - Q2Grid[j]);
1979 : }
1980 :
1981 : // Interpolate to derive gluon PDF.
1982 0 : double gl = (1. - dx) * (1. - dQ2) * gluonGrid[i][j]
1983 0 : + dx * (1. - dQ2) * gluonGrid[i + 1][j]
1984 0 : + (1. - dx) * dQ2 * gluonGrid[i][j + 1]
1985 0 : + dx * dQ2 * gluonGrid[i + 1][j + 1];
1986 :
1987 : // Interpolate to derive singlet PDF. (Sum of u, d, s, ubar, dbar, sbar.)
1988 0 : double sn = (1. - dx) * (1. - dQ2) * singletGrid[i][j]
1989 0 : + dx * (1. - dQ2) * singletGrid[i + 1][j]
1990 0 : + (1. - dx) * dQ2 * singletGrid[i][j + 1]
1991 0 : + dx * dQ2 * singletGrid[i + 1][j + 1];
1992 :
1993 : // Interpolate to derive charm PDF. (Charge-square times c and cbar.)
1994 0 : double ch = (1. - dx) * (1. - dQ2) * charmGrid[i][j]
1995 0 : + dx * (1. - dQ2) * charmGrid[i + 1][j]
1996 0 : + (1. - dx) * dQ2 * charmGrid[i][j + 1]
1997 0 : + dx * dQ2 * charmGrid[i + 1][j + 1];
1998 :
1999 : // Update values.
2000 0 : xg = rescale * gl;
2001 0 : xu = rescale * sn / 6.;
2002 0 : xd = xu;
2003 0 : xubar = xu;
2004 0 : xdbar = xu;
2005 0 : xs = xu;
2006 0 : xsbar = xu;
2007 0 : xc = rescale * ch * 9./8.;
2008 0 : xb = 0.;
2009 :
2010 : // Subdivision of valence and sea.
2011 0 : xuVal = 0.;
2012 0 : xuSea = xu;
2013 0 : xdVal = 0.;
2014 0 : xdSea = xd;
2015 :
2016 : // idSav = 9 to indicate that all flavours reset.
2017 0 : idSav = 9;
2018 :
2019 0 : }
2020 :
2021 : //==========================================================================
2022 :
2023 : // Gives electron (or muon, or tau) parton distribution.
2024 :
2025 : // Constants: alphaEM(0), m_e, m_mu, m_tau.
2026 : const double Lepton::ALPHAEM = 0.00729735;
2027 : const double Lepton::ME = 0.0005109989;
2028 : const double Lepton::MMU = 0.10566;
2029 : const double Lepton::MTAU = 1.77699;
2030 :
2031 : void Lepton::xfUpdate(int id, double x, double Q2) {
2032 :
2033 : // Squared mass of lepton species: electron, muon, tau.
2034 0 : if (!isInit) {
2035 0 : double mLep = ME;
2036 0 : if (abs(id) == 13) mLep = MMU;
2037 0 : if (abs(id) == 15) mLep = MTAU;
2038 0 : m2Lep = pow2( mLep );
2039 0 : isInit = true;
2040 0 : }
2041 :
2042 : // Electron inside electron, see R. Kleiss et al., in Z physics at
2043 : // LEP 1, CERN 89-08, p. 34
2044 0 : double xLog = log(max(1e-10,x));
2045 0 : double xMinusLog = log( max(1e-10, 1. - x) );
2046 0 : double Q2Log = log( max(3., Q2/m2Lep) );
2047 0 : double beta = (ALPHAEM / M_PI) * (Q2Log - 1.);
2048 0 : double delta = 1. + (ALPHAEM / M_PI) * (1.5 * Q2Log + 1.289868)
2049 0 : + pow2(ALPHAEM / M_PI) * (-2.164868 * Q2Log*Q2Log
2050 0 : + 9.840808 * Q2Log - 10.130464);
2051 0 : double fPrel = beta * pow(1. - x, beta - 1.) * sqrtpos( delta )
2052 0 : - 0.5 * beta * (1. + x) + 0.125 * beta*beta * ( (1. + x)
2053 0 : * (-4. * xMinusLog + 3. * xLog) - 4. * xLog / (1. - x) - 5. - x);
2054 :
2055 : // Zero distribution for very large x and rescale it for intermediate.
2056 0 : if (x > 1. - 1e-10) fPrel = 0.;
2057 0 : else if (x > 1. - 1e-7) fPrel *= pow(1000.,beta) / (pow(1000.,beta) - 1.);
2058 0 : xlepton = x * fPrel;
2059 :
2060 : // Photon inside electron (one possible scheme - primitive).
2061 0 : xgamma = (0.5 * ALPHAEM / M_PI) * Q2Log * (1. + pow2(1. - x));
2062 :
2063 : // idSav = 9 to indicate that all flavours reset.
2064 0 : idSav = 9;
2065 :
2066 0 : }
2067 :
2068 : //==========================================================================
2069 :
2070 : // The NNPDF class.
2071 : // Code for handling NNPDF2.3 QCD+QED LO
2072 : // Code provided by Juan Rojo and Stefano Carrazza.
2073 :
2074 : //--------------------------------------------------------------------------
2075 :
2076 : // Freeze PDFs below XMINGRID
2077 : const double NNPDF::fXMINGRID = 1e-9;
2078 :
2079 : //--------------------------------------------------------------------------
2080 :
2081 : // Initialize PDF: read in data grid from file.
2082 :
2083 : void NNPDF::init(int iFitIn, string xmlPath, Info* infoPtr) {
2084 :
2085 : // Choice of fit among possibilities.
2086 0 : iFit = iFitIn;
2087 :
2088 : // Select which data file to read for current fit.
2089 0 : if (xmlPath[ xmlPath.length() - 1 ] != '/') xmlPath += "/";
2090 0 : string fileName = " ";
2091 : // NNPDF2.3 LO QCD+QED, for two values of alphas
2092 0 : if (iFit == 1) fileName = "NNPDF23_lo_as_0130_qed_mem0.grid";
2093 0 : if (iFit == 2) fileName = "NNPDF23_lo_as_0119_qed_mem0.grid";
2094 : // NNPDF2.3 NLO QCD+QED
2095 0 : if (iFit == 3) fileName = "NNPDF23_nlo_as_0119_qed_mc_mem0.grid";
2096 : // NNPDF2.4 NLO QCD+QED
2097 0 : if (iFit == 4) fileName = "NNPDF23_nnlo_as_0119_qed_mc_mem0.grid";
2098 :
2099 : // Open data file.
2100 0 : fstream f;
2101 0 : f.open( (xmlPath + fileName).c_str(),ios::in);
2102 0 : if (f.fail()) {
2103 0 : if (infoPtr != 0) infoPtr->errorMsg("Error from NNPDF::init: "
2104 0 : "did not find data file ", fileName);
2105 0 : else cout << "Error: cannot open file " << (xmlPath + fileName) << endl;
2106 0 : isSet = false;
2107 0 : return;
2108 : }
2109 :
2110 : // Reading grid: removing header.
2111 0 : string tmp;
2112 0 : for (;;) {
2113 0 : getline(f,tmp);
2114 0 : if (tmp.find("NNPDF20intqed") != string::npos) {
2115 0 : getline(f,tmp);
2116 : break;
2117 : }
2118 : }
2119 :
2120 : // Get nx and x grid.
2121 0 : f >> fNX;
2122 0 : fXGrid = new double[fNX];
2123 0 : for (int ix = 0; ix < fNX; ix++) f >> fXGrid[ix];
2124 0 : fLogXGrid = new double[fNX];
2125 0 : for (int ix = 0; ix < fNX; ix++) fLogXGrid[ix] = log(fXGrid[ix]);
2126 :
2127 : // Get nQ2 and Q2 grid (ignorming first value).
2128 0 : f >> fNQ2;
2129 0 : f >> tmp;
2130 0 : fQ2Grid = new double[fNQ2];
2131 0 : for (int iq = 0; iq < fNQ2; iq++) f >> fQ2Grid[iq];
2132 0 : fLogQ2Grid = new double[fNQ2];
2133 0 : for (int iq = 0; iq < fNQ2; iq++) fLogQ2Grid[iq] = log(fQ2Grid[iq]);
2134 :
2135 : // Prepare grid array.
2136 0 : fPDFGrid = new double**[fNFL];
2137 0 : for (int i = 0; i < fNFL; i++) {
2138 0 : fPDFGrid[i] = new double*[fNX];
2139 0 : for (int j = 0; j < fNX; j++) {
2140 0 : fPDFGrid[i][j] = new double[fNQ2];
2141 0 : for (int z = 0; z < fNQ2; z++) fPDFGrid[i][j][z] = 0.0;
2142 : }
2143 : }
2144 :
2145 : // Check values of number of grid entries.
2146 0 : if (fNX<= 0 || fNX>100 || fNQ2<=0 || fNQ2>50) {
2147 0 : cout << "Error in NNPDF::init, Invalid grid values" << endl
2148 0 : << "fNX = " << fNX << endl << "fNQ2 = " << fNQ2 << endl
2149 0 : << "fNFL = " <<fNFL << endl;
2150 0 : isSet = false;
2151 0 : return;
2152 : }
2153 :
2154 : // Ignore replica number. Read PDF grid points.
2155 0 : f >> tmp;
2156 0 : for (int ix = 0; ix < fNX; ix++)
2157 0 : for (int iq = 0; iq < fNQ2; iq++)
2158 0 : for (int fl = 0; fl < fNFL; fl++)
2159 0 : f >> fPDFGrid[fl][ix][iq];
2160 0 : f.close();
2161 :
2162 : // Other vectors.
2163 0 : fRes = new double[fNFL];
2164 :
2165 0 : }
2166 :
2167 : //--------------------------------------------------------------------------
2168 :
2169 : void NNPDF::xfUpdate(int , double x, double Q2) {
2170 :
2171 : // Update using NNPDF routine, within allowed (x, q) range.
2172 0 : xfxevolve(x,Q2);
2173 :
2174 : // Then transfer to Pythia8 notation.
2175 0 : xg = fRes[6];
2176 0 : xu = fRes[8];
2177 0 : xd = fRes[7];
2178 0 : xubar = fRes[4];
2179 0 : xdbar = fRes[5];
2180 0 : xs = fRes[9];
2181 0 : xsbar = fRes[3];
2182 0 : xc = fRes[10];
2183 0 : xb = fRes[11];
2184 0 : xgamma = fRes[13];
2185 :
2186 : // Subdivision of valence and sea.
2187 0 : xuVal = xu - xubar;
2188 0 : xuSea = xubar;
2189 0 : xdVal = xd - xdbar;
2190 0 : xdSea = xdbar;
2191 :
2192 : // idSav = 9 to indicate that all flavours reset.
2193 0 : idSav = 9;
2194 :
2195 0 : }
2196 :
2197 : //--------------------------------------------------------------------------
2198 :
2199 : void NNPDF::xfxevolve(double x, double Q2) {
2200 :
2201 : // Freeze outside x-Q2 grid.
2202 0 : if (x < fXMINGRID || x > fXGrid[fNX-1]) {
2203 0 : if (x < fXMINGRID) x = fXMINGRID;
2204 0 : if (x > fXGrid[fNX-1]) x = fXGrid[fNX-1];
2205 : }
2206 0 : if (Q2 < fQ2Grid[0] || Q2 > fQ2Grid[fNQ2-1]) {
2207 0 : if (Q2 < fQ2Grid[0]) Q2 = fQ2Grid[0];
2208 0 : if (Q2 > fQ2Grid[fNQ2-1]) Q2 = fQ2Grid[fNQ2-1];
2209 : }
2210 :
2211 : // Find nearest points in the x-Q2 grid.
2212 : int minx = 0;
2213 0 : int maxx = fNX;
2214 0 : while (maxx-minx > 1) {
2215 0 : int midx = (minx+maxx)/2;
2216 0 : if (x < fXGrid[midx]) maxx = midx;
2217 : else minx = midx;
2218 : }
2219 : int ix = minx;
2220 : int minq = 0;
2221 0 : int maxq = fNQ2;
2222 0 : while (maxq-minq > 1) {
2223 0 : int midq = (minq+maxq)/2;
2224 0 : if (Q2 < fQ2Grid[midq]) maxq = midq;
2225 : else minq = midq;
2226 : }
2227 : int iq2 = minq;
2228 :
2229 : // Assign grid for interpolation. M,N -> order of polyN interpolation.
2230 0 : int ix1a[fM], ix2a[fN];
2231 0 : double x1a[fM], x2a[fN];
2232 0 : double ya[fM][fN];
2233 :
2234 0 : for (int i = 0; i < fM; i++) {
2235 0 : if (ix+1 >= fM/2 && ix+1 <= (fNX-fM/2)) ix1a[i] = ix+1 - fM/2 + i;
2236 0 : if (ix+1 < fM/2) ix1a[i] = i;
2237 0 : if (ix+1 > (fNX-fM/2)) ix1a[i] = (fNX-fM) + i;
2238 : // Check grids.
2239 0 : if (ix1a[i] < 0 || ix1a[i] >= fNX) {
2240 0 : cout << "Error in grids! i, ixia[i] = " << i << "\t" << ix1a[i] << endl;
2241 0 : return;
2242 : }
2243 : }
2244 :
2245 0 : for (int j = 0; j < fN; j++) {
2246 0 : if (iq2+1 >= fN/2 && iq2+1 <= (fNQ2-fN/2)) ix2a[j] = iq2+1 - fN/2 + j;
2247 0 : if (iq2+1 < fN/2) ix2a[j] = j;
2248 0 : if (iq2+1 > (fNQ2-fN/2)) ix2a[j] = (fNQ2-fN) + j;
2249 : // Check grids.
2250 0 : if (ix2a[j] < 0 || ix2a[j] >= fNQ2) {
2251 0 : cout << "Error in grids! j, ix2a[j] = " << j << "\t" << ix2a[j] << endl;
2252 0 : return;
2253 : }
2254 : }
2255 :
2256 : const double xch = 1e-1;
2257 : double x1;
2258 0 : if (x < xch) x1 = log(x);
2259 : else x1 = x;
2260 0 : double x2 = log(Q2);
2261 :
2262 0 : for (int ipdf = 0; ipdf < fNFL; ipdf++) {
2263 0 : fRes[ipdf] = 0.0;
2264 0 : for (int i = 0; i < fM; i++) {
2265 0 : if (x < xch) x1a[i] = fLogXGrid[ix1a[i]];
2266 0 : else x1a[i] = fXGrid[ix1a[i]];
2267 :
2268 0 : for (int j = 0; j < fN; j++) {
2269 0 : x2a[j] = fLogQ2Grid[ix2a[j]];
2270 0 : ya[i][j] = fPDFGrid[ipdf][ix1a[i]][ix2a[j]];
2271 : }
2272 : }
2273 :
2274 : // 2D polynomial interpolation.
2275 0 : double y = 0, dy = 0;
2276 0 : polin2(x1a,x2a,ya,x1,x2,y,dy);
2277 0 : fRes[ipdf] = y;
2278 0 : }
2279 :
2280 0 : }
2281 :
2282 : //--------------------------------------------------------------------------
2283 :
2284 : // 1D polynomial interpolation.
2285 :
2286 : void NNPDF::polint(double xa[], double yal[], int n, double x,
2287 : double& y, double& dy) {
2288 :
2289 : int ns = 0;
2290 0 : double dif = abs(x-xa[0]);
2291 0 : double c[fM > fN ? fM : fN];
2292 0 : double d[fM > fN ? fM : fN];
2293 :
2294 0 : for (int i = 0; i < n; i++) {
2295 0 : double dift = abs(x-xa[i]);
2296 0 : if (dift < dif) {
2297 : ns = i;
2298 : dif = dift;
2299 0 : }
2300 0 : c[i] = yal[i];
2301 0 : d[i] = yal[i];
2302 : }
2303 0 : y = yal[ns];
2304 0 : ns--;
2305 0 : for (int m = 1; m < n; m++) {
2306 0 : for (int i = 0; i < n-m; i++) {
2307 0 : double ho = xa[i]-x;
2308 0 : double hp = xa[i+m]-x;
2309 0 : double w = c[i+1]-d[i];
2310 0 : double den = ho-hp;
2311 0 : if (den == 0) {
2312 0 : cout << "NNPDF::polint, failure" << endl;
2313 0 : return;
2314 : }
2315 0 : den = w/den;
2316 0 : d[i] = hp*den;
2317 0 : c[i] = ho*den;
2318 0 : }
2319 0 : if (2*(ns+1) < n-m) dy = c[ns+1];
2320 : else {
2321 0 : dy = d[ns];
2322 0 : ns--;
2323 : }
2324 0 : y+=dy;
2325 : }
2326 0 : }
2327 :
2328 : //--------------------------------------------------------------------------
2329 :
2330 : // 2D polynomial interpolation.
2331 :
2332 : void NNPDF::polin2(double x1al[], double x2al[], double yal[][fN],
2333 : double x1, double x2, double& y, double& dy) {
2334 :
2335 0 : double yntmp[fN];
2336 0 : double ymtmp[fM];
2337 :
2338 0 : for (int j = 0; j < fM; j++) {
2339 0 : for (int k = 0; k < fN; k++) yntmp[k] = yal[j][k];
2340 0 : polint(x2al,yntmp,fN,x2,ymtmp[j],dy);
2341 : }
2342 0 : polint(x1al,ymtmp,fM,x1,y,dy);
2343 :
2344 0 : }
2345 :
2346 : //==========================================================================
2347 :
2348 : // LHAPDF plugin interface.
2349 :
2350 : //--------------------------------------------------------------------------
2351 :
2352 : // Constructor.
2353 :
2354 0 : LHAPDF::LHAPDF(int idIn, string pSet, Info* infoPtrIn) :
2355 0 : pdfPtr(0), infoPtr(infoPtrIn) {
2356 0 : isSet = false;
2357 0 : if (!infoPtr) return;
2358 :
2359 : // Determine the plugin library name.
2360 0 : if (pSet.size() < 8) {
2361 0 : infoPtr->errorMsg("Error from LHAPDF::LHAPDF: invalid pSet " + pSet);
2362 0 : return;
2363 : }
2364 0 : libName = pSet.substr(0, 7);
2365 0 : if (libName != "LHAPDF5" && libName != "LHAPDF6") {
2366 0 : infoPtr->errorMsg("Error from LHAPDF::LHAPDF: invalid pSet " + pSet);
2367 0 : return;
2368 : }
2369 0 : libName = "libpythia8lhapdf" + libName.substr(6) + ".so";
2370 :
2371 : // Determine the PDF set and member.
2372 0 : string set = pSet.substr(8);
2373 0 : int mem = 0;
2374 0 : size_t pos = set.find_last_of("/");
2375 0 : if (pos != string::npos) {
2376 0 : istringstream memStream(set.substr(pos + 1));
2377 0 : memStream >> mem;
2378 0 : }
2379 0 : set = set.substr(0, pos);
2380 :
2381 : // Load the PDF.
2382 0 : NewLHAPDF* newLHAPDF = (NewLHAPDF*)symbol("newLHAPDF");
2383 0 : if (!newLHAPDF) return;
2384 0 : pdfPtr = newLHAPDF(idIn, set, mem, infoPtr);
2385 0 : isSet = true;
2386 :
2387 0 : }
2388 :
2389 : //--------------------------------------------------------------------------
2390 :
2391 : // Destructor.
2392 :
2393 0 : LHAPDF::~LHAPDF() {
2394 0 : if (!infoPtr) return;
2395 0 : if (!isSet) return;
2396 :
2397 : // Delete the PDF.
2398 0 : DeleteLHAPDF* deleteLHAPDF = (DeleteLHAPDF*)symbol("deleteLHAPDF");
2399 0 : if (deleteLHAPDF) deleteLHAPDF(pdfPtr);
2400 :
2401 : // Close the plugin library if not needed by other instances.
2402 0 : map<string, pair<void*, int> >::iterator plugin =
2403 0 : infoPtr->plugins.find(libName);
2404 0 : if (plugin == infoPtr->plugins.end()) return;
2405 0 : --plugin->second.second;
2406 0 : if (plugin->second.first && plugin->second.second == 0) {
2407 0 : dlclose(plugin->second.first);
2408 0 : dlerror();
2409 0 : infoPtr->plugins.erase(plugin);
2410 0 : }
2411 :
2412 0 : }
2413 :
2414 : //--------------------------------------------------------------------------
2415 :
2416 : // Access a plugin library symbol.
2417 :
2418 : LHAPDF::Symbol LHAPDF::symbol(string symName) {
2419 0 : void *lib(0);
2420 : Symbol sym(0);
2421 : const char* error(0);
2422 0 : if (!infoPtr) return sym;
2423 :
2424 : // Load the library if not loaded.
2425 0 : map<string, pair<void*, int> >::iterator plugin =
2426 0 : infoPtr->plugins.find(libName);
2427 0 : if (plugin == infoPtr->plugins.end()) {
2428 0 : lib = dlopen(libName.c_str(), RTLD_LAZY);
2429 0 : error = dlerror();
2430 0 : }
2431 0 : if (error) {
2432 0 : infoPtr->errorMsg("Error from LHAPDF::symbol: " + string(error));
2433 0 : return sym;
2434 : }
2435 0 : if (plugin == infoPtr->plugins.end())
2436 0 : infoPtr->plugins[libName] = pair<void*, int>(lib, 1);
2437 : else {
2438 0 : lib = plugin->second.first;
2439 0 : ++plugin->second.second;
2440 : }
2441 0 : dlerror();
2442 :
2443 : // Load the symbol.
2444 0 : sym = (Symbol)dlsym(lib, symName.c_str());
2445 0 : error = dlerror();
2446 0 : if (error) infoPtr->errorMsg("Error from LHAPDF::symbol: " + string(error));
2447 0 : dlerror();
2448 0 : return sym;
2449 :
2450 0 : }
2451 :
2452 : //==========================================================================
2453 :
2454 : } // end namespace Pythia8
|
