source: HiSusy/trunk/Pythia8/pythia8170/src/PartonDistributions.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

first import of structure, PYTHIA8 and DELPHES

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