StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
PartonDistributions.cc
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 
13 namespace 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 
23 void 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 
54 double 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 
104 double 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 
140 double 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 
184 string LHAPDF::latestSetName = " ";
185 int LHAPDF::latestMember = -1;
186 int LHAPDF::latestNSet = 0;
187 
188 //--------------------------------------------------------------------------
189 
190 // Initialize a parton density function from LHAPDF.
191 
192 void 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 
230 void 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 
240 void 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 
284 void 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 
412 double 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 
423 double 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 
434 double 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 
458 void 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.
620 const int MSTWpdf::np = 12;
621 const int MSTWpdf::nx = 64;
622 const int MSTWpdf::nq = 48;
623 const int MSTWpdf::nqc0 = 4;
624 const int MSTWpdf::nqb0 = 14;
625 
626 // Range of (x, Q2) grid.
627 const double MSTWpdf::xmin = 1e-6;
628 const double MSTWpdf::xmax = 1.0;
629 const double MSTWpdf::qsqmin = 1.0;
630 const double MSTWpdf::qsqmax = 1e9;
631 
632 // Array of x values.
633 const 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.
642 const 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 
652 void 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 
948 void 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 
1002 double 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 
1096 double 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 
1136 double 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 
1191 int 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 
1213 double 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 
1226 double 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 
1239 double 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.
1256 const double CTEQ6pdf::EPSILON = 1e-6;
1257 
1258 // Assumed approximate power of small-x behaviour for interpolation.
1259 const double CTEQ6pdf::XPOWER = 0.3;
1260 
1261 //--------------------------------------------------------------------------
1262 
1263 // Initialize PDF: read in data grid from file.
1264 
1265 void 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 
1405 void 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 
1450 double 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 
1600 double 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:
1656 const double ProtonPoint::ALPHAEM = 0.00729735;
1657 const double ProtonPoint::Q2MAX = 2.0;
1658 const double ProtonPoint::Q20 = 0.71;
1659 const double ProtonPoint::A = 7.16;
1660 const double ProtonPoint::B = -3.96;
1661 const double ProtonPoint::C = 0.028;
1662 
1663 //--------------------------------------------------------------------------
1664 
1665 // Gives a generic Q2-independent equivalent photon spectrum.
1666 
1667 void 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 
1710 double 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 
1736 void 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 
1803 void 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 
1816 void 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 
1850 void 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 
1904 void 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 
1958 void 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 
2028 void 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.
2103 const double Lepton::ALPHAEM = 0.00729735;
2104 const double Lepton::ME = 0.0005109989;
2105 const double Lepton::MMU = 0.10566;
2106 const double Lepton::MTAU = 1.77699;
2107 
2108 void 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
ae
void ae(int tracks=-1, int hits=-1)
This function is to search for the next non-empty event and draw it by looping over StBFChain (readin...
Definition: Ed.C:102
Generated by   doxygen 1.8.5