Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / src / PartonDistributions.cxx
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