stablejet is hosted by Hepforge, IPPP Durham
StableJet
StandardModel.cc
1 // StandardModel.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2011 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 #include <cmath>
10 #include <algorithm>
11 
12 inline static double pow2(double x) { return pow(x,2); }
13 
14 
15 
16 //==========================================================================
17 
18 // The AlphaStrong class.
19 
20 //--------------------------------------------------------------------------
21 
22 // Constants: could be changed here if desired, but normally should not.
23 // These are of technical nature, as described for each.
24 
25 // Number of iterations to determine Lambda from given alpha_s.
26 const int AlphaStrong::NITER = 10;
27 
28 // Masses: m_c, m_b, m_Z. Used for flavour thresholds and normalization scale.
29 const double AlphaStrong::MC = 1.5;
30 const double AlphaStrong::MB = 4.8;
31 const double AlphaStrong::MZ = 91.188;
32 
33 // Always evaluate running alpha_s above Lambda3 to avoid disaster.
34 // Safety margin picked to freeze roughly for alpha_s = 10.
35 const double AlphaStrong::SAFETYMARGIN1 = 1.07;
36 const double AlphaStrong::SAFETYMARGIN2 = 1.33;
37 
38 
39 
40 //--------------------------------------------------------------------------
41 
42 // Initialize alpha_strong calculation by finding Lambda values etc.
43 
44 void AlphaStrong::init( double valueIn, int orderIn) {
45 
46  // Order of alpha_s evaluation.
47  valueRef = valueIn;
48  order = std::max( 0, std::min( 2, orderIn ) );
49 
50  // Fix alpha_s.
51  if (order == 0) {
52  Lambda3Save = Lambda4Save = Lambda5Save = scale2Min = 0.;
53 
54  // First order alpha_s: match at flavour thresholds.
55  } else if (order == 1) {
56  Lambda5Save = MZ * exp( -6. * M_PI / (23. * valueRef) );
57  Lambda4Save = Lambda5Save * pow(MB/Lambda5Save, 2./25.);
58  Lambda3Save = Lambda4Save * pow(MC/Lambda4Save, 2./27.);
59  scale2Min = pow2(SAFETYMARGIN1 * Lambda3Save);
60 
61  // Second order alpha_s: iterative match at flavour thresholds.
62  } else {
63  double b15 = 348. / 529.;
64  double b14 = 462. / 625.;
65  double b13 = 64. / 81.;
66  double b25 = 224687. / 242208.;
67  double b24 = 548575. / 426888.;
68  double b23 = 938709. / 663552.;
69  double logScale, loglogScale, correction, valueIter;
70 
71  // Find Lambda_5 at m_Z.
72  Lambda5Save = MZ * exp( -6. * M_PI / (23. * valueRef) );
73  for (int iter = 0; iter < NITER; ++iter) {
74  logScale = 2. * log(MZ/Lambda5Save);
75  loglogScale = log(logScale);
76  correction = 1. - b15 * loglogScale / logScale
77  + pow2(b15 / logScale) * (pow2(loglogScale - 0.5) + b25 - 1.25);
78  valueIter = valueRef / correction;
79  Lambda5Save = MZ * exp( -6. * M_PI / (23. * valueIter) );
80  }
81 
82  // Find Lambda_4 at m_b.
83  double logScaleB = 2. * log(MB/Lambda5Save);
84  double loglogScaleB = log(logScaleB);
85  double valueB = 12. * M_PI / (23. * logScaleB)
86  * (1. - b15 * loglogScaleB / logScaleB
87  + pow2(b15 / logScaleB) * (pow2(loglogScaleB - 0.5) + b25- 1.25) );
88  Lambda4Save = Lambda5Save;
89  for (int iter = 0; iter < NITER; ++iter) {
90  logScale = 2. * log(MB/Lambda4Save);
91  loglogScale = log(logScale);
92  correction = 1. - b14 * loglogScale / logScale
93  + pow2(b14 / logScale) * (pow2(loglogScale - 0.5) + b24 - 1.25);
94  valueIter = valueB / correction;
95  Lambda4Save = MB * exp( -6. * M_PI / (25. * valueIter) );
96  }
97 
98  // Find Lambda_3 at m_c.
99  double logScaleC = 2. * log(MC/Lambda4Save);
100  double loglogScaleC = log(logScaleC);
101  double valueC = 12. * M_PI / (25. * logScaleC)
102  * (1. - b14 * loglogScaleC / logScaleC
103  + pow2(b14 / logScaleC) * (pow2(loglogScaleC - 0.5) + b24 - 1.25) );
104  Lambda3Save = Lambda4Save;
105  for (int iter = 0; iter < NITER; ++iter) {
106  logScale = 2. * log(MC/Lambda3Save);
107  loglogScale = log(logScale);
108  correction = 1. - b13 * loglogScale / logScale
109  + pow2(b13 / logScale) * (pow2(loglogScale - 0.5) + b23 - 1.25);
110  valueIter = valueC / correction;
111  Lambda3Save = MC * exp( -6. * M_PI / (27. * valueIter) );
112  }
113  scale2Min = pow2(SAFETYMARGIN2 * Lambda3Save);
114  }
115 
116  // Save squares of mass and Lambda values as well.
117  mc2 = pow2(MC);
118  mb2 = pow2(MB);
119  Lambda3Save2 = pow2(Lambda3Save);
120  Lambda4Save2 = pow2(Lambda4Save);
121  Lambda5Save2 = pow2(Lambda5Save);
122  valueNow = valueIn;
123  scale2Now = MZ * MZ;
124  isInit = true;
125 
126 }
127 
128 //--------------------------------------------------------------------------
129 
130 // Calculate alpha_s value
131 
132 double AlphaStrong::alphaS( double scale2) const {
133 
134  // Check for initialization and ensure minimal scale2 value.
135  if (!isInit) return 0.;
136  if (scale2 < scale2Min) scale2 = scale2Min;
137 
138  // If equal to old scale then same answer.
139  if (scale2 == scale2Now && (order < 2 || lastCallToFull)) return valueNow;
140  scale2Now = scale2;
141  lastCallToFull = true;
142 
143  // Fix alpha_s.
144  if (order == 0) {
145  valueNow = valueRef;
146 
147  // First order alpha_s: differs by mass region.
148  } else if (order == 1) {
149  if (scale2 > mb2)
150  valueNow = 12. * M_PI / (23. * log(scale2/Lambda5Save2));
151  else if (scale2 > mc2)
152  valueNow = 12. * M_PI / (25. * log(scale2/Lambda4Save2));
153  else valueNow = 12. * M_PI / (27. * log(scale2/Lambda3Save2));
154 
155  // Second order alpha_s: differs by mass region.
156  } else {
157  double Lambda2, b0, b1, b2;
158  if (scale2 > mb2) {
159  Lambda2 = Lambda5Save2;
160  b0 = 23.;
161  b1 = 348. / 529.;
162  b2 = 224687. / 242208.;
163  } else if (scale2 > mc2) {
164  Lambda2 = Lambda4Save2;
165  b0 = 25.;
166  b1 = 462. / 625.;
167  b2 = 548575. / 426888.;
168  } else {
169  Lambda2 = Lambda3Save2;
170  b0 = 27.;
171  b1 = 64. / 81.;
172  b2 = 938709. / 663552.;
173  }
174  double logScale = log(scale2/Lambda2);
175  double loglogScale = log(logScale);
176  valueNow = 12. * M_PI / (b0 * logScale)
177  * ( 1. - b1 * loglogScale / logScale
178  + pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25) );
179  }
180 
181  // Done.
182  return valueNow;
183 
184 }
185 
186 //--------------------------------------------------------------------------
187 
188 // Calculate alpha_s value, but only use up to first-order piece.
189 // (To be combined with alphaS2OrdCorr.)
190 
191 double AlphaStrong::alphaS1Ord( double scale2) {
192 
193  // Check for initialization and ensure minimal scale2 value.
194  if (!isInit) return 0.;
195  if (scale2 < scale2Min) scale2 = scale2Min;
196 
197  // If equal to old scale then same answer.
198  if (scale2 == scale2Now && (order < 2 || !lastCallToFull)) return valueNow;
199  scale2Now = scale2;
200  lastCallToFull = false;
201 
202  // Fix alpha_S.
203  if (order == 0) {
204  valueNow = valueRef;
205 
206  // First/second order alpha_s: differs by mass region.
207  } else {
208  if (scale2 > mb2)
209  valueNow = 12. * M_PI / (23. * log(scale2/Lambda5Save2));
210  else if (scale2 > mc2)
211  valueNow = 12. * M_PI / (25. * log(scale2/Lambda4Save2));
212  else valueNow = 12. * M_PI / (27. * log(scale2/Lambda3Save2));
213  }
214 
215  // Done.
216  return valueNow;
217 }
218 
219 //--------------------------------------------------------------------------
220 
221 // Calculates the second-order extra factor in alpha_s.
222 // (To be combined with alphaS1Ord.)
223 
224 double AlphaStrong::alphaS2OrdCorr( double scale2) {
225 
226  // Check for initialization and ensure minimal scale2 value.
227  if (!isInit) return 1.;
228  if (scale2 < scale2Min) scale2 = scale2Min;
229 
230  // Only meaningful for second order calculations.
231  if (order < 2) return 1.;
232 
233  // Second order correction term: differs by mass region.
234  double Lambda2, b0, b1, b2;
235  if (scale2 > mb2) {
236  Lambda2 = Lambda5Save2;
237  b0 = 23.;
238  b1 = 348. / 529.;
239  b2 = 224687. / 242208.;
240  } else if (scale2 > mc2) {
241  Lambda2 = Lambda4Save2;
242  b0 = 25.;
243  b1 = 462. / 625.;
244  b2 = 548575. / 426888.;
245  } else {
246  Lambda2 = Lambda3Save2;
247  b0 = 27.;
248  b1 = 64. / 81.;
249  b2 = 938709. / 663552.;
250  }
251  double logScale = log(scale2/Lambda2);
252  double loglogScale = log(logScale);
253  return ( 1. - b1 * loglogScale / logScale
254  + pow2(b1 / logScale) * (pow2(loglogScale - 0.5) + b2 - 1.25) );
255 
256 }
257 
258 
Generated by   doxygen 1.8.6