source: HiSusy/trunk/Pythia8/pythia8170/src/StandardModel.cc @ 1

Last change on this file since 1 was 1, checked in by zerwas, 11 years ago

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