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