using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / src / PartonDistributions.cxx
CommitLineData
5ad4eb21 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
12namespace Pythia8 {
13
14//**************************************************************************
15
16// Base class for parton distribution functions.
17
18//*********
19
20// Standard parton densities.
21
22double 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
57double 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
84double 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
121void 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
249double 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
260double 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
271double 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
295void 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
452string LHAPDF::latestSetName = " ";
453int LHAPDF::latestMember = -1;
454int LHAPDF::latestNSet = 0;
455
456//*********
457
458// Initialize a parton density function from LHAPDF.
459
460void 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
498void 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
508void 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.
540const double Lepton::ALPHAEM = 0.00729735;
541
542void 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