using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / src / StandardModel.cxx
CommitLineData
5ad4eb21 1// StandardModel.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 AlphaStrong class.
7
8#include "StandardModel.h"
9
10namespace Pythia8 {
11
12//**************************************************************************
13
14// The AlphaStrong class.
15
16//*********
17
18// Constants: could be changed here if desired, but normally should not.
19// These are of technical nature, as described for each.
20
21// Number of iterations to determine Lambda from given alpha_s.
22const int AlphaStrong::NITER = 10;
23
24// Always evaluate running alpha_s above Lambda3 to avoid disaster.
25// Safety margin picked to freeze roughly for alpha_s = 10.
26const double AlphaStrong::SAFETYMARGIN1 = 1.07;
27const double AlphaStrong::SAFETYMARGIN2 = 1.33;
28
29//*********
30
31// Initialize alpha_strong calculation by finding Lambda values etc.
32
33void AlphaStrong::init( double valueIn, int orderIn) {
34
35 // Order of alpha_s evaluation. Charm, bottom and Z0 masses.
36 // Pick defaults if ParticleDataTable not properly initialized.
37 valueRef = valueIn;
38 order = max( 0, min( 2, orderIn ) );
39 mc = ParticleDataTable::m0(4);
40 if (mc < 1. || mc > 2.) mc = 1.5;
41 mb = ParticleDataTable::m0(5);
42 if (mb < 3. || mb > 6.) mb = 4.8;
43 mZ = ParticleDataTable::m0(23);
44 if (mZ < 90. || mZ > 92.) mZ = 91.188;
45
46 // Fix alpha_s.
47 if (order == 0) {
48 Lambda3Save = Lambda4Save = Lambda5Save = scale2Min = 0.;
49
50 // First order alpha_s: match at flavour thresholds.
51 } else if (order == 1) {
52 Lambda5Save = mZ * exp( -6. * M_PI / (23. * valueRef) );
53 Lambda4Save = Lambda5Save * pow(mb/Lambda5Save, 2./25.);
54 Lambda3Save = Lambda4Save * pow(mc/Lambda4Save, 2./27.);
55 scale2Min = pow2(SAFETYMARGIN1 * Lambda3Save);
56
57 // Second order alpha_s: iterative match at flavour thresholds.
58 } else {
59 double b15 = 348. / 529.;
60 double b14 = 462. / 625.;
61 double b13 = 64. / 81.;
62 double b25 = 224687. / 242208.;
63 double b24 = 548575. / 426888.;
64 double b23 = 938709. / 663552.;
65 double logScale, loglogScale, correction, valueIter;
66
67 // Find Lambda_5 at m_Z.
68 Lambda5Save = mZ * exp( -6. * M_PI / (23. * valueRef) );
69 for (int iter = 0; iter < NITER; ++iter) {
70 logScale = 2. * log(mZ/Lambda5Save);
71 loglogScale = log(logScale);
72 correction = 1. - b15 * loglogScale / logScale
73 + pow2(b15 / logScale) * (pow2(loglogScale - 0.5) + b25 - 1.25);
74 valueIter = valueRef / correction;
75 Lambda5Save = mZ * exp( -6. * M_PI / (23. * valueIter) );
76 }
77
78 // Find Lambda_4 at m_b.
79 double logScaleB = 2. * log(mb/Lambda5Save);
80 double loglogScaleB = log(logScaleB);
81 double valueB = 12. * M_PI / (23. * logScaleB)
82 * (1. - b15 * loglogScaleB / logScaleB
83 + pow2(b15 / logScaleB) * (pow2(loglogScaleB - 0.5) + b25- 1.25) );
84 Lambda4Save = Lambda5Save;
85 for (int iter = 0; iter < NITER; ++iter) {
86 logScale = 2. * log(mb/Lambda4Save);
87 loglogScale = log(logScale);
88 correction = 1. - b14 * loglogScale / logScale
89 + pow2(b14 / logScale) * (pow2(loglogScale - 0.5) + b24 - 1.25);
90 valueIter = valueB / correction;
91 Lambda4Save = mb * exp( -6. * M_PI / (25. * valueIter) );
92 }
93
94 // Find Lambda_3 at m_c.
95 double logScaleC = 2. * log(mc/Lambda4Save);
96 double loglogScaleC = log(logScaleC);
97 double valueC = 12. * M_PI / (25. * logScaleC)
98 * (1. - b14 * loglogScaleC / logScaleC
99 + pow2(b14 / logScaleC) * (pow2(loglogScaleC - 0.5) + b24 - 1.25) );
100 Lambda3Save = Lambda4Save;
101 for (int iter = 0; iter < NITER; ++iter) {
102 logScale = 2. * log(mc/Lambda3Save);
103 loglogScale = log(logScale);
104 correction = 1. - b13 * loglogScale / logScale
105 + pow2(b13 / logScale) * (pow2(loglogScale - 0.5) + b23 - 1.25);
106 valueIter = valueC / correction;
107 Lambda3Save = mc * exp( -6. * M_PI / (27. * valueIter) );
108 }
109 scale2Min = pow2(SAFETYMARGIN2 * Lambda3Save);
110 }
111
112 // Save squares of mass and Lambda values as well.
113 mc2 = pow2(mc);
114 mb2 = pow2(mb);
115 Lambda3Save2 = pow2(Lambda3Save);
116 Lambda4Save2 = pow2(Lambda4Save);
117 Lambda5Save2 = pow2(Lambda5Save);
118 valueNow = valueIn;
119 scale2Now = mZ * mZ;
120 isInit = true;
121
122}
123
124//*********
125
126// Calculate alpha_s value
127
128double AlphaStrong::alphaS( double scale2) {
129
130 // Check for initialization and ensure minimal scale2 value.
131 if (!isInit) return 0.;
132 if (scale2 < scale2Min) scale2 = scale2Min;
133
134 // If equal to old scale then same answer.
135 if (scale2 == scale2Now && (order < 2 || lastCallToFull)) return valueNow;
136 scale2Now = scale2;
137 lastCallToFull = true;
138
139 // Fix alpha_s.
140 if (order == 0) {
141 valueNow = valueRef;
142
143 // First order alpha_s: differs by mass region.
144 } else if (order == 1) {
145 if (scale2 > mb2)
146 valueNow = 12. * M_PI / (23. * log(scale2/Lambda5Save2));
147 else if (scale2 > mc2)
148 valueNow = 12. * M_PI / (25. * log(scale2/Lambda4Save2));
149 else valueNow = 12. * M_PI / (27. * log(scale2/Lambda3Save2));
150
151 // Second order alpha_s: differs by mass region.
152 } else {
153 double Lambda2, b0, b1, b2;
154 if (scale2 > mb2) {
155 Lambda2 = Lambda5Save2;
156 b0 = 23.;
157 b1 = 348. / 529.;
158 b2 = 224687. / 242208.;
159 } else if (scale2 > mc2) {
160 Lambda2 = Lambda4Save2;
161 b0 = 25.;
162 b1 = 462. / 625.;
163 b2 = 548575. / 426888.;
164 } else {
165 Lambda2 = Lambda3Save2;
166 b0 = 27.;
167 b1 = 64. / 81.;
168 b2 = 938709. / 663552.;
169 }
170 double logScale = log(scale2/Lambda2);
171 double loglogScale = log(logScale);
172 valueNow = 12. * M_PI / (b0 * logScale)
173 * ( 1. - b1 * loglogScale / logScale
174 + pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25) );
175 }
176
177 // Done.
178 return valueNow;
179
180}
181
182//*********
183
184// Calculate alpha_s value, but only use up to first-order piece.
185// (To be combined with alphaS2OrdCorr.)
186
187double AlphaStrong::alphaS1Ord( double scale2) {
188
189 // Check for initialization and ensure minimal scale2 value.
190 if (!isInit) return 0.;
191 if (scale2 < scale2Min) scale2 = scale2Min;
192
193 // If equal to old scale then same answer.
194 if (scale2 == scale2Now && (order < 2 || !lastCallToFull)) return valueNow;
195 scale2Now = scale2;
196 lastCallToFull = false;
197
198 // Fix alpha_S.
199 if (order == 0) {
200 valueNow = valueRef;
201
202 // First/second order alpha_s: differs by mass region.
203 } else {
204 if (scale2 > mb2)
205 valueNow = 12. * M_PI / (23. * log(scale2/Lambda5Save2));
206 else if (scale2 > mc2)
207 valueNow = 12. * M_PI / (25. * log(scale2/Lambda4Save2));
208 else valueNow = 12. * M_PI / (27. * log(scale2/Lambda3Save2));
209 }
210
211 // Done.
212 return valueNow;
213}
214
215//*********
216
217// Calculates the second-order extra factor in alpha_s.
218// (To be combined with alphaS1Ord.)
219
220double AlphaStrong::alphaS2OrdCorr( double scale2) {
221
222 // Check for initialization and ensure minimal scale2 value.
223 if (!isInit) return 1.;
224 if (scale2 < scale2Min) scale2 = scale2Min;
225
226 // Only meaningful for second order calculations.
227 if (order < 2) return 1.;
228
229 // Second order correction term: differs by mass region.
230 double Lambda2, b0, b1, b2;
231 if (scale2 > mb2) {
232 Lambda2 = Lambda5Save2;
233 b0 = 23.;
234 b1 = 348. / 529.;
235 b2 = 224687. / 242208.;
236 } else if (scale2 > mc2) {
237 Lambda2 = Lambda4Save2;
238 b0 = 25.;
239 b1 = 462. / 625.;
240 b2 = 548575. / 426888.;
241 } else {
242 Lambda2 = Lambda3Save2;
243 b0 = 27.;
244 b1 = 64. / 81.;
245 b2 = 938709. / 663552.;
246 }
247 double logScale = log(scale2/Lambda2);
248 double loglogScale = log(logScale);
249 return ( 1. - b1 * loglogScale / logScale
250 + pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25) );
251
252}
253
254//**************************************************************************
255
256// The AlphaEM class.
257
258//*********
259
260// Definitions of static variables.
261double AlphaEM::alpEM0 = 0.00729735;
262double AlphaEM::alpEMmZ = 0.00781751;
263double AlphaEM::mZ2 = 8315.;
264// Effective thresholds for electron, muon, light quarks, tau+c, b.
265double AlphaEM::Q2step[5] = {0.26e-6, 0.011, 0.25, 3.5, 90.};
266// Running coefficients are sum charge2 / 3 pi in pure QED, here slightly
267// enhanced for quarks to approximately account for QCD corrections.
268double AlphaEM::bRun[5] = {0.1061, 0.2122, 0.460, 0.700, 0.725};
269double AlphaEM::alpEMstep[5] = {};
270
271//*********
272
273// Initialize alpha_EM calculation.
274
275void AlphaEM::initStatic() {
276
277 // Read in alpha_EM value at 0 and m_Z, and mass of Z.
278 alpEM0 = Settings::parm("StandardModel:alphaEM0");
279 alpEMmZ = Settings::parm("StandardModel:alphaEMmZ");
280 double mZ = ParticleDataTable::m0(23);
281 mZ2 = mZ * mZ;
282
283 // AlphaEM values at matching scales and matching b value.
284
285 // Step down from mZ to tau/charm threshold.
286 alpEMstep[4] = alpEMmZ / ( 1. + alpEMmZ * bRun[4]
287 * log(mZ2 / Q2step[4]) );
288 alpEMstep[3] = alpEMstep[4] / ( 1. - alpEMstep[4] * bRun[3]
289 * log(Q2step[3] / Q2step[4]) );
290
291 // Step up from me to light-quark threshold.
292 alpEMstep[0] = alpEM0;
293 alpEMstep[1] = alpEMstep[0] / ( 1. - alpEMstep[0] * bRun[0]
294 * log(Q2step[1] / Q2step[0]) );
295 alpEMstep[2] = alpEMstep[1] / ( 1. - alpEMstep[1] * bRun[1]
296 * log(Q2step[2] / Q2step[1]) );
297
298 // Fit b in range between light-quark and tau/charm to join smoothly.
299 bRun[2] = (1./alpEMstep[3] - 1./alpEMstep[2])
300 / log(Q2step[2] / Q2step[3]);
301
302}
303
304//*********
305
306// Calculate alpha_EM value
307
308double AlphaEM::alphaEM( double scale2) {
309
310 // Fix alphaEM; for order = -1 fixed at m_Z.
311 if (order == 0) return alpEM0;
312 if (order < 0) return alpEMmZ;
313
314 // Running alphaEM.
315 for (int i = 4; i >= 0; --i) if (scale2 > Q2step[i])
316 return alpEMstep[i] / (1. - bRun[i] * alpEMstep[i]
317 * log(scale2 / Q2step[i]) );
318 return alpEM0;
319
320}
321
322//**************************************************************************
323
324// The CoupEW class.
325
326//*********
327
328// Definitions of static variables.
329
330double CoupEW::s2tW = 0.232;
331double CoupEW::c2tW = 0.768;
332double CoupEW::s2tWbar = 0.232;
333double CoupEW::efSave[20] = { 0., -1./3., 2./3., -1./3., 2./3., -1./3.,
334 2./3., -1./3., 2./3., 0., 0., -1., 0., -1., 0., -1., 0., -1., 0., 0.};
335double CoupEW::vfSave[20] = { };
336double CoupEW::afSave[20] = { 0., -1., 1., -1., 1., -1., 1., -1., 1.,
337 0., 0., -1., 1., -1., 1., -1., 1., -1., 1., 0.};
338double CoupEW::lfSave[20] = { };
339double CoupEW::rfSave[20] = { };
340double CoupEW::ef2Save[20] = { };
341double CoupEW::vf2Save[20] = { };
342double CoupEW::af2Save[20] = { };
343double CoupEW::efvfSave[20] = { };
344double CoupEW::vf2af2Save[20] = { };
345
346//*********
347
348// Initialize electroweak mixing angle and couplings.
349
350void CoupEW::initStatic() {
351
352 // Read in electroweak mixing angle.
353 s2tW = Settings::parm("StandardModel:sin2thetaW");
354 c2tW = 1. - s2tW;
355 s2tWbar = Settings::parm("StandardModel:sin2thetaWbar");
356
357 // Initialize electroweak couplings.
358 for (int i = 0; i < 20; ++i) {
359 vfSave[i] = afSave[i] - 4. * s2tWbar * efSave[i];
360 lfSave[i] = afSave[i] - 2. * s2tWbar * efSave[i];
361 rfSave[i] = - 2. * s2tWbar * efSave[i];
362 ef2Save[i] = pow2(efSave[i]);
363 vf2Save[i] = pow2(vfSave[i]);
364 af2Save[i] = pow2(afSave[i]);
365 efvfSave[i] = efSave[i] * vfSave[i];
366 vf2af2Save[i] = vf2Save[i] + af2Save[i];
367 }
368
369}
370
371//**************************************************************************
372
373// The VCKM class.
374
375//*********
376
377// Definitions of static variables. Initialize to all elements zero.
378
379double VCKM::Vsave[5][5] = { };
380double VCKM::V2save[5][5] = { };
381double VCKM::V2out[20] = { };
382
383//*********
384
385// Prepare to handle CKM matrix elements.
386
387void VCKM::initStatic() {
388
389 // Read in matrix element values and store them.
390 Vsave[1][1] = Settings::parm("StandardModel:Vud");
391 Vsave[1][2] = Settings::parm("StandardModel:Vus");
392 Vsave[1][3] = Settings::parm("StandardModel:Vub");
393 Vsave[2][1] = Settings::parm("StandardModel:Vcd");
394 Vsave[2][2] = Settings::parm("StandardModel:Vcs");
395 Vsave[2][3] = Settings::parm("StandardModel:Vcb");
396 Vsave[3][1] = Settings::parm("StandardModel:Vtd");
397 Vsave[3][2] = Settings::parm("StandardModel:Vts");
398 Vsave[3][3] = Settings::parm("StandardModel:Vtb");
399
400 // Also allow for the potential existence of a fourth generation.
401 Vsave[1][4] = Settings::parm("FourthGeneration:VubPrime");
402 Vsave[2][4] = Settings::parm("FourthGeneration:VcbPrime");
403 Vsave[3][4] = Settings::parm("FourthGeneration:VtbPrime");
404 Vsave[4][1] = Settings::parm("FourthGeneration:VtPrimed");
405 Vsave[4][2] = Settings::parm("FourthGeneration:VtPrimes");
406 Vsave[4][3] = Settings::parm("FourthGeneration:VtPrimeb");
407 Vsave[4][4] = Settings::parm("FourthGeneration:VtPrimebPrime");
408
409 // Calculate squares of matrix elements.
410 for(int i = 1; i < 5; ++i) for(int j = 1; j < 5; ++j)
411 V2save[i][j] = pow2(Vsave[i][j]);
412
413 // Sum VCKM^2_out sum for given incoming flavour, excluding top as partner.
414 V2out[1] = V2save[1][1] + V2save[2][1];
415 V2out[2] = V2save[1][1] + V2save[1][2] + V2save[1][3];
416 V2out[3] = V2save[1][2] + V2save[2][2];
417 V2out[4] = V2save[2][1] + V2save[2][2] + V2save[2][3];
418 V2out[5] = V2save[1][3] + V2save[2][3];
419 V2out[6] = V2save[3][1] + V2save[3][2] + V2save[3][3];
420 V2out[7] = V2save[1][4] + V2save[2][4];
421 V2out[8] = V2save[4][1] + V2save[4][2] + V2save[4][3];
422 for (int i = 11; i <= 18; ++i) V2out[i] = 1.;
423
424}
425
426//*********
427
428// Return CKM value for incoming flavours (sign irrelevant).
429
430double VCKM::Vid(int id1, int id2) {
431
432 // Use absolute sign (want to cover both f -> f' W and f fbar' -> W).
433 int id1Abs = abs(id1);
434 int id2Abs = abs(id2);
435 if (id1Abs == 0 || id2Abs == 0 || (id1Abs + id2Abs)%2 != 1) return 0.;
436
437 // Ensure proper order before reading out from Vsave or lepton match.
438 if (id1Abs%2 == 1) swap(id1Abs, id2Abs);
439 if (id1Abs <= 8 && id2Abs <= 8) return Vsave[id1Abs/2][(id2Abs + 1)/2];
440 if ( (id1Abs == 12 || id1Abs == 14 || id1Abs == 16 || id1Abs == 18)
441 && id2Abs == id1Abs - 1 ) return 1.;
442
443 // No more valid cases.
444 return 0.;
445
446}
447
448//*********
449
450// Return squared CKM value for incoming flavours (sign irrelevant).
451
452double VCKM::V2id(int id1, int id2) {
453
454 // Use absolute sign (want to cover both f -> f' W and f fbar' -> W).
455 int id1Abs = abs(id1);
456 int id2Abs = abs(id2);
457 if (id1Abs == 0 || id2Abs == 0 || (id1Abs + id2Abs)%2 != 1) return 0.;
458
459 // Ensure proper order before reading out from V2save or lepton match.
460 if (id1Abs%2 == 1) swap(id1Abs, id2Abs);
461 if (id1Abs <= 8 && id2Abs <= 8) return V2save[id1Abs/2][(id2Abs + 1)/2];
462 if ( (id1Abs == 12 || id1Abs == 14 || id1Abs == 16 || id1Abs == 18)
463 && id2Abs == id1Abs - 1 ) return 1.;
464
465 // No more valid cases.
466 return 0.;
467
468}
469
470//*********
471
472// Pick an outgoing flavour for given incoming one, given CKM mixing.
473
474int VCKM::V2pick(int id) {
475
476 // Initial values.
477 int idIn = abs(id);
478 int idOut = 0;
479
480 // Quarks: need to make random choice.
481 if (idIn >= 1 && idIn <= 8) {
482 double V2rndm = Rndm::flat() * V2out[idIn];
483 if (idIn == 1) idOut = (V2rndm < V2save[1][1]) ? 2 : 4;
484 else if (idIn == 2) idOut = (V2rndm < V2save[1][1]) ? 1
485 : ( (V2rndm < V2save[1][1] + V2save[1][2]) ? 3 : 5 );
486 else if (idIn == 3) idOut = (V2rndm < V2save[1][2]) ? 2 : 4;
487 else if (idIn == 4) idOut = (V2rndm < V2save[2][1]) ? 1
488 : ( (V2rndm < V2save[2][1] + V2save[2][2]) ? 3 : 5 );
489 else if (idIn == 5) idOut = (V2rndm < V2save[1][3]) ? 2 : 4;
490 else if (idIn == 6) idOut = (V2rndm < V2save[3][1]) ? 1
491 : ( (V2rndm < V2save[3][1] + V2save[3][2]) ? 3 : 5 );
492 else if (idIn == 7) idOut = (V2rndm < V2save[1][4]) ? 2 : 4;
493 else if (idIn == 8) idOut = (V2rndm < V2save[4][1]) ? 1
494 : ( (V2rndm < V2save[4][1] + V2save[4][2]) ? 3 : 5 );
495
496 // Leptons: unambiguous.
497 } else if (idIn >= 11 && idIn <= 18) {
498 if (idIn%2 == 1) idOut = idIn + 1;
499 else idOut = idIn - 1;
500 }
501
502 // Done. Return with sign.
503 return ( (id > 0) ? idOut : -idOut );
504
505}
506
507//**************************************************************************
508
509} // end namespace Pythia8