PYTHIA8/pythia8130/src/PartonDistributions.cxx

   1 // PartonDistributions.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2008 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,
   7 // GRV94L, CTEQ5L, LHAPDF  and Lepton classes.
   8
   9 #include "PartonDistributions.h"
  10 #include "LHAPDFInterface.h"
  11
  12 namespace Pythia8 {
  13
  14 //**************************************************************************
  15
  16 // Base class for parton distribution functions.
  17
  18 //*********
  19
  20 // Standard parton densities.
  21
  22 double PDF::xf(int id, double x, double Q2) {
  23
  24   // Need to update if flavour, x or Q2 changed.
  25   // Use idSav = 9 to indicate that ALL flavours are up-to-date.
  26   // Assume that flavour and antiflavour always updated simultaneously.
  27   if ( (abs(idSav) != abs(id) && idSav != 9) || x != xSav || Q2 != Q2Sav)
  28     {idSav = id; xfUpdate(id, x, Q2); xSav = x; Q2Sav = Q2;}
  29
  30   // Baryon beam.
  31   if (abs(idBeam) > 100) {
  32     int idNow = (idBeam > 0) ? id : -id;
  33     int idAbs = abs(id);
  34     if (idNow == 0 || idAbs == 21) return max(0., xg);
  35     if (idNow == 1) return max(0., xd);
  36     if (idNow == -1) return max(0., xdbar);
  37     if (idNow == 2) return max(0., xu);
  38     if (idNow == -2) return max(0., xubar);
  39     if (idAbs == 3) return max(0., xs);
  40     if (idAbs == 4) return max(0., xc);
  41     if (idAbs == 5) return max(0., xb);
  42     return 0.;
  43
  44   // Lepton beam.
  45   } else {
  46     if (id == idBeam ) return max(0., xlepton);
  47     if (abs(id) == 22) return max(0., xgamma);
  48     return 0.;
  49   }
  50
  51 }
  52
  53 //*********
  54
  55 // Only valence part of parton densities.
  56
  57 double PDF::xfVal(int id, double x, double Q2) {
  58
  59   // Need to update if flavour, x or Q2 changed.
  60   // Use idSav = 9 to indicate that ALL flavours are up-to-date.
  61   // Assume that flavour and antiflavour always updated simultaneously.
  62   if ( (abs(idSav) != abs(id) && idSav != 9) || x != xSav || Q2 != Q2Sav)
  63     {idSav = id; xfUpdate(id, x, Q2); xSav = x; Q2Sav = Q2;}
  64
  65   // Baryon beam.
  66   if (abs(idBeam) > 100) {
  67     int idNow = (idBeam > 0) ? id : -id;
  68     if (idNow == 1) return max(0., xdVal);
  69     if (idNow == 2) return max(0., xuVal);
  70     return 0.;
  71
  72   // Lepton beam.
  73   } else {
  74     if (id == idBeam) return max(0., xlepton);
  75     return 0.;
  76   }
  77
  78 }
  79
  80 //*********
  81
  82 // Only sea part of parton densities.
  83
  84 double PDF::xfSea(int id, double x, double Q2) {
  85
  86   // Need to update if flavour, x or Q2 changed.
  87   // Use idSav = 9 to indicate that ALL flavours are up-to-date.
  88   // Assume that flavour and antiflavour always updated simultaneously.
  89   if ( (abs(idSav) != abs(id) && idSav != 9) || x != xSav || Q2 != Q2Sav)
  90     {idSav = id; xfUpdate(id, x, Q2); xSav = x; Q2Sav = Q2;}
  91
  92   // Baryon beam.
  93   if (abs(idBeam) > 100) {
  94     int idNow = (idBeam > 0) ? id : -id;
  95     int idAbs = abs(id);
  96     if (idNow == 0 || idAbs == 21) return max(0., xg);
  97     if (idNow == 1) return max(0., xdSea);
  98     if (idNow == -1) return max(0., xdbar);
  99     if (idNow == 2) return max(0., xuSea);
 100     if (idNow == -2) return max(0., xubar);
 101     if (idAbs == 3) return max(0., xs);
 102     if (idAbs == 4) return max(0., xc);
 103     if (idAbs == 5) return max(0., xb);
 104     if (idAbs == 5) return max(0., xb);
 105     return 0.;
 106
 107   // Lepton beam.
 108   } else {
 109     if (abs(id) == 22) return max(0., xgamma);
 110     return 0.;
 111   }
 112
 113 }
 114
 115 //**************************************************************************
 116
 117 // Gives the GRV 94 L (leading order) parton distribution function set
 118 // in parametrized form. Authors: M. Glueck, E. Reya and A. Vogt.
 119 // Ref: M. Glueck, E. Reya and A. Vogt, Z.Phys. C67 (1995) 433.
 120
 121 void GRV94L::xfUpdate(int id, double x, double Q2) {
 122
 123   // Common expressions.
 124   double mu2  = 0.23;
 125   double lam2 = 0.2322 * 0.2322;
 126   double s    = log (log(Q2/lam2) / log(mu2/lam2));
 127   double ds   = sqrt(s);
 128   double s2   = s * s;
 129   double s3   = s2 * s;
 130
 131   // uv :
 132   double nu  =  2.284 + 0.802 * s + 0.055 * s2;
 133   double aku =  0.590 - 0.024 * s;
 134   double bku =  0.131 + 0.063 * s;
 135   double au  = -0.449 - 0.138 * s - 0.076 * s2;
 136   double bu  =  0.213 + 2.669 * s - 0.728 * s2;
 137   double cu  =  8.854 - 9.135 * s + 1.979 * s2;
 138   double du  =  2.997 + 0.753 * s - 0.076 * s2;
 139   double uv  = grvv (x, nu, aku, bku, au, bu, cu, du);
 140
 141   // dv :
 142   double nd  =  0.371 + 0.083 * s + 0.039 * s2;
 143   double akd =  0.376;
 144   double bkd =  0.486 + 0.062 * s;
 145   double ad  = -0.509 + 3.310 * s - 1.248 * s2;
 146   double bd  =  12.41 - 10.52 * s + 2.267 * s2;
 147   double cd  =  6.373 - 6.208 * s + 1.418 * s2;
 148   double dd  =  3.691 + 0.799 * s - 0.071 * s2;
 149   double dv  = grvv (x, nd, akd, bkd, ad, bd, cd, dd);
 150
 151   // udb :
 152   double alx =  1.451;
 153   double bex =  0.271;
 154   double akx =  0.410 - 0.232 * s;
 155   double bkx =  0.534 - 0.457 * s;
 156   double agx =  0.890 - 0.140 * s;
 157   double bgx = -0.981;
 158   double cx  =  0.320 + 0.683 * s;
 159   double dx  =  4.752 + 1.164 * s + 0.286 * s2;
 160   double ex  =  4.119 + 1.713 * s;
 161   double esx =  0.682 + 2.978 * s;
 162   double udb = grvw (x, s, alx, bex, akx, bkx, agx, bgx, cx,
 163     dx, ex, esx);
 164
 165   // del :
 166   double ne  =  0.082 + 0.014 * s + 0.008 * s2;
 167   double ake =  0.409 - 0.005 * s;
 168   double bke =  0.799 + 0.071 * s;
 169   double ae  = -38.07 + 36.13 * s - 0.656 * s2;
 170   double be  =  90.31 - 74.15 * s + 7.645 * s2;
 171   double ce  =  0.;
 172   double de  =  7.486 + 1.217 * s - 0.159 * s2;
 173   double del = grvv (x, ne, ake, bke, ae, be, ce, de);
 174
 175   // sb :
 176   double sts =  0.;
 177   double als =  0.914;
 178   double bes =  0.577;
 179   double aks =  1.798 - 0.596 * s;
 180   double as  = -5.548 + 3.669 * ds - 0.616 * s;
 181   double bs  =  18.92 - 16.73 * ds + 5.168 * s;
 182   double dst =  6.379 - 0.350 * s  + 0.142 * s2;
 183   double est =  3.981 + 1.638 * s;
 184   double ess =  6.402;
 185   double sb  = grvs (x, s, sts, als, bes, aks, as, bs, dst, est, ess);
 186
 187   // cb :
 188   double stc =  0.888;
 189   double alc =  1.01;
 190   double bec =  0.37;
 191   double akc =  0.;
 192   double ac  =  0.;
 193   double bc  =  4.24  - 0.804 * s;
 194   double dct =  3.46  - 1.076 * s;
 195   double ect =  4.61  + 1.49  * s;
 196   double esc =  2.555 + 1.961 * s;
 197   double chm = grvs (x, s, stc, alc, bec, akc, ac, bc, dct, ect, esc);
 198
 199   // bb :
 200   double stb =  1.351;
 201   double alb =  1.00;
 202   double beb =  0.51;
 203   double akb =  0.;
 204   double ab  =  0.;
 205   double bb  =  1.848;
 206   double dbt =  2.929 + 1.396 * s;
 207   double ebt =  4.71  + 1.514 * s;
 208   double esb =  4.02  + 1.239 * s;
 209   double bot = grvs (x, s, stb, alb, beb, akb, ab, bb, dbt, ebt, esb);
 210
 211   // gl :
 212   double alg =  0.524;
 213   double beg =  1.088;
 214   double akg =  1.742 - 0.930 * s;
 215   double bkg =                     - 0.399 * s2;
 216   double ag  =  7.486 - 2.185 * s;
 217   double bg  =  16.69 - 22.74 * s  + 5.779 * s2;
 218   double cg  = -25.59 + 29.71 * s  - 7.296 * s2;
 219   double dg  =  2.792 + 2.215 * s  + 0.422 * s2 - 0.104 * s3;
 220   double eg  =  0.807 + 2.005 * s;
 221   double esg =  3.841 + 0.316 * s;
 222   double gl  = grvw (x, s, alg, beg, akg, bkg, ag, bg, cg,
 223     dg, eg, esg);
 224
 225   // Update values
 226   xg    = gl;
 227   xu    = uv + 0.5*(udb - del);
 228   xd    = dv + 0.5*(udb + del);
 229   xubar = 0.5*(udb - del);
 230   xdbar = 0.5*(udb + del);
 231   xs    = sb;
 232   xc    = chm;
 233   xb    = bot;
 234
 235   // Subdivision of valence and sea.
 236   xuVal = uv;
 237   xuSea = xubar;
 238   xdVal = dv;
 239   xdSea = xdbar;
 240
 241   // idSav = 9 to indicate that all flavours reset. id change dummy.
 242   idSav = 9;
 243   id   = 0;
 244
 245 }
 246
 247 //*********
 248
 249 double GRV94L::grvv (double x, double n, double ak, double bk, double a,
 250    double b, double c, double d) {
 251
 252   double dx = sqrt(x);
 253   return n * pow(x, ak) * (1. + a * pow(x, bk) + x * (b + c * dx)) *
 254     pow(1. - x, d);
 255
 256 }
 257
 258 //*********
 259
 260 double GRV94L::grvw (double x, double s, double al, double be, double ak,
 261   double bk, double a, double b, double c, double d, double e, double es) {
 262
 263   double lx = log(1./x);
 264   return (pow(x, ak) * (a + x * (b + x * c)) * pow(lx, bk) + pow(s, al)
 265     * exp(-e + sqrt(es * pow(s, be) * lx))) * pow(1. - x, d);
 266
 267 }
 268
 269 //*********
 270
 271 double GRV94L::grvs (double x, double s, double sth, double al, double be,
 272   double ak, double ag, double b, double d, double e, double es) {
 273
 274   if(s <= sth) {
 275     return 0.;
 276   } else {
 277     double dx = sqrt(x);
 278     double lx = log(1./x);
 279     return pow(s - sth, al) / pow(lx, ak) * (1. + ag * dx + b * x) *
 280       pow(1. - x, d) * exp(-e + sqrt(es * pow(s, be) * lx));
 281   }
 282
 283 }
 284
 285 //**************************************************************************
 286
 287 // Gives the CTEQ 5 L (leading order) parton distribution function set
 288 // in parametrized form. Parametrization by J. Pumplin.
 289 // Ref: CTEQ Collaboration, H.L. Lai et al., Eur.Phys.J. C12 (2000) 375.
 290
 291 // The range of (x, Q) covered by this parametrization of the QCD
 292 // evolved parton distributions is 1E-6 < x < 1, 1.1 GeV < Q < 10 TeV.
 293 // In the current implementation, densities are frozen at borders.
 294
 295 void CTEQ5L::xfUpdate(int id, double x, double Q2) {
 296
 297   // Constrain x and Q2 to range for which parametrization is valid.
 298   double Q = sqrt( max( 1., min( 1e8, Q2) ) );
 299   x = max( 1e-6, min( 1.-1e-10, x) );
 300
 301   // Derived kinematical quantities.
 302   double y = - log(x);
 303   double u = log( x / 0.00001);
 304   double x1 = 1. - x;
 305   double x1L = log(1. - x);
 306   double sumUbarDbar = 0.;
 307
 308   // Parameters of parametrizations.
 309   const double Qmin[8] = { 0., 0., 0., 0., 0., 0., 1.3, 4.5};
 310   const double alpha[8] = { 0.2987216, 0.3407552, 0.4491863, 0.2457668,
 311     0.5293999, 0.3713141, 0.03712017, 0.004952010 };
 312   const double ut1[8] = { 4.971265, 2.612618, -0.4656819, 3.862583,
 313     0.1895615, 3.753257, 4.400772, 5.562568 };
 314   const double ut2[8] = { -1.105128, -1.258304e5, -274.2390, -1.265969,
 315     -3.069097, -1.113085, -1.356116, -1.801317 };
 316   const double am[8][9][3] = {
 317     // d.
 318     { {  0.5292616E+01, -0.2751910E+01, -0.2488990E+01 },
 319       {  0.9714424E+00,  0.1011827E-01, -0.1023660E-01 },
 320       { -0.1651006E+02,  0.7959721E+01,  0.8810563E+01 },
 321       { -0.1643394E+02,  0.5892854E+01,  0.9348874E+01 },
 322       {  0.3067422E+02,  0.4235796E+01, -0.5112136E+00 },
 323       {  0.2352526E+02, -0.5305168E+01, -0.1169174E+02 },
 324       { -0.1095451E+02,  0.3006577E+01,  0.5638136E+01 },
 325       { -0.1172251E+02, -0.2183624E+01,  0.4955794E+01 },
 326       {  0.1662533E-01,  0.7622870E-02, -0.4895887E-03 } },
 327     // u.
 328     { {  0.9905300E+00, -0.4502235E+00,  0.1624441E+00 },
 329       {  0.8867534E+00,  0.1630829E-01, -0.4049085E-01 },
 330       {  0.8547974E+00,  0.3336301E+00,  0.1371388E+00 },
 331       {  0.2941113E+00, -0.1527905E+01,  0.2331879E+00 },
 332       {  0.3384235E+02,  0.3715315E+01,  0.8276930E+00 },
 333       {  0.6230115E+01,  0.3134639E+01, -0.1729099E+01 },
 334       { -0.1186928E+01, -0.3282460E+00,  0.1052020E+00 },
 335       { -0.8545702E+01, -0.6247947E+01,  0.3692561E+01 },
 336       {  0.1724598E-01,  0.7120465E-02,  0.4003646E-04 } },
 337     // g.
 338     { {  0.1193572E+03, -0.3886845E+01, -0.1133965E+01 },
 339       { -0.9421449E+02,  0.3995885E+01,  0.1607363E+01 },
 340       {  0.4206383E+01,  0.2485954E+00,  0.2497468E+00 },
 341       {  0.1210557E+03, -0.3015765E+01, -0.1423651E+01 },
 342       { -0.1013897E+03, -0.7113478E+00,  0.2621865E+00 },
 343       { -0.1312404E+01, -0.9297691E+00, -0.1562531E+00 },
 344       {  0.1627137E+01,  0.4954111E+00, -0.6387009E+00 },
 345       {  0.1537698E+00, -0.2487878E+00,  0.8305947E+00 },
 346       {  0.2496448E-01,  0.2457823E-02,  0.8234276E-03 } },
 347     // ubar + dbar.
 348     { {  0.2647441E+02,  0.1059277E+02, -0.9176654E+00 },
 349       {  0.1990636E+01,  0.8558918E-01,  0.4248667E-01 },
 350       { -0.1476095E+02, -0.3276255E+02,  0.1558110E+01 },
 351       { -0.2966889E+01, -0.3649037E+02,  0.1195914E+01 },
 352       { -0.1000519E+03, -0.2464635E+01,  0.1964849E+00 },
 353       {  0.3718331E+02,  0.4700389E+02, -0.2772142E+01 },
 354       { -0.1872722E+02, -0.2291189E+02,  0.1089052E+01 },
 355       { -0.1628146E+02, -0.1823993E+02,  0.2537369E+01 },
 356       { -0.1156300E+01, -0.1280495E+00,  0.5153245E-01 } },
 357     // dbar/ubar.
 358     { { -0.6556775E+00,  0.2490190E+00,  0.3966485E-01 },
 359       {  0.1305102E+01, -0.1188925E+00, -0.4600870E-02 },
 360       { -0.2371436E+01,  0.3566814E+00, -0.2834683E+00 },
 361       { -0.6152826E+01,  0.8339877E+00, -0.7233230E+00 },
 362       { -0.8346558E+01,  0.2892168E+01,  0.2137099E+00 },
 363       {  0.1279530E+02,  0.1021114E+00,  0.5787439E+00 },
 364       {  0.5858816E+00, -0.1940375E+01, -0.4029269E+00 },
 365       { -0.2795725E+02, -0.5263392E+00,  0.1290229E+01 },
 366       {  0.0000000E+00,  0.0000000E+00,  0.0000000E+00 } },
 367     // sbar.
 368     { {  0.1580931E+01, -0.2273826E+01, -0.1822245E+01 },
 369       {  0.2702644E+01,  0.6763243E+00,  0.7231586E-02 },
 370       { -0.1857924E+02,  0.3907500E+01,  0.5850109E+01 },
 371       { -0.3044793E+02,  0.2639332E+01,  0.5566644E+01 },
 372       { -0.4258011E+01, -0.5429244E+01,  0.4418946E+00 },
 373       {  0.3465259E+02, -0.5532604E+01, -0.4904153E+01 },
 374       { -0.1658858E+02,  0.2923275E+01,  0.2266286E+01 },
 375       { -0.1149263E+02,  0.2877475E+01, -0.7999105E+00 },
 376       {  0.0000000E+00,  0.0000000E+00,  0.0000000E+00 } },
 377     // cbar.
 378     { { -0.8293661E+00, -0.3982375E+01, -0.6494283E-01 },
 379       {  0.2754618E+01,  0.8338636E+00, -0.6885160E-01 },
 380       { -0.1657987E+02,  0.1439143E+02, -0.6887240E+00 },
 381       { -0.2800703E+02,  0.1535966E+02, -0.7377693E+00 },
 382       { -0.6460216E+01, -0.4783019E+01,  0.4913297E+00 },
 383       {  0.3141830E+02, -0.3178031E+02,  0.7136013E+01 },
 384       { -0.1802509E+02,  0.1862163E+02, -0.4632843E+01 },
 385       { -0.1240412E+02,  0.2565386E+02, -0.1066570E+02 },
 386       {  0.0000000E+00,  0.0000000E+00,  0.0000000E+00 } },
 387     // bbar.
 388     { { -0.6031237E+01,  0.1992727E+01, -0.1076331E+01 },
 389       {  0.2933912E+01,  0.5839674E+00,  0.7509435E-01 },
 390       { -0.8284919E+01,  0.1488593E+01, -0.8251678E+00 },
 391       { -0.1925986E+02,  0.2805753E+01, -0.3015446E+01 },
 392       { -0.9480483E+01, -0.9767837E+00, -0.1165544E+01 },
 393       {  0.2193195E+02, -0.1788518E+02,  0.9460908E+01 },
 394       { -0.1327377E+02,  0.1201754E+02, -0.6277844E+01 },
 395       {  0.0000000E+00,  0.0000000E+00,  0.0000000E+00 },
 396       {  0.0000000E+00,  0.0000000E+00,  0.0000000E+00 } } };
 397
 398   // Loop over 8 different parametrizations. Check if inside allowed region.
 399   for (int i = 0; i < 8; ++i) {
 400     double answer = 0.;
 401     if (Q > max(Qmin[i], alpha[i])) {
 402
 403       // Evaluate answer.
 404       double tmp = log(Q / alpha[i]);
 405       double sb = log(tmp);
 406       double sb1 = sb - 1.2;
 407       double sb2 = sb1*sb1;
 408       double af[9];
 409       for (int j = 0; j < 9; ++j)
 410         af[j] = am[i][j][0] + sb1 * am[i][j][1] + sb2 * am[i][j][2];
 411       double part1 = af[1] * pow( y, 1. + 0.01 * af[4]) * (1. + af[8] * u);
 412       double part2 = af[0] * x1 + af[3] * x;
 413       double part3 = x * x1 * (af[5] + af[6] * x1 + af[7] * x * x1);
 414       double part4 = (ut2[i] < -100.) ? ut1[i] * x1L + af[2] * x1L
 415                    : ut1[i] * x1L + af[2] * log(x1 + exp(ut2[i]));
 416       answer       = x * exp( part1 + part2 + part3 + part4);
 417       answer      *= 1. - Qmin[i] / Q;
 418     }
 419
 420     // Store results.
 421     if (i == 0) xd = x * answer;
 422     else if (i == 1) xu = x * answer;
 423     else if (i == 2) xg = x * answer;
 424     else if (i == 3) sumUbarDbar = x * answer;
 425     else if (i == 4) { xubar = sumUbarDbar / (1. + answer);
 426       xdbar = sumUbarDbar * answer / (1. + answer); }
 427     else if (i == 5) xs = x * answer;
 428     else if (i == 6) xc = x * answer;
 429     else if (i == 7) xb = x * answer;
 430   }
 431
 432   // Subdivision of valence and sea.
 433   xuVal = xu - xubar;
 434   xuSea = xubar;
 435   xdVal = xd - xdbar;
 436   xdSea = xdbar;
 437
 438   // idSav = 9 to indicate that all flavours reset. id change is dummy here.
 439   idSav = 9;
 440   id    = 0;
 441
 442 }
 443
 444 //**************************************************************************
 445
 446 // Interface to the LHAPDF library.
 447
 448 //*********
 449
 450 // Definitions of static variables.
 451
 452 string LHAPDF::latestSetName = " ";
 453 int    LHAPDF::latestMember  = -1;
 454 int    LHAPDF::latestNSet    = 0;
 455
 456 //*********
 457
 458 // Initialize a parton density function from LHAPDF.
 459
 460 void LHAPDF::init(string setName, int member, Info* infoPtr) {
 461
 462   // If already initialized then need not do anything.
 463   if (setName == latestSetName && member == latestMember
 464     && nSet == latestNSet) return;
 465
 466   // Initialize set. If first character is '/' then assume that name
 467   // is given with path, else not.
 468   if (setName[0] == '/') LHAPDFInterface::initPDFsetM( nSet, setName);
 469   else LHAPDFInterface::initPDFsetByNameM( nSet, setName);
 470
 471   // Check that not dummy library was linked and put nSet negative.
 472   isSet = (nSet >= 0);
 473   if (!isSet) {
 474     if (infoPtr > 0) infoPtr->errorMsg("Error from LHAPDF::init: "
 475       "you try to use LHAPDF but did not link it");
 476     else cout << "Error from LHAPDF::init: you try to use LHAPDF "
 477       << "but did not link it" << endl;
 478   }
 479
 480   // Initialize member.
 481   LHAPDFInterface::initPDFM(nSet, member);
 482
 483   // Do not collect statistics on under/overflow to save time and space.
 484    LHAPDFInterface::setPDFparm( "NOSTAT" );
 485    LHAPDFInterface::setPDFparm( "LOWKEY" );
 486
 487   // Save values to avoid unnecessary reinitializations.
 488   latestSetName = setName;
 489   latestMember  = member;
 490   latestNSet    = nSet;
 491
 492 }
 493
 494 //*********
 495
 496 // Allow optional extrapolation beyond boundaries.
 497
 498 void LHAPDF::setExtrapolate(bool extrapol) {
 499
 500    LHAPDFInterface::setPDFparm( (extrapol) ? "EXTRAPOLATE" : "18" );
 501
 502 }
 503
 504 //*********
 505
 506 // Give the parton distribution function set from LHAPDF.
 507
 508 void LHAPDF::xfUpdate(int , double x, double Q2) {
 509
 510   // Let LHAPDF do the evaluation of parton densities.
 511   double Q = sqrt( max( 0., Q2));
 512   LHAPDFInterface::evolvePDFM( nSet, x, Q, xfArray);
 513
 514   // Update values.
 515   xg    = xfArray[6];
 516   xu    = xfArray[8];
 517   xd    = xfArray[7];
 518   xubar = xfArray[4];
 519   xdbar = xfArray[5];
 520   xs    = xfArray[9];
 521   xc    = xfArray[10];
 522   xb    = xfArray[11];
 523
 524   // Subdivision of valence and sea.
 525   xuVal = xu - xubar;
 526   xuSea = xubar;
 527   xdVal = xd - xdbar;
 528   xdSea = xdbar;
 529
 530   // idSav = 9 to indicate that all flavours reset.
 531   idSav = 9;
 532
 533 }
 534
 535 //**************************************************************************
 536
 537 // Gives electron (or muon, or tau) parton distribution.
 538
 539 // Constant. alphaEM(0) could be taken from settings but safer this way.
 540 const double Lepton::ALPHAEM = 0.00729735;
 541
 542 void Lepton::xfUpdate(int id, double x, double Q2) {
 543
 544   // Squared mass of lepton species: electron, muon, tau.
 545   if (!isInit) {
 546     m2Lep = pow2( ParticleDataTable::m0(idBeam) );
 547     isInit = true;
 548   }
 549
 550   // Electron inside electron, see R. Kleiss et al., in Z physics at
 551   // LEP 1, CERN 89-08, p. 34
 552   double xLog = log(max(1e-10,x));
 553   double xMinusLog = log( max(1e-10, 1. - x) );
 554   double Q2Log = log( max(3., Q2/m2Lep) );
 555   double beta = (ALPHAEM / M_PI) * (Q2Log - 1.);
 556   double delta = 1. + (ALPHAEM / M_PI) * (1.5 * Q2Log + 1.289868)
 557     + pow2(ALPHAEM / M_PI) * (-2.164868 * Q2Log*Q2Log
 558     + 9.840808 * Q2Log - 10.130464);
 559   double fPrel =  beta * pow(1. - x, beta - 1.) * sqrtpos( delta )
 560      - 0.5 * beta * (1. + x) + 0.125 * beta*beta * ( (1. + x)
 561      * (-4. * xMinusLog + 3. * xLog) - 4. * xLog / (1. - x) - 5. - x);
 562
 563   // Zero distribution for very large x and rescale it for intermediate.
 564   if (x > 1. - 1e-10) fPrel = 0.;
 565   else if (x > 1. - 1e-7) fPrel *= pow(1000.,beta) / (pow(1000.,beta) - 1.);
 566   xlepton = x * fPrel;
 567
 568   // Photon inside electron (one possible scheme - primitive).
 569   xgamma = (0.5 * ALPHAEM / M_PI) * Q2Log * (1. + pow2(1. - x));
 570
 571   // idSav = 9 to indicate that all flavours reset. id change is dummy here.
 572   idSav = 9;
 573   id    = 0;
 574
 575 }
 576
 577 //**************************************************************************
 578
 579 } // end namespace Pythia8