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