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