a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
EWSMApproximateFormulae.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2012 HEPfit Collaboration
3  *
4  *
5  * For the licensing terms see doc/COPYING.
6  */
7 
8 #include "EWSMApproximateFormulae.h"
9 #include <stdexcept>
10 #include <sstream>
11 
12 #define UpperBoundForApproximateFormulae 1000.0
13 //#define UpperBoundForApproximateFormulae 1500.0 // for test
14 
15 EWSMApproximateFormulae::EWSMApproximateFormulae(const EWSMcache& cache_i)
16 : mycache(cache_i)
17 {
18 }
19 
20 
22 
23 double EWSMApproximateFormulae::Mw() const
24 {
25  // Parametrization from arXiv:hep-ph/0311148v2 (updates from the journal version)
26  double Mw0, c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11;
27  if (mycache.getSM().getMHl() >= 100.0 && mycache.getSM().getMHl() <= UpperBoundForApproximateFormulae) {
28  // applicable for 100 GeV <= mHl <= 1 TeV
29  Mw0 = 80.3779;
30  c1 = 0.05263;
31  c2 = 0.010239;
32  c3 = 0.000954;
33  c4 = -0.000054;
34  c5 = 1.077;
35  c6 = 0.5252;
36  c7 = 0.0700;
37  c8 = 0.004102;
38  c9 = 0.000111;
39  c10 = 0.0774;
40  c11 = 115.0;
41  } else if (mycache.getSM().getMHl() >= 10.0 && mycache.getSM().getMHl() <= 1000.0) {
42  // applicable for 10 GeV <= mHl <= 1 TeV
43  Mw0 = 80.3799;
44  c1 = 0.05427;
45  c2 = 0.008931;
46  c3 = 0.0000882;
47  c4 = 0.000161;
48  c5 = 1.070;
49  c6 = 0.5237;
50  c7 = 0.0679;
51  c8 = 0.00179;
52  c9 = 0.0000664;
53  c10 = 0.0795;
54  c11 = 114.9;
55  } else {
56  std::stringstream out;
57  out << mycache.getSM().getMHl();
58  throw std::runtime_error("ApproximateFormulae::Mw(): mh=" + out.str() + " is out of range");
59  }
60 
61  double dH = log(mycache.getSM().getMHl() / 100.0);
62  double dh = pow((mycache.getSM().getMHl() / 100.0), 2.0);
63  double dt = pow((mycache.getSM().getMtpole() / 174.3), 2.0) - 1.0;
64  double dZ = mycache.getSM().getMz() / 91.1875 - 1.0;
65  double dalphae = mycache.getSM().DeltaAlphaL5q() / 0.05907 - 1.0;
66  double dalphas = mycache.getSM().getAlsMz() / 0.119 - 1.0;
67 
68  return (Mw0 - c1 * dH - c2 * dH * dH + c3 * pow(dH, 4.0)
69  + c4 * (dh - 1.0) - c5 * dalphae + c6 * dt - c7 * dt * dt
70  - c8 * dH * dt + c9 * dh * dt - c10 * dalphas + c11 * dZ
71  + mycache.getSM().getDelMw());
72 }
73 
74 double EWSMApproximateFormulae::sin2thetaEff_l(const QCD::lepton l) const
75 {
76  // applicable for 10 GeV <= mHl <= 1 TeV
77  if (mycache.getSM().getMHl() < 10.0 || mycache.getSM().getMHl() > UpperBoundForApproximateFormulae) {
78  std::stringstream out;
79  out << mycache.getSM().getMHl();
80  throw std::runtime_error("ApproximateFormulae::sin2thetaEff_l(): mh=" + out.str() + " is out of range");
81  }
82 
83  double s0, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10;
84  switch (l) {
85  case StandardModel::NEUTRINO_1:
86  case StandardModel::NEUTRINO_2:
87  case StandardModel::NEUTRINO_3:
88  s0 = 0.2308772;
89  d1 = 4.713 * 0.0001;
90  d2 = 2.05 * 0.00001;
91  d3 = 3.85 * 0.000001;
92  d4 = -1.85 * 0.000001;
93  d5 = 2.06 * 0.01;
94  d6 = -2.850 * 0.001;
95  d7 = 1.82 * 0.0001;
96  d8 = -9.71 * 0.000001;
97  d9 = 3.96 * 0.0001;
98  d10 = -6.54 * 0.1;
99  break;
100  case StandardModel::ELECTRON:
101  case StandardModel::MU:
102  case StandardModel::TAU:
103  s0 = 0.2312527;
104  d1 = 4.729 * 0.0001;
105  d2 = 2.07 * 0.00001;
106  d3 = 3.85 * 0.000001;
107  d4 = -1.85 * 0.000001;
108  d5 = 2.07 * 0.01;
109  d6 = -2.851 * 0.001;
110  d7 = 1.82 * 0.0001;
111  d8 = -9.74 * 0.000001;
112  d9 = 3.98 * 0.0001;
113  d10 = -6.55 * 0.1;
114  break;
115  default:
116  throw std::runtime_error("Error in ApproximateFormulae::sin2thetaEff_l()");
117  }
118 
119  double L_H = log(mycache.getSM().getMHl() / 100.0);
120  double Delta_H = mycache.getSM().getMHl() / 100.0;
121  double Delta_ale = mycache.getSM().DeltaAlphaL5q() / 0.05907 - 1.0;
122  double Delta_t = pow((mycache.getSM().getMtpole() / 178.0), 2.0) - 1.0;
123  double Delta_alphas = mycache.getSM().getAlsMz() / 0.117 - 1.0;
124  double Delta_Z = mycache.getSM().getMz() / 91.1876 - 1.0;
125 
126  return (s0 + d1 * L_H + d2 * L_H * L_H + d3 * pow(L_H, 4.0)
127  + d4 * (Delta_H * Delta_H - 1.0) + d5 * Delta_ale + d6 * Delta_t
128  + d7 * Delta_t * Delta_t + d8 * Delta_t * (Delta_H - 1.0)
129  + d9 * Delta_alphas + d10 * Delta_Z
130  + mycache.getSM().getDelSin2th_l());
131 }
132 
133 double EWSMApproximateFormulae::sin2thetaEff_q(const QCD::quark q) const
134 {
135  // applicable for 10 GeV <= mHl <= 1 TeV
136  if (mycache.getSM().getMHl() < 10.0 || mycache.getSM().getMHl() > UpperBoundForApproximateFormulae) {
137  std::stringstream out;
138  out << mycache.getSM().getMHl();
139  throw std::runtime_error("ApproximateFormulae::sin2thetaEff_q(): mh=" + out.str() + " is out of range");
140  }
141 
142  double s0, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10;
143  double ThError = 0.0; // Theoretical uncertainty
144 
145  switch (q) {
146  case QCD::UP:
147  case QCD::CHARM:
148  s0 = 0.2311395;
149  d1 = 4.726 * 0.0001;
150  d2 = 2.07 * 0.00001;
151  d3 = 3.85 * 0.000001;
152  d4 = -1.85 * 0.000001;
153  d5 = 2.07 * 0.01;
154  d6 = -2.853 * 0.001;
155  d7 = 1.83 * 0.0001;
156  d8 = -9.73 * 0.000001;
157  d9 = 3.98 * 0.0001;
158  d10 = -6.55 * 0.1;
159  ThError = mycache.getSM().getDelSin2th_q();
160  break;
161  case QCD::DOWN:
162  case QCD::STRANGE:
163  s0 = 0.2310286;
164  d1 = 4.720 * 0.0001;
165  d2 = 2.06 * 0.00001;
166  d3 = 3.85 * 0.000001;
167  d4 = -1.85 * 0.000001;
168  d5 = 2.07 * 0.01;
169  d6 = -2.848 * 0.001;
170  d7 = 1.81 * 0.0001;
171  d8 = -9.73 * 0.000001;
172  d9 = 3.97 * 0.0001;
173  d10 = -6.55 * 0.1;
174  ThError = mycache.getSM().getDelSin2th_q();
175  break;
176  case QCD::BOTTOM:
177 
178  if (mycache.getSM().getMHl() < 120.1 || mycache.getSM().getMHl() > 130.1) {
179  s0 = 0.2327580;
180  d1 = 4.749 * 0.0001;
181  d2 = 2.03 * 0.00001;
182  d3 = 3.94 * 0.000001;
183  d4 = -1.84 * 0.000001;
184  d5 = 2.08 * 0.01;
185  d6 = -0.993 * 0.001;
186  d7 = 0.708 * 0.0001;
187  d8 = -7.61 * 0.000001;
188  d9 = 4.03 * 0.0001;
189  d10 = 6.61 * 0.1;
190  ThError = mycache.getSM().getDelSin2th_b();
191  break;
192 
193  } else {
194  return sin2thetaEff_b();
195  }
196  case QCD::TOP:
197  return 0.0;
198  default:
199  throw std::runtime_error("Error in ApproximateFormulae::sin2thetaEff_q()");
200  }
201 
202  double L_H = log(mycache.getSM().getMHl() / 100.0);
203  double Delta_H = mycache.getSM().getMHl() / 100.0;
204  double Delta_ale = mycache.getSM().DeltaAlphaL5q() / 0.05907 - 1.0;
205  double Delta_t = pow((mycache.getSM().getMtpole() / 178.0), 2.0) - 1.0;
206  double Delta_alphas = mycache.getSM().getAlsMz() / 0.117 - 1.0;
207  double Delta_Z = mycache.getSM().getMz() / 91.1876 - 1.0;
208 
209  return (s0 + d1 * L_H + d2 * L_H * L_H + d3 * pow(L_H, 4.0)
210  + d4 * (Delta_H * Delta_H - 1.0) + d5 * Delta_ale + d6 * Delta_t
211  + d7 * Delta_t * Delta_t + d8 * Delta_t * (Delta_H - 1.0)
212  + d9 * Delta_alphas + d10 * Delta_Z
213  + ThError);
214 }
215 
216 double EWSMApproximateFormulae::sin2thetaEff_b() const
217 {
218  // applicable for 120.1 GeV <= mHl <= 130.1 GeV
219  if (mycache.getSM().getMHl() < 120.1 || mycache.getSM().getMHl() > 130.1) {
220  std::stringstream out;
221  out << mycache.getSM().getMHl();
222  throw std::runtime_error("ApproximateFormulae::sin2thetaEff_b(): mh=" + out.str() + " is out of range");
223  }
224 
225  double s0, d1, d2, d3, d4, d5, d6, d7, d8, d9;
226 
227  s0 = 0.232704;
228  d1 = 4.723 * 0.0001;
229  d2 = 1.97 * 0.0001;
230  d3 = 2.07 * 0.01;
231  d4 = -9.733 * 0.0001;
232  d5 = 3.93 * 0.0001;
233  d6 = -1.38 * 0.0001;
234  d7 = 2.42 * 0.0001;
235  d8 = -8.10 * 0.0001;
236  d9 = -0.664;
237 
238  double L_H = log(mycache.getSM().getMHl() / 125.7);
239  double Delta_ale = mycache.getSM().DeltaAlphaL5q() / 0.059 - 1.0;
240  double Delta_t = pow((mycache.getSM().getMtpole() / 173.2), 2.0) - 1.0;
241  double Delta_alphas = mycache.getSM().getAlsMz() / 0.1184 - 1.0;
242  double Delta_Z = mycache.getSM().getMz() / 91.1876 - 1.0;
243 
244  return (s0 + d1 * L_H + d2 * L_H * L_H
245  + d3 * Delta_ale
246  + d4 * Delta_t + d5 * Delta_t * Delta_t + d6 * Delta_t * L_H
247  + d7 * Delta_alphas + d8 * Delta_t * Delta_alphas
248  + d9 * Delta_Z
249  + mycache.getSM().getDelSin2th_b());
250 }
251 
252 double EWSMApproximateFormulae::DeltaR_TwoLoopEW_rem(const double Mw_i) const
253 {
254  // applicable for 10 GeV <= mHl <= 1 TeV
255  if (mycache.getSM().getMHl() < 10.0 || mycache.getSM().getMHl() > UpperBoundForApproximateFormulae) {
256  std::stringstream out;
257  out << mycache.getSM().getMHl();
258  throw std::runtime_error("ApproximateFormulae::DeltaR_TwoLoopEW_rem(): mh=" + out.str() + " is out of range");
259  }
260 
261  double r0 = 0.003354;
262  double r1 = -0.000209;
263  double r2 = 0.0000254;
264  double r3 = -0.00000785;
265  double r4 = -0.00000233;
266  double r5 = 0.00783;
267  double r6 = 0.00338;
268  double r7 = -0.00000989;
269  double r8 = 0.0939;
270  double r9 = 0.204;
271  double r10 = -0.103;
272 
273  //double Mw = Mw(DeltaAlphaL5q_i); /* for test */
274  double Mw = Mw_i;
275 
276  double L_H = log(mycache.getSM().getMHl() / 100.0);
277  double Delta_H = mycache.getSM().getMHl() / 100.0;
278  double Delta_t = pow((mycache.getSM().getMtpole() / 178.0), 2.0) - 1.0;
279  double Delta_Z = mycache.getSM().getMz() / 91.1876 - 1.0;
280  double Delta_W = Mw / 80.404 - 1.0;
281 
282  return ( r0 + r1 * L_H + r2 * L_H * L_H + r3 * pow(L_H, 4.0)
283  + r4 * (Delta_H * Delta_H - 1.0) + r5 * Delta_t
284  + r6 * Delta_t * Delta_t + r7 * Delta_t * L_H + r8 * Delta_W
285  + r9 * Delta_W * Delta_t + r10 * Delta_Z);
286 }
287 
288 double EWSMApproximateFormulae::DeltaKappa_l_TwoLoopEW_rem(const double Mw_i) const
289 {
290  // applicable for 10 GeV <= mHl <= 1 TeV
291  if (mycache.getSM().getMHl() < 10.0 || mycache.getSM().getMHl() > UpperBoundForApproximateFormulae) {
292  std::stringstream out;
293  out << mycache.getSM().getMHl();
294  throw std::runtime_error("ApproximateFormulae::DeltaKappa_l_TwoLoopEW_rem(): mh=" + out.str() + " is out of range");
295  }
296 
297  double k0 = -0.002711;
298  double k1 = -0.0000312;
299  double k2 = -0.0000412;
300  double k3 = 0.00000528;
301  double k4 = 0.00000375;
302  double k5 = -0.00516;
303  double k6 = -0.00206;
304  double k7 = -0.000232;
305  double k8 = -0.0647;
306  double k9 = -0.129;
307  double k10 = 0.0712;
308 
309  double L_H = log(mycache.getSM().getMHl() / 100.0);
310  double Delta_H = mycache.getSM().getMHl() / 100.0;
311  double Delta_t = pow((mycache.getSM().getMtpole() / 178.0), 2.0) - 1.0;
312  double Delta_Z = mycache.getSM().getMz() / 91.1876 - 1.0;
313  double Delta_W = Mw_i / 80.404 - 1.0;
314 
315  return ( k0 + k1 * L_H + k2 * L_H * L_H + k3 * pow(L_H, 4.0)
316  + k4 * (Delta_H * Delta_H - 1.0) + k5 * Delta_t
317  + k6 * Delta_t * Delta_t + k7 * Delta_t * L_H
318  + k8 * Delta_W + k9 * Delta_W * Delta_t + k10 * Delta_Z);
319 }
320 
321 double EWSMApproximateFormulae::DeltaKappa_b_TwoLoopEW_rem(const double Mw_i) const
322 {
323  // applicable for 10 GeV <= mHl <= 1 TeV
324  if (mycache.getSM().getMHl() < 10.0 || mycache.getSM().getMHl() > UpperBoundForApproximateFormulae) {
325  std::stringstream out;
326  out << mycache.getSM().getMHl();
327  throw std::runtime_error("ApproximateFormulae::DeltaKappa_b_TwoLoopEW_rem(): mh=" + out.str() + " is out of range");
328  }
329 
330  double k0 = -0.002666;
331  double k1 = -0.0000592;
332  double k2 = -0.00000329;
333  double k3 = 0.00000349;
334  double k4 = 0.00000283;
335  double k5 = -0.00534;
336  double k6 = -0.00210;
337  double k7 = -0.000219;
338  double k8 = -0.0631;
339  double k9 = -0.126;
340  double k10 = 0.0647;
341 
342  double L_H = log(mycache.getSM().getMHl() / 100.0);
343  double Delta_H = mycache.getSM().getMHl() / 100.0;
344  double Delta_t = pow((mycache.getSM().getMtpole() / 178.0), 2.0) - 1.0;
345  double Delta_Z = mycache.getSM().getMz() / 91.1876 - 1.0;
346  double Delta_W = Mw_i / 80.404 - 1.0;
347 
348  return ( k0 + k1 * L_H + k2 * L_H * L_H + k3 * pow(L_H, 4.0)
349  + k4 * (Delta_H * Delta_H - 1.0) + k5 * Delta_t
350  + k6 * Delta_t * Delta_t + k7 * Delta_t * L_H
351  + k8 * Delta_W + k9 * Delta_W * Delta_t + k10 * Delta_Z);
352 }
353 
354 double EWSMApproximateFormulae::R0_bottom_OLD() const
355 {
356  // applicable for 10 GeV <= mHl <= 1 TeV
357  if (mycache.getSM().getMHl() < 10.0 || mycache.getSM().getMHl() > UpperBoundForApproximateFormulae) {
358  std::stringstream out;
359  out << mycache.getSM().getMHl();
360  throw std::runtime_error("ApproximateFormulae::R0_bottom(): mh=" + out.str() + " is out of range");
361  }
362 
363  /*-----------------------------------------------*/
364  /* arXiv:1205.0299v1 by Freitas and Huang */
365  /*
366  double Rb00 = 0.2147464;
367  double c1 = 0.0000221;
368  double c2 = 0.0000026;
369  double c3 = -0.00000067;
370  double c4 = 0.0000000911;
371  double c5 = 0.000647;
372  double c6 = -0.003239;
373  double c7 = 0.0000673;
374  double c8 = -0.000324;
375  double c9 = 0.0610;
376 
377  double L_H = log(mycache.getSM().getMHl()/100.0);
378  double Delta_H = mycache.getSM().getMHl()/100.0;
379  double Delta_ale = mycache.getSM().mycache.getSM().DeltaAlphaL5q()/0.05900 - 1.0;
380  double Delta_t = pow((mycache.getSM().getMtpole()/173.2), 2.0) - 1.0;
381  double Delta_alphas = mycache.getSM().getAlsMz()/0.1184 - 1.0;
382  double Delta_Z = mycache.getSM().getMz()/91.1876 - 1.0;
383 
384  return (Rb00 + c1*L_H + c2*L_H*L_H + c3*pow(L_H, 4.0)
385  + c4*(Delta_H*Delta_H - 1.0) + c5*Delta_ale + c6*Delta_t
386  + c7*Delta_t*L_H + c8*Delta_alphas + c9*Delta_Z );
387  */
388 
389  /*-----------------------------------------------*/
390  /* arXiv:1205.0299v2 by Freitas and Huang */
391  /*
392  double Rb00 = 0.2149246;
393  double c1 = 2.23 * pow(10.0, -5.);
394  double c2 = 2.6 * pow(10.0, -6.);
395  double c3 = -6.8 * pow(10.0, -7.);
396  double c4 = 9.19 *pow(10.0, -8.);
397  double c5 = 6.58 * pow(10.0, -4.);
398  double c6 = -3.363 * pow(10.0, -3.);
399  double c7 = 6.74 * pow(10.0, -5.);
400  double c8 = -1.688 * pow(10.0, -3.);
401  double c9 = -9.26 * pow(10.0, -4.);
402  double c10 = 5.93 * pow(10.0, -2.);
403 
404  double L_H = log(mycache.getSM().getMHl()/100.0);
405  double Delta_H = mycache.getSM().getMHl()/100.0;
406  double Delta_ale = mycache.getSM().mycache.getSM().DeltaAlphaL5q()/0.05900 - 1.0;
407  double Delta_t = pow((mycache.getSM().getMtpole()/173.2), 2.0) - 1.0;
408  double Delta_alphas = mycache.getSM().getAlsMz()/0.1184 - 1.0;
409  double Delta_Z = mycache.getSM().getMz()/91.1876 - 1.0;
410 
411  return (Rb00 + c1*L_H + c2*L_H*L_H + c3*pow(L_H, 4.0)
412  + c4*(Delta_H*Delta_H - 1.0) + c5*Delta_ale + c6*Delta_t
413  + c7*Delta_t*L_H + c8*Delta_alphas + c9*Delta_alphas*Delta_alphas
414  + c10*Delta_Z );
415  */
416 
417  /*-----------------------------------------------*/
418  /* arXiv:1205.0299v3 by Freitas and Huang */
419 
420  double Rb00 = 0.2154940;
421  double c1 = 1.88 * pow(10.0, -5.);
422  double c2 = 2.0 * pow(10.0, -6.);
423  double c3 = -6.0 * pow(10.0, -7.);
424  double c4 = 8.53 * pow(10.0, -8.);
425  double c5 = 7.05 * pow(10.0, -4.);
426  double c6 = -3.159 * pow(10.0, -3.);
427  double c7 = 6.65 * pow(10.0, -5.);
428  double c8 = -1.704 * pow(10.0, -3.);
429  double c9 = -9.30 * pow(10.0, -4.);
430  double c10 = 6.26 * pow(10.0, -2.);
431 
432  double L_H = log(mycache.getSM().getMHl() / 100.0);
433  double Delta_H = mycache.getSM().getMHl() / 100.0;
434  double Delta_ale = mycache.getSM().DeltaAlphaL5q() / 0.05900 - 1.0;
435  double Delta_t = pow((mycache.getSM().getMtpole() / 173.2), 2.0) - 1.0;
436  double Delta_alphas = mycache.getSM().getAlsMz() / 0.1184 - 1.0;
437  double Delta_Z = mycache.getSM().getMz() / 91.1876 - 1.0;
438 
439  /* Debug (parameters in arXiv:1205.0299v3) */
440  //double mHpaper = 100.0;
441  //double mHpaper = 200.0;
442  //double mHpaper = 400.0;
443  //double mHpaper = 600.0;
444  //double mHpaper = 1000.0;
445  //L_H = log(mHpaper/100.0);
446  //Delta_H = mHpaper/100.0;
447  //Delta_ale = 0.0;//0.05900/0.05900 - 1.0;
448  //Delta_t = 0.0;//pow((173.2/173.2), 2.0) - 1.0;
449  //Delta_alphas = 0.0;//0.1184/0.1184 - 1.0;
450  //Delta_Z = 0.0;//91.1876/91.1876 - 1.0;
451 
452  return (Rb00 + c1 * L_H + c2 * L_H * L_H + c3 * pow(L_H, 4.0)
453  + c4 * (Delta_H * Delta_H - 1.0) + c5 * Delta_ale + c6 * Delta_t
454  + c7 * Delta_t * L_H + c8 * Delta_alphas + c9 * Delta_alphas * Delta_alphas
455  + c10 * Delta_Z);
456 
457 }
458 
459 double EWSMApproximateFormulae::Gu_over_Gb_OLD() const
460 {
461  // applicable for 10 GeV <= mHl <= 1 TeV
462  if (mycache.getSM().getMHl() < 10.0 || mycache.getSM().getMHl() > UpperBoundForApproximateFormulae) {
463  std::stringstream out;
464  out << mycache.getSM().getMHl();
465  throw std::runtime_error("ApproximateFormulae::Gu_over_Gb(): mh=" + out.str() + " is out of range");
466  }
467 
468  // obtained from Freitas on Apr. 23, 2013
469  /*
470  double R = 0.8024769;
471  double c1 = -1.9007e-4;
472  double c2 = -2.112e-5;
473  double c3 = 6.63e-6;
474  double c4 = -1.0284e-6;
475  double c5 = -8.081e-3;
476  double c6 = 1.830e-4;
477  double c7 = 1.7522e-2;
478  double c8 = 4.440e-3;
479  double c9 = -3.245e-4;
480  double c10 = 1.8079e-2;
481  double c11 = 1.0720e-2;
482  double c12 = -0.129;
483  */
484 
485  // obtained from Freitas on Sep. 21, 2013
486  double R = 0.7997930;
487  double c1 = -1.7991e-4;
488  double c2 = -1.980e-5;
489  double c3 = 6.24e-6;
490  double c4 = -0.9829e-6;
491  double c5 = -8.200e-3;
492  double c6 = 1.657e-4;
493  double c7 = 1.6476e-2;
494  double c8 = 4.463e-3;
495  double c9 = -3.187e-4;
496  double c10 = 1.8113e-2;
497  double c11 = 1.0720e-2;
498  double c12 = -0.144;
499 
500  double LH = log(mycache.getSM().getMHl() / 100.0);
501  double DH = pow((mycache.getSM().getMHl() / 100.0), 2.0) - 1.0;
502  double Dal = mycache.getSM().DeltaAlphaL5q() / 0.059 - 1.0;
503  double Dt = pow((mycache.getSM().getMtpole() / 173.2), 2.0) - 1.0;
504  double Das = mycache.getSM().getAlsMz() / 0.1184 - 1.0;
505  double Dmz = mycache.getSM().getMz() / 91.1876 - 1.0;
506 
507  return ( R + c1 * LH + c2 * LH * LH + c3 * pow(LH, 4.0) + c4 * DH + c5 * Dal + c6 * LH * Dal
508  + c7 * Dt + c8 * Dt * Dt + c9 * LH * Dt + c10 * Das + c11 * Das * Das + c12 * Dmz);
509 }
510 
511 double EWSMApproximateFormulae::Gd_over_Gb_OLD() const
512 {
513  // applicable for 10 GeV <= mHl <= 1 TeV
514  if (mycache.getSM().getMHl() < 10.0 || mycache.getSM().getMHl() > UpperBoundForApproximateFormulae) {
515  std::stringstream out;
516  out << mycache.getSM().getMHl();
517  throw std::runtime_error("ApproximateFormulae::Gd_over_Gb(): mh=" + out.str() + " is out of range");
518  }
519 
520  // obtained from Freitas on Apr. 23, 2013
521  /*
522  double R = 1.0239191;
523  double c1 = -5.093e-5;
524  double c2 = -7.08e-6;
525  double c3 = 7.4e-7;
526  double c4 = 3.27e-8;
527  double c5 = 8.68e-4;
528  double c6 = 1.064e-4;
529  double c7 = 1.8875e-2;
530  double c8 = 7.093e-3;
531  double c9 = -4.128e-4;
532  double c10 = 1.898e-4;
533  double c11 = -8.0e-6;
534  double c12 = -0.513;
535  */
536 
537  // obtained from Freitas on Sep. 21, 2013
538  double R = 1.0204024;
539  double c1 = -2.242e-5;
540  double c2 = -1.70e-6;
541  double c3 = 2.1e-7;
542  double c4 = 6.38e-8;
543  double c5 = 5.28e-4;
544  double c6 = 0.999e-4;
545  double c7 = 1.7539e-2;
546  double c8 = 7.138e-3;
547  double c9 = -4.041e-4;
548  double c10 = 2.290e-4;
549  double c11 = -8.0e-6;
550  double c12 = -0.530;
551 
552  double LH = log(mycache.getSM().getMHl() / 100.0);
553  double DH = pow((mycache.getSM().getMHl() / 100.0), 2.0) - 1.0;
554  double Dal = mycache.getSM().DeltaAlphaL5q() / 0.059 - 1.0;
555  double Dt = pow((mycache.getSM().getMtpole() / 173.2), 2.0) - 1.0;
556  double Das = mycache.getSM().getAlsMz() / 0.1184 - 1.0;
557  double Dmz = mycache.getSM().getMz() / 91.1876 - 1.0;
558 
559  return ( R + c1 * LH + c2 * LH * LH + c3 * pow(LH, 4.0) + c4 * DH + c5 * Dal + c6 * LH * Dal
560  + c7 * Dt + c8 * Dt * Dt + c9 * LH * Dt + c10 * Das + c11 * Das * Das + c12 * Dmz);
561 }
562 
563 double EWSMApproximateFormulae::X(const std::string observable) const
564 {
565  double LH = log(mycache.getSM().getMHl() / 125.7);
566  double Dt = pow(mycache.getSM().getMtpole() / 173.2, 2.0) - 1.0;
567  double Das = mycache.getSM().getAlsMz() / 0.1184 - 1.0;
568  double Dal = mycache.getSM().DeltaAlphaL5q() / 0.059 - 1.0;
569  double DZ = mycache.getSM().getMz() / 91.1876 - 1.0;
570 
571  double X0, c1, c2, c3, c4, c5, c6, c7;
572  if (observable.compare("Gamma_nu") == 0) {
573  X0 = 167.157;
574  c1 = -0.055;
575  c2 = 1.26;
576  c3 = -0.19;
577  c4 = -0.02;
578  c5 = 0.36;
579  c6 = -0.1;
580  c7 = 503.0;
581  } else if (observable.compare("Gamma_e_mu") == 0) {
582  X0 = 83.966;
583  c1 = -0.047;
584  c2 = 0.807;
585  c3 = -0.095;
586  c4 = -0.01;
587  c5 = 0.25;
588  c6 = -1.1;
589  c7 = 285.0;
590  } else if (observable.compare("Gamma_tau") == 0) {
591  X0 = 83.776;
592  c1 = -0.047;
593  c2 = 0.806;
594  c3 = -0.095;
595  c4 = -0.01;
596  c5 = 0.25;
597  c6 = -1.1;
598  c7 = 285.0;
599  } else if (observable.compare("Gamma_u") == 0) {
600  X0 = 299.936;
601  c1 = -0.34;
602  c2 = 4.07;
603  c3 = 14.27;
604  c4 = 1.6;
605  c5 = 1.8;
606  c6 = -11.1;
607  c7 = 1253.0;
608  } else if (observable.compare("Gamma_c") == 0) {
609  X0 = 299.860;
610  c1 = -0.34;
611  c2 = 4.07;
612  c3 = 14.27;
613  c4 = 1.6;
614  c5 = 1.8;
615  c6 = -11.1;
616  c7 = 1253.0;
617  } else if (observable.compare("Gamma_d_s") == 0) {
618  X0 = 382.770;
619  c1 = -0.34;
620  c2 = 3.83;
621  c3 = 10.20;
622  c4 = -2.4;
623  c5 = 0.67;
624  c6 = -10.1;
625  c7 = 1469.0;
626  } else if (observable.compare("Gamma_b") == 0) {
627  X0 = 375.724;
628  c1 = -0.30;
629  c2 = -2.28;
630  c3 = 10.53;
631  c4 = -2.4;
632  c5 = 1.2;
633  c6 = -10.0;
634  c7 = 1458.0;
635  } else if (observable.compare("GammaZ") == 0) {
636  X0 = 2494.24;
637  c1 = -2.0;
638  c2 = 19.7;
639  c3 = 58.60;
640  c4 = -4.0;
641  c5 = 8.0;
642  c6 = -55.9;
643  c7 = 9267.0;
644  } else if (observable.compare("sigmaHadron") == 0) {
645  X0 = 41488.4;
646  c1 = 3.0;
647  c2 = 60.9;
648  c3 = -579.4;
649  c4 = 38.0;
650  c5 = 7.3;
651  c6 = 85.0;
652  c7 = -86027.0;
653  } else if (observable.compare("R0_lepton") == 0) {
654  X0 = 20750.9;
655  c1 = -8.1;
656  c2 = -39.0;
657  c3 = 732.1;
658  c4 = -44.0;
659  c5 = 5.5;
660  c6 = -358.0;
661  c7 = 11702.0;
662  } else if (observable.compare("R0_charm") == 0) {
663  X0 = 172.23;
664  c1 = -0.029;
665  c2 = 1.0;
666  c3 = 2.3;
667  c4 = 1.3;
668  c5 = 0.38;
669  c6 = -1.2;
670  c7 = 37.0;
671  } else if (observable.compare("R0_bottom") == 0) {
672  X0 = 215.80;
673  c1 = 0.031;
674  c2 = -2.98;
675  c3 = -1.32;
676  c4 = -0.84;
677  c5 = 0.035;
678  c6 = 0.73;
679  c7 = -18.0;
680  } else
681  throw std::runtime_error("ApproximateFormulae::X(): " + observable + " is not defined");
682 
683  return ( 0.001
684  * (X0 + c1 * LH + c2 * Dt + c3 * Das + c4 * Das * Das + c5 * Das * Dt + c6 * Dal + c7 * DZ));
685 }
686 
687 double EWSMApproximateFormulae::X_extended(const std::string observable) const
688 {
689  double LH = log(mycache.getSM().getMHl() / 125.7);
690  double DH = mycache.getSM().getMHl() / 125.7 - 1.0;
691  double Dt = pow(mycache.getSM().getMtpole() / 173.2, 2.0) - 1.0;
692  double Das = mycache.getSM().getAlsMz() / 0.1184 - 1.0;
693  double Dal = mycache.getSM().DeltaAlphaL5q() / 0.059 - 1.0;
694  double DZ = mycache.getSM().getMz() / 91.1876 - 1.0;
695 
696  double X0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15;
697  double ThError = 0.0; // Theoretical uncertainty
698  if (observable.compare("Gamma_nu") == 0) {
699  X0 = 167.157;
700  a1 = -0.1567;
701  a2 = -0.1194;
702  a3 = 0.1031;
703  a4 = -0.00269;
704  a5 = 1.258;
705  a6 = -0.13;
706  a7 = -0.020;
707  a8 = 0.0133;
708  a9 = -0.19;
709  a10 = -0.018;
710  a11 = -0.021;
711  a12 = 0.34;
712  a13 = -0.084;
713  a14 = 0.064;
714  a15 = 503.0;
715  } else if (observable.compare("Gamma_e_mu") == 0) {
716  X0 = 83.966;
717  a1 = -0.1017;
718  a2 = -0.06352;
719  a3 = 0.05500;
720  a4 = -0.00145;
721  a5 = 0.8051;
722  a6 = -0.027;
723  a7 = -0.017;
724  a8 = 0.0066;
725  a9 = -0.095;
726  a10 = -0.010;
727  a11 = -0.015;
728  a12 = 0.23;
729  a13 = -1.1;
730  a14 = 0.064;
731  a15 = 285.0;
732  } else if (observable.compare("Gamma_tau") == 0) {
733  X0 = 83.776;
734  a1 = -0.1016;
735  a2 = -0.06339;
736  a3 = 0.05488;
737  a4 = -0.00145;
738  a5 = 0.8036;
739  a6 = -0.026;
740  a7 = -0.017;
741  a8 = 0.0066;
742  a9 = -0.095;
743  a10 = -0.010;
744  a11 = -0.015;
745  a12 = 0.23;
746  a13 = -1.1;
747  a14 = 0.064;
748  a15 = 285.0;
749  } else if (observable.compare("Gamma_u") == 0) {
750  X0 = 299.936;
751  a1 = -0.5681;
752  a2 = -0.2636;
753  a3 = 0.2334;
754  a4 = -0.00592;
755  a5 = 4.057;
756  a6 = -0.50;
757  a7 = -0.058;
758  a8 = 0.0352;
759  a9 = 14.26;
760  a10 = 1.6;
761  a11 = -0.081;
762  a12 = 1.7;
763  a13 = -11.1;
764  a14 = 0.19;
765  a15 = 1251.0;
766  } else if (observable.compare("Gamma_c") == 0) {
767  X0 = 299.859;
768  a1 = -0.5680;
769  a2 = -0.2635;
770  a3 = 0.2334;
771  a4 = -0.00592;
772  a5 = 4.056;
773  a6 = -0.50;
774  a7 = -0.058;
775  a8 = 0.0352;
776  a9 = 14.26;
777  a10 = 1.6;
778  a11 = -0.081;
779  a12 = 1.7;
780  a13 = -11.1;
781  a14 = 0.19;
782  a15 = 1251.0;
783  } else if (observable.compare("Gamma_d_s") == 0) {
784  X0 = 382.770;
785  a1 = -0.6199;
786  a2 = -0.3182;
787  a3 = 0.2800;
788  a4 = -0.00711;
789  a5 = 3.810;
790  a6 = -0.25;
791  a7 = -0.060;
792  a8 = 0.0420;
793  a9 = 10.20;
794  a10 = -2.4;
795  a11 = -0.083;
796  a12 = 0.65;
797  a13 = -10.1;
798  a14 = 0.19;
799  a15 = 1468.0;
800  } else if (observable.compare("Gamma_b") == 0) {
801  X0 = 375.723;
802  a1 = -0.5744;
803  a2 = -0.3074;
804  a3 = 0.2725;
805  a4 = -0.00703;
806  a5 = -2.292;
807  a6 = -0.027;
808  a7 = -0.013;
809  a8 = 0.0428;
810  a9 = 10.53;
811  a10 = -2.4;
812  a11 = -0.088;
813  a12 = 1.2;
814  a13 = -10.1;
815  a14 = 0.19;
816  a15 = 1456.0;
817  } else if (observable.compare("Gamma_had") == 0) {
818 // Removing leptonic contributions from GammaZ
819 // X0 = 1741.06;
820 // a1 = -2.9501;
821 // a2 = -1.47064;
822 // a3 = 1.29906;
823 // a4 = -0.033105;
824 // a5 = 13.4416;
825 // a6 = -1.5285;
826 // a7 = -0.249;
827 // a8 = 0.19725;
828 // a9 = 59.4525;
829 // a10 = -4.008;
830 // a11 = -0.419;
831 // a12 = 5.895;
832 // a13 = -52.474;
833 // a14 = 0.933;
834 // a15 = 6893.;
835 
836 // Summing all hadronic contributions
837  X0 = 1741.058;
838  a1 = -2.9503;
839  a2 = -1.4708999;
840  a3 = 1.2993;
841  a4 = -0.03308999;
842  a5 = 13.440999;
843  a6 = -1.527;
844  a7 = -0.249;
845  a8 = 0.1972;
846  a9 = 59.44999;
847  a10 = -3.9999;
848  a11 = -0.416;
849  a12 = 5.9;
850  a13 = -52.5;
851  a14 = 0.95;
852  a15 = 6894.;
853 
854  } else if (observable.compare("GammaZ") == 0) {
855  X0 = 2494.24;
856  a1 = -3.725;
857  a2 = -2.019;
858  a3 = 1.773;
859  a4 = -0.04554;
860  a5 = 19.63;
861  a6 = -2.0;
862  a7 = -0.36;
863  a8 = 0.257;
864  a9 = 58.60;
865  a10 = -4.1;
866  a11 = -0.53;
867  a12 = 7.6;
868  a13 = -56.0;
869  a14 = 1.3;
870  a15 = 9256.0;
871  ThError = mycache.getSM().getDelGammaZ();
872  } else if (observable.compare("sigmaHadron") == 0) {
873  X0 = 41488.4;
874  a1 = 3.88;
875  a2 = 0.829;
876  a3 = -0.911;
877  a4 = 0.0076;
878  a5 = 61.10;
879  a6 = 16.0;
880  a7 = -2.0;
881  a8 = -0.59;
882  a9 = -579.4;
883  a10 = 38.0;
884  a11 = -0.26;
885  a12 = 6.5;
886  a13 = 84.0;
887  a14 = 9.5;
888  a15 = -86152.0;
889  ThError = mycache.getSM().getDelSigma0H();
890  } else if (observable.compare("R0_lepton") == 0) {
891  X0 = 20750.9;
892  a1 = -10.00;
893  a2 = -1.83;
894  a3 = 1.878;
895  a4 = -0.0343;
896  a5 = -38.8;
897  a6 = -11.0;
898  a7 = 1.2;
899  a8 = 0.72;
900  a9 = 732.1;
901  a10 = -44.0;
902  a11 = -0.64;
903  a12 = 5.6;
904  a13 = -357.0;
905  a14 = -4.7;
906  a15 = 11771.0;
907  ThError = mycache.getSM().getDelR0l();
908  } else if (observable.compare("R0_electron") == 0) {
909  ThError = mycache.getSM().getDelR0l();
910  return X_extended("Gamma_had")/X_extended("Gamma_e_mu") + ThError;
911 
912  } else if (observable.compare("R0_muon") == 0) {
913  ThError = mycache.getSM().getDelR0l();
914  return X_extended("Gamma_had")/X_extended("Gamma_e_mu") + ThError;
915 
916  } else if (observable.compare("R0_tau") == 0) {
917  ThError = mycache.getSM().getDelR0l();
918  return X_extended("Gamma_had")/X_extended("Gamma_tau") + ThError;
919 
920  } else if (observable.compare("R0_up") == 0) {
921  ThError = 0.0; // Set to zero for the moment
922  return X_extended("Gamma_u")/X_extended("Gamma_had") + ThError;
923 
924  } else if (observable.compare("R0_strange") == 0) {
925  ThError = 0.0; // Set to zero for the moment
926  return X_extended("Gamma_d_s")/X_extended("Gamma_had") + ThError;
927 
928  } else if (observable.compare("R0_charm") == 0) {
929  X0 = 172.23;
930  a1 = -0.034;
931  a2 = -0.0058;
932  a3 = 0.0054;
933  a4 = -0.00012;
934  a5 = 1.00;
935  a6 = -0.15;
936  a7 = -0.0074;
937  a8 = 0.00091;
938  a9 = 2.3;
939  a10 = 1.3;
940  a11 = -0.0013;
941  a12 = 0.35;
942  a13 = -1.2;
943  a14 = 0.014;
944  a15 = 37.0;
945  ThError = mycache.getSM().getDelR0c();
946  } else if (observable.compare("R0_bottom") == 0) {
947  X0 = 215.80;
948  a1 = 0.036;
949  a2 = 0.0057;
950  a3 = -0.0044;
951  a4 = 0.000062;
952  a5 = -2.98;
953  a6 = 0.20;
954  a7 = 0.020;
955  a8 = -0.00036;
956  a9 = -1.3;
957  a10 = -0.84;
958  a11 = -0.0019;
959  a12 = 0.054;
960  a13 = 0.73;
961  a14 = -0.011;
962  a15 = -18.0;
963  ThError = mycache.getSM().getDelR0b();
964  } else
965  throw std::runtime_error("ApproximateFormulae::X_extended(): " + observable + " is not defined");
966 
967  return ( 0.001
968  * (X0 + a1 * LH + a2 * LH * LH + a3 * DH + a4 * DH * DH + a5 * Dt + a6 * Dt * Dt
969  + a7 * Dt * LH + a8 * Dt * LH * LH + a9 * Das + a10 * Das * Das + a11 * Das * LH
970  + a12 * Das * Dt + a13 * Dal + a14 * Dal * LH + a15 * DZ) + ThError);
971 }
972 
973 
974 double EWSMApproximateFormulae::X_full_2_loop(const std::string observable) const
975 {
976 
977 // For MH not in [85,165] GeV there are significant differences with some predicions
978 // of X_extended, which go well beyond the expected size of the bosonic corrections (>~2x).
979 // Use EWSMApproximateFormulae::X_extended in that case
980  if (mycache.getSM().getMHl() < 85.0 || mycache.getSM().getMHl() > 165.0) {
981  return X_extended(observable);
982  }
983 
984 // Otherwise proceed with the full 2-loop code
985  double LH = log(mycache.getSM().getMHl() / 125.7);
986  double Dt = pow(mycache.getSM().getMtpole() / 173.2, 2.0) - 1.0;
987  double Das = mycache.getSM().getAlsMz() / 0.1184 - 1.0;
988  double Dal = mycache.getSM().DeltaAlphaL5q() / 0.059 - 1.0;
989  double DZ = mycache.getSM().getMz() / 91.1876 - 1.0;
990 
991  double X0, c1, c2, c3, c4, c5, c6, c7;
992  double ThError = 0.0; // Theoretical uncertainty
993  if (observable.compare("Gamma_nu") == 0) {
994  X0 = 167.176;
995  c1 = -0.071;
996  c2 = 1.26;
997  c3 = -0.19;
998  c4 = -0.02;
999  c5 = 0.36;
1000  c6 = -0.1;
1001  c7 = 504.0;
1002 
1003  } else if (observable.compare("Gamma_e_mu") == 0) {
1004  X0 = 83.983;
1005  c1 = -0.061;
1006  c2 = 0.810;
1007  c3 = -0.096;
1008  c4 = -0.01;
1009  c5 = 0.25;
1010  c6 = -1.1;
1011  c7 = 286.0;
1012 
1013  } else if (observable.compare("Gamma_tau") == 0) {
1014  X0 = 83.793;
1015  c1 = -0.060;
1016  c2 = 0.810;
1017  c3 = -0.095;
1018  c4 = -0.01;
1019  c5 = 0.25;
1020  c6 = -1.1;
1021  c7 = 285.0;
1022 
1023  } else if (observable.compare("Gamma_u") == 0) {
1024  X0 = 299.993;
1025  c1 = -0.38;
1026  c2 = 4.08;
1027  c3 = 14.27;
1028  c4 = 1.6;
1029  c5 = 1.8;
1030  c6 = -11.1;
1031  c7 = 1253.0;
1032 
1033  } else if (observable.compare("Gamma_c") == 0) {
1034  X0 = 299.916;
1035  c1 = -0.38;
1036  c2 = 4.08;
1037  c3 = 14.27;
1038  c4 = 1.6;
1039  c5 = 1.8;
1040  c6 = -11.1;
1041  c7 = 1253.0;
1042 
1043  } else if (observable.compare("Gamma_d_s") == 0) {
1044  X0 = 382.828;
1045  c1 = -0.39;
1046  c2 = 3.83;
1047  c3 = 10.20;
1048  c4 = -2.4;
1049  c5 = 0.67;
1050  c6 = -10.1;
1051  c7 = 1470.0;
1052 
1053  } else if (observable.compare("Gamma_b") == 0) {
1054  X0 = 375.889;
1055  c1 = -0.36;
1056  c2 = -2.14;
1057  c3 = 10.53;
1058  c4 = -2.4;
1059  c5 = 1.2;
1060  c6 = -10.1;
1061  c7 = 1459.0;
1062 
1063  } else if (observable.compare("Gamma_had") == 0) {
1064 
1065 // Summing all hadronic contributions
1066  X0 = 1741.454;
1067  c1 = -1.9;
1068  c2 = 13.68;
1069  c3 = 59.47;
1070  c4 = -4.0;
1071  c5 = 6.14;
1072  c6 = -52.5;
1073  c7 = 6905.0;
1074 
1075  } else if (observable.compare("GammaZ") == 0) {
1076  X0 = 2494.74;
1077  c1 = -2.3;
1078  c2 = 19.9;
1079  c3 = 58.61;
1080  c4 = -4.0;
1081  c5 = 8.0;
1082  c6 = -56.0;
1083  c7 = 9273.0;
1084 
1085  ThError = mycache.getSM().getDelGammaZ();
1086  } else if (observable.compare("sigmaHadron") == 0) {
1087  X0 = 41489.6;
1088  c1 = 1.6;
1089  c2 = 60.0;
1090  c3 = -579.6;
1091  c4 = 38.0;
1092  c5 = 7.3;
1093  c6 = 85.0;
1094  c7 = -86011.0;
1095 
1096  ThError = mycache.getSM().getDelSigma0H();
1097  } else if (observable.compare("R0_lepton") == 0) {
1098  X0 = 20751.6;
1099  c1 = -7.8;
1100  c2 = -37.0;
1101  c3 = 732.3;
1102  c4 = -44.0;
1103  c5 = 5.5;
1104  c6 = -358.0;
1105  c7 = 11696.0;
1106 
1107  ThError = mycache.getSM().getDelR0l();
1108  } else if (observable.compare("R0_electron") == 0) {
1109  ThError = mycache.getSM().getDelR0l();
1110  return X_full_2_loop("Gamma_had")/X_full_2_loop("Gamma_e_mu") + ThError;
1111 
1112  } else if (observable.compare("R0_muon") == 0) {
1113  ThError = mycache.getSM().getDelR0l();
1114  return X_full_2_loop("Gamma_had")/X_full_2_loop("Gamma_e_mu") + ThError;
1115 
1116  } else if (observable.compare("R0_tau") == 0) {
1117  ThError = mycache.getSM().getDelR0l();
1118  return X_full_2_loop("Gamma_had")/X_full_2_loop("Gamma_tau") + ThError;
1119 
1120  } else if (observable.compare("R0_neutrino") == 0) {
1121  ThError = 0.0;
1122  return X_full_2_loop("Gamma_nu")/X_full_2_loop("Gamma_had") + ThError;
1123 
1124  } else if (observable.compare("R0_up") == 0) {
1125  ThError = 0.0; // Set to zero for the moment
1126  return X_full_2_loop("Gamma_u")/X_full_2_loop("Gamma_had") + ThError;
1127 
1128  } else if (observable.compare("R0_strange") == 0) {
1129  ThError = 0.0; // Set to zero for the moment
1130  return X_full_2_loop("Gamma_d_s")/X_full_2_loop("Gamma_had") + ThError;
1131 
1132  } else if (observable.compare("R0_charm") == 0) {
1133  X0 = 172.22;
1134  c1 = -0.031;
1135  c2 = 1.0;
1136  c3 = 2.3;
1137  c4 = 1.3;
1138  c5 = 0.38;
1139  c6 = -1.2;
1140  c7 = 37.0;
1141 
1142  ThError = mycache.getSM().getDelR0c();
1143  } else if (observable.compare("R0_bottom") == 0) {
1144  X0 = 215.85;
1145  c1 = 0.029;
1146  c2 = -2.92;
1147  c3 = -1.32;
1148  c4 = -0.84;
1149  c5 = 0.032;
1150  c6 = 0.72;
1151  c7 = -18.0;
1152 
1153  ThError = mycache.getSM().getDelR0b();
1154  } else
1155  throw std::runtime_error("ApproximateFormulae::X_full_2_loop(): " + observable + " is not defined");
1156 
1157  return ( 0.001
1158  * (X0 + c1 * LH + c2 * Dt
1159  + c3 * Das + c4 * Das * Das
1160  + c5 * Das * Dt + c6 * Dal + c7 * DZ) + ThError);
1161 }
1162 
1163 
1164 
1165 double EWSMApproximateFormulae::X_full(const std::string observable) const
1166 {
1167 
1168 // Full 2-loop implementation
1169 
1170  double LH = log(mycache.getSM().getMHl() / 125.7);
1171  double LH2 = LH * LH;
1172 
1173  double DH = pow(125.7 / (mycache.getSM().getMHl()), 4.0) - 1.0;
1174 
1175  double Dt = pow(mycache.getSM().getMtpole() / 173.2, 2.0) - 1.0;
1176 
1177  double Das = mycache.getSM().getAlsMz() / 0.1184 - 1.0;
1178 
1179  double Dal = mycache.getSM().DeltaAlphaL5q() / 0.059 - 1.0;
1180 
1181  double DZ = mycache.getSM().getMz() / 91.1876 - 1.0;
1182 
1183  double X0, a1, a2, a3, a4, a5, a6, a7, a8, a9, a10, a11, a12, a13, a14, a15, a16;
1184 
1185  double ThError = 0.0; // Theoretical uncertainty
1186 
1187  if (observable.compare("Gamma_e_mu") == 0) {
1188  X0 = 83.983;
1189  a1 = -0.1202;
1190  a2 = -0.06919;
1191  a3 = 0.00383;
1192  a4 = 0.0597;
1193  a5 = 0.8037;
1194  a6 = -0.015;
1195  a7 = -0.0195;
1196  a8 = 0.0032;
1197  a9 = -0.0956;
1198  a10 = -0.0078;
1199  a11 = -0.0095;
1200  a12 = 0.25;
1201  a13 = -1.08;
1202  a14 = 0.056;
1203  a15 = -0.37;
1204  a16 = 286.0;
1205 
1206  } else if (observable.compare("Gamma_tau") == 0) {
1207  X0 = 83.793;
1208  a1 = -0.1200;
1209  a2 = -0.06905;
1210  a3 = 0.00382;
1211  a4 = 0.0596;
1212  a5 = 0.8023;
1213  a6 = -0.015;
1214  a7 = -0.0195;
1215  a8 = 0.0032;
1216  a9 = -0.0954;
1217  a10 = -0.0078;
1218  a11 = -0.0094;
1219  a12 = 0.25;
1220  a13 = -1.08;
1221  a14 = 0.056;
1222  a15 = -0.37;
1223  a16 = 285.0;
1224 
1225  } else if (observable.compare("Gamma_nu") == 0) {
1226  X0 = 167.176;
1227  a1 = -0.1752;
1228  a2 = -0.1249;
1229  a3 = 0.00595;
1230  a4 = 0.1046;
1231  a5 = 1.253;
1232  a6 = -0.110;
1233  a7 = -0.0232;
1234  a8 = 0.0064;
1235  a9 = -0.187;
1236  a10 = -0.014;
1237  a11 = -0.014;
1238  a12 = 0.37;
1239  a13 = -0.085;
1240  a14 = 0.054;
1241  a15 = -0.30;
1242  a16 = 503.0;
1243 
1244  } else if (observable.compare("Gamma_u") == 0) {
1245  X0 = 299.994;
1246  a1 = -0.6152;
1247  a2 = -0.2771;
1248  a3 = 0.0174;
1249  a4 = 0.2341;
1250  a5 = 4.051;
1251  a6 = -0.467;
1252  a7 = -0.0676;
1253  a8 = 0.017;
1254  a9 = 14.26;
1255  a10 = 1.6;
1256  a11 = -0.046;
1257  a12 = 1.82;
1258  a13 = -11.1;
1259  a14 = 0.16;
1260  a15 = -1.0;
1261  a16 = 1253.0;
1262 
1263  } else if (observable.compare("Gamma_c") == 0) {
1264  X0 = 299.918;
1265  a1 = -0.6152;
1266  a2 = -0.2771;
1267  a3 = 0.0174;
1268  a4 = 0.2340;
1269  a5 = 4.051;
1270  a6 = -0.467;
1271  a7 = -0.0676;
1272  a8 = 0.017;
1273  a9 = 14.26;
1274  a10 = 1.6;
1275  a11 = -0.046;
1276  a12 = 1.82;
1277  a13 = -11.1;
1278  a14 = 0.16;
1279  a15 = -1.0;
1280  a16 = 1252.0;
1281 
1282  } else if (observable.compare("Gamma_d_s") == 0) {
1283  X0 = 382.829;
1284  a1 = -0.6685;
1285  a2 = -0.3322;
1286  a3 = 0.0193;
1287  a4 = 0.2792;
1288  a5 = 3.792;
1289  a6 = -0.18;
1290  a7 = -0.0706;
1291  a8 = 0.020;
1292  a9 = 10.20;
1293  a10 = -2.4;
1294  a11 = -0.052;
1295  a12 = 0.71;
1296  a13 = -10.1;
1297  a14 = 0.16;
1298  a15 = -0.92;
1299  a16 = 1469.0;
1300 
1301  } else if (observable.compare("Gamma_b") == 0) {
1302  X0 = 375.890;
1303  a1 = -0.6017;
1304  a2 = -0.3158;
1305  a3 = 0.0190;
1306  a4 = 0.227;
1307  a5 = -2.174;
1308  a6 = 0.042;
1309  a7 = -0.027;
1310  a8 = 0.021;
1311  a9 = 10.53;
1312  a10 = -2.4;
1313  a11 = -0.056;
1314  a12 = 1.2;
1315  a13 = -10.1;
1316  a14 = 0.15;
1317  a15 = -0.95;
1318  a16 = 1458.0;
1319 
1320  } else if (observable.compare("Gamma_had") == 0) {
1321 
1322 // Removing leptonic contributions from the total Z witdh
1323  X0 = 1741.46;
1324  a1 = -3.169;
1325  a2 = -1.53487;
1326  a3 = 0.09267;
1327  a4 = 1.2532;
1328  a5 = 13.5113;
1329  a6 = -1.255;
1330  a7 = -0.3039;
1331  a8 = 0.0912;
1332  a9 = 59.4576;
1333  a10 = -3.9346;
1334  a11 = -0.2496;
1335  a12 = 6.24;
1336  a13 = -52.605;
1337  a14 = 0.77;
1338  a15 = -4.79;
1339  a16 = 6901.0;
1340 
1341  } else if (observable.compare("GammaZ") == 0) {
1342  X0 = 2494.75;
1343  a1 = -4.055;
1344  a2 = -2.117;
1345  a3 = 0.122;
1346  a4 = 1.746;
1347  a5 = 19.68;
1348  a6 = -1.63;
1349  a7 = -0.432;
1350  a8 = 0.12;
1351  a9 = 58.61;
1352  a10 = -4.0;
1353  a11 = -0.32;
1354  a12 = 8.1;
1355  a13 = -56.1;
1356  a14 = 1.1;
1357  a15 = -6.8;
1358  a16 = 9267.0;
1359 
1360  ThError = mycache.getSM().getDelGammaZ();
1361  } else if (observable.compare("sigmaHadron") == 0) {
1362  X0 = 41489.6;
1363  a1 = 0.408;
1364  a2 = -0.320;
1365  a3 = 0.0424;
1366  a4 = 1.32;
1367  a5 = 60.17;
1368  a6 = 16.3;
1369  a7 = -2.31;
1370  a8 = -0.19;
1371  a9 = -579.58;
1372  a10 = 38.0;
1373  a11 = 0.010;
1374  a12 = 7.5;
1375  a13 = 85.2;
1376  a14 = 9.1;
1377  a15 = -68.0;
1378  a16 = -85957.0;
1379 
1380  ThError = mycache.getSM().getDelSigma0H();
1381  } else if (observable.compare("R0_lepton") == 0) {
1382  X0 = 20751.6;
1383  a1 = -8.112;
1384  a2 = -1.174;
1385  a3 = 0.155;
1386  a4 = 0.16;
1387  a5 = -37.59;
1388  a6 = -10.9;
1389  a7 = 1.27;
1390  a8 = 0.29;
1391  a9 = 732.30;
1392  a10 = -44.0;
1393  a11 = -0.61;
1394  a12 = 5.7;
1395  a13 = -358.0;
1396  a14 = -4.7;
1397  a15 = 37;
1398  a16 = 11649.0;
1399 
1400  ThError = mycache.getSM().getDelR0l();
1401  } else if (observable.compare("R0_electron") == 0) {
1402  ThError = mycache.getSM().getDelR0l();
1403  return X_full("Gamma_had")/X_full("Gamma_e_mu") + ThError;
1404 
1405  } else if (observable.compare("R0_muon") == 0) {
1406  ThError = mycache.getSM().getDelR0l();
1407  return X_full("Gamma_had")/X_full("Gamma_e_mu") + ThError;
1408 
1409  } else if (observable.compare("R0_tau") == 0) {
1410  ThError = mycache.getSM().getDelR0l();
1411  return X_full("Gamma_had")/X_full("Gamma_tau") + ThError;
1412 
1413  } else if (observable.compare("R0_neutrino") == 0) {
1414  ThError = 0.0;
1415  return X_full("Gamma_nu")/X_full("Gamma_had") + ThError;
1416 
1417  } else if (observable.compare("R0_up") == 0) {
1418  ThError = 0.0; // Set to zero for the moment
1419  return X_full("Gamma_u")/X_full("Gamma_had") + ThError;
1420 
1421  } else if (observable.compare("R0_strange") == 0) {
1422  ThError = 0.0; // Set to zero for the moment
1423  return X_full("Gamma_d_s")/X_full("Gamma_had") + ThError;
1424 
1425  } else if (observable.compare("R0_charm") == 0) {
1426  X0 = 172.222;
1427  a1 = -0.04049;
1428  a2 = -0.00749;
1429  a3 = 0.000832;
1430  a4 = 0.0108;
1431  a5 = 0.98956;
1432  a6 = -0.151;
1433  a7 = -0.00761;
1434  a8 = 0.00080;
1435  a9 = 2.309;
1436  a10 = 1.25;
1437  a11 = 0.00045;
1438  a12 = 0.369;
1439  a13 = -1.20;
1440  a14 = 0.012;
1441  a15 = -0.062;
1442  a16 = 36.67;
1443 
1444  ThError = mycache.getSM().getDelR0c();
1445  } else if (observable.compare("R0_bottom") == 0) {
1446  X0 = 215.850;
1447  a1 = 0.04904;
1448  a2 = 0.009149;
1449  a3 = -0.000535;
1450  a4 = -0.02676;
1451  a5 = -2.9221;
1452  a6 = 0.200;
1453  a7 = 0.0197;
1454  a8 = -0.0011;
1455  a9 = -1.319;
1456  a10 = -0.84;
1457  a11 = -0.0027;
1458  a12 = 0.044;
1459  a13 = 0.719;
1460  a14 = -0.0077;
1461  a15 = -0.044;
1462  a16 = -17.90;
1463 
1464  ThError = mycache.getSM().getDelR0b();
1465  } else
1466  throw std::runtime_error("ApproximateFormulae::X_full(): " + observable + " is not defined");
1467 
1468  return ( 0.001
1469  * ( X0 + a1 * LH + a2 * LH2 + a3 * LH2 * LH2 + a4 *DH
1470  + a5 * Dt + a6 * Dt*Dt + a7 * Dt*LH + a8 * Dt * LH2
1471  + a9 * Das + a10 * Das * Das + a11 * Das*DH + a12*Das*Dt
1472  + a13 * Dal + a14 * Dal * DH + a15 * Dal * Dt
1473  + a16 * DZ ) + ThError);
1474 }
1475 
1476 
1477 double EWSMApproximateFormulae::sin2thetaEff_b_full() const
1478 {
1479  // applicable for 25 GeV <= mHl <= 225 GeV. Remove boundaries for the moment
1480  //if (mycache.getSM().getMHl() < 25.0 || mycache.getSM().getMHl() > 225.0) {
1481  // std::stringstream out;
1482  // out << mycache.getSM().getMHl();
1483  // throw std::runtime_error("ApproximateFormulae::sin2thetaEff_b_full(): mh=" + out.str() + " is out of range");
1484  //}
1485 
1486 // Full 2-loop implementation
1487 
1488  double LH = log(mycache.getSM().getMHl() / 125.7);
1489  double LH2 = LH * LH;
1490 
1491 // double DH = pow(125.7 / (mycache.getSM().getMHl()), 4.0) - 1.0;
1492 
1493  double Dt = pow(mycache.getSM().getMtpole() / 173.2, 2.0) - 1.0;
1494 
1495  double Das = mycache.getSM().getAlsMz() / 0.1184 - 1.0;
1496 
1497  double Dal = mycache.getSM().DeltaAlphaL5q() / 0.059 - 1.0;
1498 
1499  double DZ = mycache.getSM().getMz() / 91.1876 - 1.0;
1500 
1501  double X0, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10;
1502 
1503  double ThError = 0.0; // Theoretical uncertainty
1504 
1505  X0 = 2327.04;
1506  d1 = 4.638;
1507  d2 = 0.558;
1508  d3 = -0.0700;
1509  d4 = 207.0;
1510  d5 = -9.554;
1511  d6 = 3.83;
1512  d7 = 0.179;
1513  d8 = 2.41;
1514  d9 = -8.24;
1515  d10 = -6630.0;
1516 
1517  ThError = mycache.getSM().getDelSin2th_b();
1518 
1519  return ( 0.0001*( X0 + d1 * LH + d2 * LH2 + d3 * LH2 * LH2
1520  + d4 * Dal
1521  + d5 * Dt + d6 * Dt * Dt
1522  + d7 * Dt * LH
1523  + d8 * Das + d9 * Das * Dt
1524  + d10 * DZ ) + ThError);
1525 }
1526 
1527 
1528 double EWSMApproximateFormulae::sin2thetaEff_l_full() const
1529 {
1530  // applicable for 25 GeV <= mHl <= 225 GeV. Remove boundaries for the moment
1531  //if (mycache.getSM().getMHl() < 25.0 || mycache.getSM().getMHl() > 225.0) {
1532  // std::stringstream out;
1533  // out << mycache.getSM().getMHl();
1534  // throw std::runtime_error("ApproximateFormulae::sin2thetaEff_l_full(): mh=" + out.str() + " is out of range");
1535  //}
1536 
1537 // Full 2-loop implementation
1538 
1539  double LH = log(mycache.getSM().getMHl() / 125.7);
1540  double LH2 = LH * LH;
1541 
1542 // double DH = pow(125.7 / (mycache.getSM().getMHl()), 4.0) - 1.0;
1543 
1544  double Dt = pow(mycache.getSM().getMtpole() / 173.2, 2.0) - 1.0;
1545 
1546  double Das = mycache.getSM().getAlsMz() / 0.1184 - 1.0;
1547 
1548  double Dal = mycache.getSM().DeltaAlphaL5q() / 0.059 - 1.0;
1549 
1550  double DZ = mycache.getSM().getMz() / 91.1876 - 1.0;
1551 
1552  double X0, d1, d2, d3, d4, d5, d6, d7, d8, d9, d10;
1553 
1554  double ThError = 0.0; // Theoretical uncertainty
1555 
1556  X0 = 2314.64;
1557  d1 = 4.616;
1558  d2 = 0.539;
1559  d3 = -0.0737;
1560  d4 = 206.0;
1561  d5 = -25.71;
1562  d6 = 4.00;
1563  d7 = 0.288;
1564  d8 = 3.88;
1565  d9 = -6.49;
1566  d10 = -6560.0;
1567 
1568  ThError = mycache.getSM().getDelSin2th_l();
1569 
1570  return ( 0.0001*( X0 + d1 * LH + d2 * LH2 + d3 * LH2 * LH2
1571  + d4 * Dal
1572  + d5 * Dt + d6 * Dt * Dt
1573  + d7 * Dt * LH
1574  + d8 * Das + d9 * Das * Dt
1575  + d10 * DZ ) + ThError);
1576 }
QCD::TAU
Definition: QCD.h:316
QCD::NEUTRINO_3
Definition: QCD.h:315
EWSMApproximateFormulae::X_full
double X_full(const std::string observable) const
, , , , , , , , , , , or .
Definition: EWSMApproximateFormulae.cpp:1165
QCD::BOTTOM
Definition: QCD.h:329
EWSMApproximateFormulae::sin2thetaEff_q
double sin2thetaEff_q(const QCD::quark q) const
with the full two-loop EW corrections (bosonic two-loop EW corrections are missing for ).
Definition: EWSMApproximateFormulae.cpp:133
StandardModel::getDelSigma0H
double getDelSigma0H() const
A get method to retrieve the theoretical uncertainty in , denoted as .
Definition: StandardModel.h:828
StandardModel::getAlsMz
double getAlsMz() const
A get method to access the value of .
Definition: StandardModel.h:727
StandardModel::getDelR0l
double getDelR0l() const
A get method to retrieve the theoretical uncertainty in , denoted as .
Definition: StandardModel.h:838
QCD::UP
Definition: QCD.h:324
EWSMApproximateFormulae::DeltaKappa_b_TwoLoopEW_rem
double DeltaKappa_b_TwoLoopEW_rem(const double Mw_i) const
.
Definition: EWSMApproximateFormulae.cpp:321
EWSMApproximateFormulae::DeltaKappa_l_TwoLoopEW_rem
double DeltaKappa_l_TwoLoopEW_rem(const double Mw_i) const
.
Definition: EWSMApproximateFormulae.cpp:288
QCD::CHARM
Definition: QCD.h:326
StandardModel::DeltaAlphaL5q
double DeltaAlphaL5q() const
The sum of the leptonic and the five-flavour hadronic corrections to the electromagnetic coupling at...
Definition: StandardModel.cpp:830
EWSMApproximateFormulae::Gu_over_Gb_OLD
double Gu_over_Gb_OLD() const
.
Definition: EWSMApproximateFormulae.cpp:459
QCD::NEUTRINO_2
Definition: QCD.h:313
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
QCD::ELECTRON
Definition: QCD.h:312
StandardModel::getDelMw
double getDelMw() const
A get method to retrieve the theoretical uncertainty in , denoted as .
Definition: StandardModel.h:775
EWSMApproximateFormulae::EWSMApproximateFormulae
EWSMApproximateFormulae(const EWSMcache &cache_i)
Constructor.
Definition: EWSMApproximateFormulae.cpp:15
EWSMApproximateFormulae::sin2thetaEff_b
double sin2thetaEff_b() const
with the full two-loop EW corrections.
Definition: EWSMApproximateFormulae.cpp:216
EWSMApproximateFormulae::mycache
const EWSMcache & mycache
A reference to an object of type StandardModel.
Definition: EWSMApproximateFormulae.h:429
EWSMApproximateFormulae::X_extended
double X_extended(const std::string observable) const
, , , , , , , , , , , or .
Definition: EWSMApproximateFormulae.cpp:687
EWSMcache::getSM
const StandardModel & getSM() const
Definition: EWSMcache.h:56
EWSMApproximateFormulae::R0_bottom_OLD
double R0_bottom_OLD() const
.
Definition: EWSMApproximateFormulae.cpp:354
EWSMApproximateFormulae::X
double X(const std::string observable) const
, , , , , , , , , , , or .
Definition: EWSMApproximateFormulae.cpp:563
EWSMApproximateFormulae::Mw
double Mw() const
The -boson mass with the full two-loop EW corrections.
Definition: EWSMApproximateFormulae.cpp:23
EWSMApproximateFormulae::sin2thetaEff_b_full
double sin2thetaEff_b_full() const
with the full two-loop EW corrections.
Definition: EWSMApproximateFormulae.cpp:1477
QCD::TOP
Definition: QCD.h:328
gslpp::pow
complex pow(const complex &z1, const complex &z2)
Definition: gslpp_complex.cpp:395
EWSMApproximateFormulae::X_full_2_loop
double X_full_2_loop(const std::string observable) const
, , , , , , , , , , , or .
Definition: EWSMApproximateFormulae.cpp:974
EWSMApproximateFormulae.h
EWSMcache
A class for cache variables used in computing radiative corrections to the EW precision observables.
Definition: EWSMcache.h:40
QCD::getMtpole
double getMtpole() const
A get method to access the pole mass of the top quark.
Definition: QCD.h:588
StandardModel::getDelSin2th_l
double getDelSin2th_l() const
A get method to retrieve the theoretical uncertainty in , denoted as .
Definition: StandardModel.h:786
EWSMApproximateFormulae::sin2thetaEff_l
double sin2thetaEff_l(const QCD::lepton l) const
with the full two-loop EW corrections.
Definition: EWSMApproximateFormulae.cpp:74
EWSMApproximateFormulae::DeltaR_TwoLoopEW_rem
double DeltaR_TwoLoopEW_rem(const double Mw_i) const
.
Definition: EWSMApproximateFormulae.cpp:252
QCD::quark
quark
An enum type for quarks.
Definition: QCD.h:323
StandardModel::getDelSin2th_q
double getDelSin2th_q() const
A get method to retrieve the theoretical uncertainty in , denoted as .
Definition: StandardModel.h:797
StandardModel::getMz
double getMz() const
A get method to access the mass of the boson .
Definition: StandardModel.h:718
EWSMApproximateFormulae::sin2thetaEff_l_full
double sin2thetaEff_l_full() const
with the full two-loop EW corrections.
Definition: EWSMApproximateFormulae.cpp:1528
StandardModel::getDelR0b
double getDelR0b() const
A get method to retrieve the theoretical uncertainty in , denoted as .
Definition: StandardModel.h:858
Mw
An observable class for the -boson mass.
Definition: Mw.h:22
StandardModel::getDelR0c
double getDelR0c() const
A get method to retrieve the theoretical uncertainty in , denoted as .
Definition: StandardModel.h:848
QCD::STRANGE
Definition: QCD.h:327
StandardModel::getMHl
virtual double getMHl() const
A get method to retrieve the Higgs mass .
Definition: StandardModel.h:765
StandardModel::getDelSin2th_b
double getDelSin2th_b() const
A get method to retrieve the theoretical uncertainty in , denoted as .
Definition: StandardModel.h:808
StandardModel::getDelGammaZ
double getDelGammaZ() const
A get method to retrieve the theoretical uncertainty in , denoted as .
Definition: StandardModel.h:818
QCD::DOWN
Definition: QCD.h:325
QCD::NEUTRINO_1
Definition: QCD.h:311
EWSMApproximateFormulae::Gd_over_Gb_OLD
double Gd_over_Gb_OLD() const
.
Definition: EWSMApproximateFormulae.cpp:511
QCD::MU
Definition: QCD.h:314
QCD::lepton
lepton
An enum type for leptons.
Definition: QCD.h:310