a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
EWSMTwoLoopQCD.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 <stdexcept>
9 #include <gsl/gsl_sf.h>
10 #include "EWSMTwoLoopQCD.h"
11 
12 EWSMTwoLoopQCD::EWSMTwoLoopQCD(const EWSMcache& cache_i)
13 : cache(cache_i)
14 {
15 }
16 
17 
19 
20 double EWSMTwoLoopQCD::DeltaAlpha_l(const double s) const
21 {
22  return (0.0);
23 }
24 
25 double EWSMTwoLoopQCD::DeltaAlpha_t(const double s) const
26 {
27  double xt = s / cache.getSM().getMtpole() / cache.getSM().getMtpole();
28  double als;
29  if (s == cache.getSM().getMz() * cache.getSM().getMz())
30  als = cache.getSM().getAlsMz();
31  else
32  als = cache.Als(sqrt(s), FULLNNLO);
33  double tmp = (5.062 + xt * 0.8315) * als / M_PI;
34  tmp *= -4.0 / 45.0 * cache.getSM().getAle() / M_PI*xt;
35  return tmp;
36 }
37 
38 double EWSMTwoLoopQCD::DeltaRho(const double Mw_i) const
39 {
40  double Mw = Mw_i;
41  return ( 3.0 * cache.Xt_alpha(Mw) * cache.alsMt() / M_PI * deltaQCD_2());
42 }
43 
44 double EWSMTwoLoopQCD::DeltaR_rem(const double Mw_i) const
45 {
46  double Mw = Mw_i;
47  return ( (2.0 * DeltaR_ud(Mw) + DeltaR_tb(Mw))
48  + cache.getSM().cW2(Mw) / cache.getSM().sW2(Mw) * DeltaRho(Mw));
49 }
50 
51 gslpp::complex EWSMTwoLoopQCD::deltaRho_rem_f(const Particle f, const double Mw_i) const
52 {
53  if (f.is("TOP")) return ( gslpp::complex(0.0, 0.0, false));
54  double Mw = Mw_i;
55  return ( (2.0 * DeltaRho_ud(Mw) + DeltaRho_tb(Mw)) - DeltaRho(Mw));
56 }
57 
58 gslpp::complex EWSMTwoLoopQCD::deltaKappa_rem_f(const Particle f, const double Mw_i) const
59 {
60  if (f.is("TOP")) return ( gslpp::complex(0.0, 0.0, false));
61  double Mw = Mw_i;
62  return ( (2.0 * DeltaKappa_ud(Mw) + DeltaKappa_tb(Mw))
63  - cache.getSM().cW2(Mw) / cache.getSM().sW2(Mw) * DeltaRho(Mw));
64 }
65 
66 
68 
69 double EWSMTwoLoopQCD::deltaQCD_2() const
70 {
71  return ( -2.0 / 3.0 * (1.0 + 2.0 * cache.getZeta2()));
72 }
73 
74 double EWSMTwoLoopQCD::F1(const double x, const double Mw_i) const
75 {
76  if (x < 0.0 || x >= 1.0) throw std::runtime_error("x is out of range in EWSMTwoLoopQCD::F1");
77 
78  /* Zeta functions */
79  double zeta_2 = cache.getZeta2();
80  double zeta_3 = cache.getZeta3();
81 
82  if (x == 0.0) return (23.0 / 16.0 - zeta_2 / 2.0 - 3.0 / 2.0 * zeta_3);
83 
84  /* Dilogarithm and Trilogarithm */
85  double Li2_x, Li3_x, Li3_mx_1mx;
86  double Mw = Mw_i;
87  double Mt = cache.getSM().getMtpole();
88  if (x == Mw * Mw / Mt / Mt) {
89  Li2_x = cache.Li2_MW2toMTOP2(Mw);
90  Li3_x = cache.Li3_MW2toMTOP2(Mw);
91  Li3_mx_1mx = cache.Li3_for_F1(Mw);
92  } else {
93  Li2_x = cache.getPolyLog().Li2(x).real(); // x <= 1.0
94  Li3_x = cache.getPolyLog().Li3(x);
95  Li3_mx_1mx = cache.getPolyLog().Li3(-x / (1.0 - x));
96  }
97 
98  double b = log(1.0 - x);
99 
100  double F1;
101  F1 = (x - 3.0 / 2.0 + 1.0 / 2.0 / x / x)
102  *(-Li3_x - Li3_mx_1mx + b / 3.0 * (2.0 * Li2_x - zeta_2) + b * b * b / 6.0)
103  + 1.0 / 3.0 * (x + 1.0 / 2.0 - 1.0 / 2.0 / x) * Li2_x
104  + b * b / 6.0 * (x - 3.0 / 4.0 - 3.0 / 2.0 / x + 5.0 / 4.0 / x / x)
105  - b / 4.0 * (x - 5.0 / 2.0 + 2.0 / 3.0 / x + 5.0 / 6.0 / x / x)
106  + zeta_3 * (x - 3.0 / 2.0) + zeta_2 / 3.0 * (x - 7.0 / 4.0 - 1.0 / 2.0 / x)
107  + 13.0 / 12.0 - 5.0 / 24.0 / x;
108  return F1;
109 }
110 
111 double EWSMTwoLoopQCD::V1(const double r) const
112 {
113  if (r < 0.0 || r >= 1.0) throw std::runtime_error("r is out of range in EWSMTwoLoopQCD::V1");
114 
115  /* Zeta functions */
116  double zeta_3 = cache.getZeta3();
117 
118  if (r == 0.0) return (0.0);
119 
120  double Mz = cache.getSM().getMz();
121  double Mt = cache.getSM().getMtpole();
122 
123  /* Logarithms etc */
124  double Phi, gamma, h;
125  if (r == Mz * Mz / 4.0 / Mt / Mt) {
126  Phi = cache.Phi_QCD2();
127  gamma = cache.gamma_QCD2();
128  h = cache.h_QCD2();
129  } else {
130  Phi = asin(sqrt(r));
131  gamma = log(2.0 * sqrt(r));
132  h = log(2.0 * sqrt(1.0 - r));
133  }
134 
135  /* Clausen functions */
136  double Cl3_2Phi, Cl3_4Phi, Cl2_2Phi, Cl2_4Phi;
137  if (r == Mz * Mz / 4.0 / Mt / Mt) {
138  Cl3_2Phi = cache.Cl3_2Phi();
139  Cl3_4Phi = cache.Cl3_4Phi();
140  Cl2_2Phi = cache.Cl2_2Phi();
141  Cl2_4Phi = cache.Cl2_4Phi();
142  } else {
143  ClausenFunctions* myClausen;
144  myClausen = new ClausenFunctions();
145  Cl3_2Phi = myClausen->Cl3(2.0 * Phi);
146  Cl3_4Phi = myClausen->Cl3(4.0 * Phi);
147  Cl2_2Phi = myClausen->Cl2(2.0 * Phi);
148  Cl2_4Phi = myClausen->Cl2(4.0 * Phi);
149  delete myClausen;
150  }
151 
152  double V1;
153  V1 = 4.0 * (r - 1.0 / 4.0 / r)
154  *(2.0 * Cl3_2Phi - Cl3_4Phi + 8.0 / 3.0 * Phi * (Cl2_2Phi - Cl2_4Phi)
155  - 4.0 / 3.0 * Phi * Phi * (gamma + 2.0 * h))
156  + sqrt(1.0 / r - 1.0)*(8.0 / 3.0 * (r + 1.0 / 2.0)
157  *(-Cl2_2Phi + Cl2_4Phi + 2.0 * Phi * (gamma + 2.0 * h))
158  - 2.0 * Phi * (r + 3.0 / 2.0))
159  + 8.0 * Phi * Phi * (r - 1.0 / 6.0 - 7.0 / 48.0 / r) + 13.0 / 6.0 + zeta_3 / r;
160  return V1;
161 }
162 
163 double EWSMTwoLoopQCD::A1(const double r) const
164 {
165  if (r < 0.0 || r >= 1.0) throw std::runtime_error("r is out of range in EWSMTwoLoopQCD::A1");
166 
167  /* Zeta functions */
168  double zeta_2 = cache.getZeta2();
169  double zeta_3 = cache.getZeta3();
170 
171  if (r == 0.0) return (3.0 * (7.0 / 4.0 - zeta_2 - 2.0 * zeta_3));
172 
173  double Mz = cache.getSM().getMz();
174  double Mt = cache.getSM().getMtpole();
175 
176  /* Logarithms etc */
177  double Phi, gamma, h;
178  if (r == Mz * Mz / 4.0 / Mt / Mt) {
179  Phi = cache.Phi_QCD2();
180  gamma = cache.gamma_QCD2();
181  h = cache.h_QCD2();
182  } else {
183  Phi = asin(sqrt(r));
184  gamma = log(2.0 * sqrt(r));
185  h = log(2.0 * sqrt(1.0 - r));
186  }
187 
188  /* Clausen functions */
189  double Cl3_2Phi, Cl3_4Phi, Cl2_2Phi, Cl2_4Phi;
190  if (r == Mz * Mz / 4.0 / Mt / Mt) {
191  Cl3_2Phi = cache.Cl3_2Phi();
192  Cl3_4Phi = cache.Cl3_4Phi();
193  Cl2_2Phi = cache.Cl2_2Phi();
194  Cl2_4Phi = cache.Cl2_4Phi();
195  } else {
196  ClausenFunctions* myClausen;
197  myClausen = new ClausenFunctions();
198  Cl3_2Phi = myClausen->Cl3(2.0 * Phi);
199  Cl3_4Phi = myClausen->Cl3(4.0 * Phi);
200  Cl2_2Phi = myClausen->Cl2(2.0 * Phi);
201  Cl2_4Phi = myClausen->Cl2(4.0 * Phi);
202  delete myClausen;
203  }
204 
205  double A1;
206  A1 = 4.0 * (r - 3.0 / 2.0 + 1.0 / 2.0 / r)*(2.0 * Cl3_2Phi - Cl3_4Phi
207  + 8.0 / 3.0 * Phi * (Cl2_2Phi - Cl2_4Phi) - 4.0 / 3.0 * Phi * Phi * (gamma + 2.0 * h))
208  + sqrt(1.0 / r - 1.0)*(8.0 / 3.0 * (r - 1.0)
209  *(-Cl2_2Phi + Cl2_4Phi + 2.0 * Phi * (gamma + 2.0 * h))
210  - 2.0 * Phi * (r - 3.0 + 1.0 / 4.0 / r))
211  + 8.0 * Phi * Phi * (r - 11.0 / 12.0 + 5.0 / 48.0 / r + 1.0 / 32.0 / r / r)
212  - 3.0 * zeta_2 + 13.0 / 6.0 + (-2.0 * zeta_3 + 1.0 / 4.0) / r;
213  return A1;
214 }
215 
216 double EWSMTwoLoopQCD::V1prime(const double r) const
217 {
218  if (r < 0.0 || r >= 1.0) throw std::runtime_error("r is out of range in EWSMTwoLoopQCD::V1prime");
219 
220  /* Zeta functions */
221  double zeta_3 = cache.getZeta3();
222 
223  if (r == 0.0) return (4.0 * zeta_3 - 5.0 / 6.0);
224 
225  double Mz = cache.getSM().getMz();
226  double Mt = cache.getSM().getMtpole();
227 
228  /* Logarithms etc */
229  double Phi, gamma, h;
230  if (r == Mz * Mz / 4.0 / Mt / Mt) {
231  Phi = cache.Phi_QCD2();
232  gamma = cache.gamma_QCD2();
233  h = cache.h_QCD2();
234  } else {
235  Phi = asin(sqrt(r));
236  gamma = log(2.0 * sqrt(r));
237  h = log(2.0 * sqrt(1.0 - r));
238  }
239 
240  /* Clausen functions */
241  double Cl3_2Phi, Cl3_4Phi, Cl2_2Phi, Cl2_4Phi;
242  if (r == Mz * Mz / 4.0 / Mt / Mt) {
243  Cl3_2Phi = cache.Cl3_2Phi();
244  Cl3_4Phi = cache.Cl3_4Phi();
245  Cl2_2Phi = cache.Cl2_2Phi();
246  Cl2_4Phi = cache.Cl2_4Phi();
247  } else {
248  ClausenFunctions* myClausen;
249  myClausen = new ClausenFunctions();
250  Cl3_2Phi = myClausen->Cl3(2.0 * Phi);
251  Cl3_4Phi = myClausen->Cl3(4.0 * Phi);
252  Cl2_2Phi = myClausen->Cl2(2.0 * Phi);
253  Cl2_4Phi = myClausen->Cl2(4.0 * Phi);
254  delete myClausen;
255  }
256 
257  /* for Bprime and Dprime below */
258  gsl_complex OneMinusE2Iphi = gsl_complex_rect(1.0 - cos(2.0 * Phi), -sin(2.0 * Phi));
259  gsl_complex OneMinusE4Iphi = gsl_complex_rect(1.0 - cos(4.0 * Phi), -sin(4.0 * Phi));
260  double log_real;
261  if (r == Mz * Mz / 4.0 / Mt / Mt) {
262  log_real = cache.logV1primeAndA1prime();
263  } else {
264  log_real = GSL_REAL(gsl_complex_log(OneMinusE2Iphi))
265  - 2.0 * GSL_REAL(gsl_complex_log(OneMinusE4Iphi));
266  }
267 
268  double PhiPrime = 1.0 / 2.0 / sqrt(r * (1.0 - r));
269  double Phi2Prime = Phi / sqrt(r * (1.0 - r));
270  double gammaPrime = 1.0 / 2.0 / r;
271  double hPrime = -1.0 / 2.0 / (1.0 - r);
272 
273  // V1(r) = 4.0*A*B + C*D + E
274  double A = r - 1.0 / 4.0 / r;
275  double Aprime = 1.0 + 1.0 / 4.0 / r / r;
276  double B = 2.0 * Cl3_2Phi - Cl3_4Phi + 8.0 / 3.0 * Phi * (Cl2_2Phi - Cl2_4Phi)
277  - 4.0 / 3.0 * Phi * Phi * (gamma + 2.0 * h);
278  double Bprime = -2.0 / sqrt(r * (1.0 - r))*(Cl2_2Phi - Cl2_4Phi)
279  + 8.0 / 3.0 * PhiPrime * (Cl2_2Phi - Cl2_4Phi)
280  + 8.0 / 3.0 * Phi * (-log_real / sqrt(r * (1.0 - r)))
281  - 4.0 / 3.0 * Phi2Prime * (gamma + 2.0 * h)
282  - 4.0 / 3.0 * Phi * Phi * (gammaPrime + 2.0 * hPrime);
283  double C = sqrt(1.0 / r - 1.0);
284  double Cprime = -1.0 / 2.0 / r / sqrt(r * (1.0 - r));
285  double D = 8.0 / 3.0 * (r + 1.0 / 2.0)
286  *(-Cl2_2Phi + Cl2_4Phi + 2.0 * Phi * (gamma + 2.0 * h))
287  - 2.0 * Phi * (r + 3.0 / 2.0);
288  double Dprime = 8.0 / 3.0 * (-Cl2_2Phi + Cl2_4Phi + 2.0 * Phi * (gamma + 2.0 * h))
289  + 8.0 / 3.0 * (r + 1.0 / 2.0)
290  *(log_real / sqrt(r * (1.0 - r))
291  + 2.0 * PhiPrime * (gamma + 2.0 * h)
292  + 2.0 * Phi * (gammaPrime + 2.0 * hPrime))
293  - 2.0 * PhiPrime * (r + 3.0 / 2.0) - 2.0 * Phi;
294  double Eprime = 8.0 * Phi2Prime * (r - 1.0 / 6.0 - 7.0 / 48.0 / r)
295  + 8.0 * Phi * Phi * (1.0 + 7.0 / 48.0 / r / r) - zeta_3 / r / r;
296  return (4.0 * Aprime * B + 4.0 * A * Bprime + Cprime * D + C * Dprime + Eprime);
297 
298 
299  /* TEST: Exact - Expansion */
300  //std::cout << "V1: "
301  // << (4.0*Aprime*B + 4.0*A*Bprime + Cprime*D + C*Dprime + Eprime)
302  // - (4.0*zeta_3 - 5.0/6.0 + 656.0/81.0*r)
303  // << std::endl;
304 
305  /* Expansion for r << 1 */
306  //return (4.0*zeta_3 - 5.0/6.0 + 656.0/81.0*r);
307 }
308 
309 double EWSMTwoLoopQCD::A1prime(const double r) const
310 {
311  if (r < 0.0 || r >= 1.0) throw std::runtime_error("r is out of range in EWSMTwoLoopQCD::A1prime");
312 
313  /* Zeta functions */
314  double zeta_2 = cache.getZeta2();
315  double zeta_3 = cache.getZeta3();
316 
317  if (r == 0.0) return (3.0 * (7.0 / 4.0 - zeta_2 - 2.0 * zeta_3));
318 
319  double Mz = cache.getSM().getMz();
320  double Mt = cache.getSM().getMtpole();
321 
322  /* Logarithms etc */
323  double Phi, gamma, h;
324  if (r == Mz * Mz / 4.0 / Mt / Mt) {
325  Phi = cache.Phi_QCD2();
326  gamma = cache.gamma_QCD2();
327  h = cache.h_QCD2();
328  } else {
329  Phi = asin(sqrt(r));
330  gamma = log(2.0 * sqrt(r));
331  h = log(2.0 * sqrt(1.0 - r));
332  }
333 
334  /* Clausen functions */
335  double Cl3_2Phi, Cl3_4Phi, Cl2_2Phi, Cl2_4Phi;
336  if (r == Mz * Mz / 4.0 / Mt / Mt) {
337  Cl3_2Phi = cache.Cl3_2Phi();
338  Cl3_4Phi = cache.Cl3_4Phi();
339  Cl2_2Phi = cache.Cl2_2Phi();
340  Cl2_4Phi = cache.Cl2_4Phi();
341  } else {
342  ClausenFunctions* myClausen;
343  myClausen = new ClausenFunctions();
344  Cl3_2Phi = myClausen->Cl3(2.0 * Phi);
345  Cl3_4Phi = myClausen->Cl3(4.0 * Phi);
346  Cl2_2Phi = myClausen->Cl2(2.0 * Phi);
347  Cl2_4Phi = myClausen->Cl2(4.0 * Phi);
348  delete myClausen;
349  }
350 
351  /* for Bprime and Dprime below */
352  gsl_complex OneMinusE2Iphi = gsl_complex_rect(1.0 - cos(2.0 * Phi), -sin(2.0 * Phi));
353  gsl_complex OneMinusE4Iphi = gsl_complex_rect(1.0 - cos(4.0 * Phi), -sin(4.0 * Phi));
354  double log_real;
355  if (r == Mz * Mz / 4.0 / Mt / Mt) {
356  log_real = cache.logV1primeAndA1prime();
357  } else {
358  log_real = GSL_REAL(gsl_complex_log(OneMinusE2Iphi))
359  - 2.0 * GSL_REAL(gsl_complex_log(OneMinusE4Iphi));
360  }
361 
362  double PhiPrime = 1.0 / 2.0 / sqrt(r * (1.0 - r));
363  double Phi2Prime = Phi / sqrt(r * (1.0 - r));
364  double gammaPrime = 1.0 / 2.0 / r;
365  double hPrime = -1.0 / 2.0 / (1.0 - r);
366 
367  // A1(r) = 4.0*A*B + C*D + E
368  double A = r - 3.0 / 2.0 + 1.0 / 2.0 / r;
369  double Aprime = 1.0 - 1.0 / 2.0 / r / r;
370  double B = 2.0 * Cl3_2Phi - Cl3_4Phi + 8.0 / 3.0 * Phi * (Cl2_2Phi - Cl2_4Phi)
371  - 4.0 / 3.0 * Phi * Phi * (gamma + 2.0 * h);
372  double Bprime = -2.0 / sqrt(r * (1.0 - r))*(Cl2_2Phi - Cl2_4Phi)
373  + 8.0 / 3.0 * PhiPrime * (Cl2_2Phi - Cl2_4Phi)
374  + 8.0 / 3.0 * Phi * (-log_real / sqrt(r * (1.0 - r)))
375  - 4.0 / 3.0 * Phi2Prime * (gamma + 2.0 * h)
376  - 4.0 / 3.0 * Phi * Phi * (gammaPrime + 2.0 * hPrime);
377  double C = sqrt(1.0 / r - 1.0);
378  double Cprime = -1.0 / 2.0 / r / sqrt(r * (1.0 - r));
379  double D = 8.0 / 3.0 * (r - 1.0)*(-Cl2_2Phi + Cl2_4Phi + 2.0 * Phi * (gamma + 2.0 * h))
380  - 2.0 * Phi * (r - 3.0 + 1.0 / 4.0 / r);
381  double Dprime = 8.0 / 3.0 * (-Cl2_2Phi + Cl2_4Phi + 2.0 * Phi * (gamma + 2.0 * h))
382  + 8.0 / 3.0 * (r - 1.0)*(log_real / sqrt(r * (1.0 - r))
383  + 2.0 * PhiPrime * (gamma + 2.0 * h)
384  + 2.0 * Phi * (gammaPrime + 2.0 * hPrime))
385  - 2.0 * PhiPrime * (r - 3.0 + 1.0 / 4.0 / r)
386  - 2.0 * Phi * (1.0 - 1.0 / 4.0 / r / r);
387  double Eprime = 8.0 * Phi2Prime * (r - 11.0 / 12.0 + 5.0 / 48.0 / r + 1.0 / 32.0 / r / r)
388  + 8.0 * Phi * Phi * (1.0 - 5.0 / 48.0 / r / r - 1.0 / 16.0 / r / r / r)
389  - (-2.0 * zeta_3 + 1.0 / 4.0) / r / r;
390  return (4.0 * Aprime * B + 4.0 * A * Bprime + Cprime * D + C * Dprime + Eprime);
391 
392  /* TEST: Exact - Expansion */
393  //std::cout << "A1:"
394  // << (4.0*Aprime*B + 4.0*A*Bprime + Cprime*D + C*Dprime + Eprime)
395  // - (4.0*zeta_3 - 49.0/18.0 + 1378.0/405.0*r )
396  // << std::endl;
397 
398  /* Expansion for r << 1 */
399  //return (4.0*zeta_3 - 49.0/18.0 + 1378.0/405.0*r);
400 }
401 
402 double EWSMTwoLoopQCD::DeltaR_ud(const double Mw_i) const
403 {
404  double Mw = Mw_i;
405  double sW2 = cache.getSM().sW2(Mw);
406  double cW2 = cache.getSM().cW2(Mw);
407 
408  /* Logarithm */
409  double log_cW2 = cache.log_cW2(Mw);
410 
411  double DeltaR;
412  DeltaR = -log_cW2;
413  DeltaR *= (cW2 - sW2) / 4.0 / sW2 / sW2;
414  DeltaR *= cache.getSM().getAle() * cache.getSM().getAlsMz() / M_PI / M_PI;
415  return DeltaR;
416 }
417 
418 double EWSMTwoLoopQCD::DeltaR_tb(const double Mw_i) const
419 {
420  double Mw = Mw_i;
421  double sW2 = cache.getSM().sW2(Mw);
422  double cW2 = cache.getSM().cW2(Mw);
423  double Mz = cache.getSM().getMz();
424  double Mt = cache.getSM().getMtpole();
425  double wt = Mt * Mt / Mw / Mw;
426  double zt = Mt * Mt / Mz / Mz;
427  double rZ4t = Mz * Mz / 4.0 / Mt / Mt;
428  double xWt = Mw * Mw / Mt / Mt;
429 
430  double vt = cache.v_f(cache.getSM().getQuarks(QCD::TOP), Mw);
431  double at = cache.a_f(cache.getSM().getQuarks(QCD::TOP));
432  double vb = cache.v_f(cache.getSM().getQuarks(QCD::BOTTOM), Mw);
433  double ab = cache.a_f(cache.getSM().getQuarks(QCD::BOTTOM));
434 
435  /* Zeta functions */
436  double zeta_2 = cache.getZeta2();
437 
438  /* Logarithm */
439  double log_zt = -2.0 * cache.logMZtoMTOP();
440 
441  double DeltaR;
442  DeltaR = pow(cache.Q_f(cache.getSM().getQuarks(QCD::TOP)), 2.0) * V1prime(0.0)
443  + cW2 / sW2 / sW2 * wt / 4.0 * (zeta_2 + 1.0 / 2.0)
444  - zt / sW2 / sW2 * (vt * vt * V1(rZ4t) + at * at * (A1(rZ4t) - A1(0.0)))
445  + (cW2 - sW2) / sW2 / sW2 * wt * (F1(xWt, Mw) - F1(0.0, Mw))
446  - vb * ab / 2.0 / sW2 / sW2*log_zt;
447  DeltaR *= cache.getSM().getAle() * cache.alsMt() / M_PI / M_PI;
448  return DeltaR;
449 }
450 
451 double EWSMTwoLoopQCD::DeltaRho_ud(const double Mw_i) const
452 {
453  double Mw = Mw_i;
454  double sW2 = cache.getSM().sW2(Mw);
455  double cW2 = cache.getSM().cW2(Mw);
456 
457  double DeltaRho;
458  DeltaRho = pow(cache.v_f(cache.getSM().getQuarks(QCD::UP), Mw), 2.0)
459  + pow(cache.v_f(cache.getSM().getQuarks(QCD::DOWN), Mw), 2.0)
460  + pow(cache.a_f(cache.getSM().getQuarks(QCD::UP)), 2.0)
461  + pow(cache.a_f(cache.getSM().getQuarks(QCD::TOP)), 2.0);
462  DeltaRho /= 4.0 * sW2*cW2;
463  DeltaRho *= cache.getSM().getAle() * cache.getSM().getAlsMz() / M_PI / M_PI;
464  return DeltaRho;
465 }
466 
467 double EWSMTwoLoopQCD::DeltaRho_tb(const double Mw_i) const
468 {
469  double Mw = Mw_i;
470  double Mz = cache.getSM().getMz();
471  double sW2 = cache.getSM().sW2(Mw);
472  double cW2 = cache.getSM().cW2(Mw);
473  double Mt = cache.getSM().getMtpole();
474  double zt = Mt * Mt / Mz / Mz;
475  double rZ4t = Mz * Mz / 4.0 / Mt / Mt;
476 
477  double vt = cache.v_f(cache.getSM().getQuarks(QCD::TOP), Mw);
478  double at = cache.a_f(cache.getSM().getQuarks(QCD::TOP));
479  double vb = cache.v_f(cache.getSM().getQuarks(QCD::BOTTOM), Mw);
480  double ab = cache.a_f(cache.getSM().getQuarks(QCD::BOTTOM));
481 
482  double DeltaRho;
483  DeltaRho = -(vt * vt * V1prime(rZ4t) + at * at * A1prime(rZ4t))
484  + 4.0 * zt * (vt * vt * V1(rZ4t) + at * at * A1(rZ4t))
485  + vb * vb + ab * ab - 4.0 * zt * F1(0.0, Mw);
486  DeltaRho /= 4.0 * sW2*cW2;
487  DeltaRho *= cache.getSM().getAle() * cache.alsMt() / M_PI / M_PI;
488  return DeltaRho;
489 }
490 
491 gslpp::complex EWSMTwoLoopQCD::DeltaKappa_ud(const double Mw_i) const
492 {
493  double Mw = Mw_i;
494  double sW2 = cache.getSM().sW2(Mw);
495  double cW2 = cache.getSM().cW2(Mw);
496 
497  /* Logarithm */
498  double log_cW2 = cache.log_cW2(Mw);
499 
500  gslpp::complex DeltaKappa(0.0, 0.0, false);
501  DeltaKappa = cW2 / 4.0 / sW2 / sW2 * log_cW2
502  + M_PI / 4.0 / sW2 * (1.0 - 20.0 / 9.0 * sW2)*(gslpp::complex::i());
503  DeltaKappa *= cache.getSM().getAle() * cache.getSM().getAlsMz() / M_PI / M_PI;
504  return DeltaKappa;
505 }
506 
507 gslpp::complex EWSMTwoLoopQCD::DeltaKappa_tb(const double Mw_i) const
508 {
509  double Mw = Mw_i;
510  double Mz = cache.getSM().getMz();
511  double sW2 = cache.getSM().sW2(Mw);
512  double cW2 = cache.getSM().cW2(Mw);
513  double Mt = cache.getSM().getMtpole();
514  double wt = Mt * Mt / Mw / Mw;
515  double zt = Mt * Mt / Mz / Mz;
516  double rZ4t = Mz * Mz / 4.0 / Mt / Mt;
517  double xWt = Mw * Mw / Mt / Mt;
518 
519  double vt = cache.v_f(cache.getSM().getQuarks(QCD::TOP), Mw);
520  double at = cache.a_f(cache.getSM().getQuarks(QCD::TOP));
521  double Qt = cache.Q_f(cache.getSM().getQuarks(QCD::TOP));
522  double vb = cache.v_f(cache.getSM().getQuarks(QCD::BOTTOM), Mw);
523  double ab = cache.a_f(cache.getSM().getQuarks(QCD::BOTTOM));
524  double Qb = cache.Q_f(cache.getSM().getQuarks(QCD::BOTTOM));
525 
526  /* Logarithm */
527  double log_zt = -2.0 * cache.logMZtoMTOP();
528 
529  gslpp::complex DeltaKappa(0.0, 0.0, false);
530  DeltaKappa = 4.0 * cW2 * wt * (vt * vt * V1(rZ4t) + at * at * A1(rZ4t) - F1(xWt, Mw))
531  + 4.0 * sW2 * (fabs(Qt) - 4.0 * sW2 * Qt * Qt) * zt * V1(rZ4t)
532  + (vb * vb + ab * ab + sW2 * (fabs(Qb) - 4.0 * sW2 * Qb * Qb)) * log_zt;
533  DeltaKappa += M_PI * sW2 * (1.0 / 3.0 - 4.0 / 9.0 * sW2)*(gslpp::complex::i());
534  DeltaKappa /= 4.0 * sW2*sW2;
535  DeltaKappa *= cache.getSM().getAle() * cache.alsMt() / M_PI / M_PI;
536  return DeltaKappa;
537 }
538 
539 
540 
541 
EWSMcache::gamma_QCD2
double gamma_QCD2() const
The constant for two-loop QCD contribution.
Definition: EWSMcache.h:400
gslpp::cos
complex cos(const complex &z)
Definition: gslpp_complex.cpp:429
StandardModel::cW2
virtual double cW2(const double Mw_i) const
The square of the cosine of the weak mixing angle in the on-shell scheme, denoted as .
Definition: StandardModel.cpp:989
QCD::BOTTOM
Definition: QCD.h:329
Particle::is
bool is(std::string name_i) const
Definition: Particle.cpp:23
Particle
A class for particles.
Definition: Particle.h:26
EWSMTwoLoopQCD::cache
const EWSMcache & cache
A reference to an object of type EWSMcache.
Definition: EWSMTwoLoopQCD.h:376
EWSMTwoLoopQCD::deltaQCD_2
double deltaQCD_2() const
The function .
Definition: EWSMTwoLoopQCD.cpp:69
gslpp::sin
complex sin(const complex &z)
Definition: gslpp_complex.cpp:420
StandardModel::getAlsMz
double getAlsMz() const
A get method to access the value of .
Definition: StandardModel.h:727
EWSMcache::getZeta3
double getZeta3() const
A get method to access the value of the zeta function .
Definition: EWSMcache.h:146
QCD::UP
Definition: QCD.h:324
Polylogarithms::Li3
double Li3(const double x) const
The trilogarithm .
Definition: Polylogarithms.cpp:36
Polylogarithms::Li2
gslpp::complex Li2(const double x) const
The dilogarithm with a real argument, .
Definition: Polylogarithms.cpp:22
EWSMTwoLoopQCD::DeltaRho_ud
double DeltaRho_ud(const double Mw_i) const
Light-quark contribution to , denoted as .
Definition: EWSMTwoLoopQCD.cpp:451
gslpp::complex
A class for defining operations on and functions of complex numbers.
Definition: gslpp_complex.h:35
EWSMcache::Cl2_2Phi
double Cl2_2Phi() const
The constant .
Definition: EWSMcache.h:438
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
EWSMcache::Li3_for_F1
double Li3_for_F1(const double Mw_i) const
A cache method.
Definition: EWSMcache.cpp:182
EWSMTwoLoopQCD::DeltaRho_tb
double DeltaRho_tb(const double Mw_i) const
Heavy-quark contribution to , denoted as .
Definition: EWSMTwoLoopQCD.cpp:467
StandardModel::sW2
virtual double sW2(const double Mw_i) const
The square of the sine of the weak mixing angle in the on-shell scheme, denoted as .
Definition: StandardModel.cpp:1000
EWSMTwoLoopQCD::DeltaR_rem
double DeltaR_rem(const double Mw_i) const
Remainder contribution of to , denoted as .
Definition: EWSMTwoLoopQCD.cpp:44
EWSMcache::Li3_MW2toMTOP2
double Li3_MW2toMTOP2(const double Mw_i) const
A cache method.
Definition: EWSMcache.cpp:168
EWSMcache::getSM
const StandardModel & getSM() const
Definition: EWSMcache.h:56
ClausenFunctions::Cl3
double Cl3(const double phi) const
The Clausen function of index 3, .
Definition: ClausenFunctions.cpp:26
EWSMcache::Cl3_2Phi
double Cl3_2Phi() const
The constant .
Definition: EWSMcache.h:462
EWSMcache::Phi_QCD2
double Phi_QCD2() const
The constant for two-loop QCD contribution.
Definition: EWSMcache.h:390
EWSMTwoLoopQCD::DeltaKappa_ud
gslpp::complex DeltaKappa_ud(const double Mw_i) const
Light-quark contribution to , denoted as .
Definition: EWSMTwoLoopQCD.cpp:491
ClausenFunctions::Cl2
double Cl2(const double phi) const
The Clausen function of index 2, .
Definition: ClausenFunctions.cpp:21
EWSMcache::h_QCD2
double h_QCD2() const
The constant for two-loop QCD contribution.
Definition: EWSMcache.h:410
QCD::TOP
Definition: QCD.h:328
EWSMTwoLoopQCD::V1
double V1(const double r) const
The function .
Definition: EWSMTwoLoopQCD.cpp:111
gslpp::pow
complex pow(const complex &z1, const complex &z2)
Definition: gslpp_complex.cpp:395
gslpp::sqrt
complex sqrt(const complex &z)
Definition: gslpp_complex.cpp:385
gslpp::complex::i
static const complex & i()
Definition: gslpp_complex.cpp:154
EWSMTwoLoopQCD::DeltaAlpha_l
double DeltaAlpha_l(const double s) const
Leptonic contribution of to the electromagnetic coupling , denoted as .
Definition: EWSMTwoLoopQCD.cpp:20
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
EWSMTwoLoopQCD::DeltaRho
double DeltaRho(const double Mw_i) const
Leading two-loop QCD contribution of to , denoted as .
Definition: EWSMTwoLoopQCD.cpp:38
EWSMTwoLoopQCD::F1
double F1(const double x, const double Mw_i) const
The function .
Definition: EWSMTwoLoopQCD.cpp:74
EWSMTwoLoopQCD::DeltaR_tb
double DeltaR_tb(const double Mw_i) const
Heavy-quark contribution to , denoted as .
Definition: EWSMTwoLoopQCD.cpp:418
EWSMcache::a_f
double a_f(const Particle f) const
The tree-level axial-vector coupling for , denoted as .
Definition: EWSMcache.h:301
QCD::getQuarks
Particle getQuarks(const QCD::quark q) const
A get method to access a quark as an object of the type Particle.
Definition: QCD.h:534
EWSMTwoLoopQCD::deltaRho_rem_f
gslpp::complex deltaRho_rem_f(const Particle f, const double Mw_i) const
Remainder contribution of to the effective couplings , denoted as .
Definition: EWSMTwoLoopQCD.cpp:51
EWSMcache::getPolyLog
const Polylogarithms getPolyLog() const
A get method to access the member object of type Polylogarithms.
Definition: EWSMcache.h:115
EWSMTwoLoopQCD::DeltaR_ud
double DeltaR_ud(const double Mw_i) const
Light-quark contribution to , not including , denoted as .
Definition: EWSMTwoLoopQCD.cpp:402
EWSMcache::Als
double Als(const double mu, const orders order) const
The strong coupling .
Definition: EWSMcache.h:366
EWSMcache::Li2_MW2toMTOP2
double Li2_MW2toMTOP2(const double Mw_i) const
A cache method.
Definition: EWSMcache.cpp:154
EWSMcache::logMZtoMTOP
double logMZtoMTOP() const
A cache method.
Definition: EWSMcache.cpp:112
EWSMcache::logV1primeAndA1prime
double logV1primeAndA1prime() const
A logarithm appearing in the functions and for two-loop QCD contribution.
Definition: EWSMcache.h:422
EWSMcache::getZeta2
double getZeta2() const
A get method to access the value of the zeta function .
Definition: EWSMcache.h:137
EWSMcache::log_cW2
double log_cW2(const double Mw_i) const
A cache method.
Definition: EWSMcache.cpp:140
EWSMcache::Cl2_4Phi
double Cl2_4Phi() const
The constant .
Definition: EWSMcache.h:450
StandardModel::getMz
double getMz() const
A get method to access the mass of the boson .
Definition: StandardModel.h:718
EWSMTwoLoopQCD::EWSMTwoLoopQCD
EWSMTwoLoopQCD(const EWSMcache &cache_i)
Constructor.
Definition: EWSMTwoLoopQCD.cpp:12
EWSMTwoLoopQCD::DeltaKappa_tb
gslpp::complex DeltaKappa_tb(const double Mw_i) const
Heavy-quark contribution to , denoted as .
Definition: EWSMTwoLoopQCD.cpp:507
Mw
An observable class for the -boson mass.
Definition: Mw.h:22
ClausenFunctions
A class for the Clausen functions.
Definition: ClausenFunctions.h:22
EWSMTwoLoopQCD::V1prime
double V1prime(const double r) const
The derivative of the function .
Definition: EWSMTwoLoopQCD.cpp:216
gslpp::complex::real
const double & real() const
Definition: gslpp_complex.cpp:53
EWSMTwoLoopQCD::A1
double A1(const double r) const
The function .
Definition: EWSMTwoLoopQCD.cpp:163
EWSMTwoLoopQCD::A1prime
double A1prime(const double r) const
The derivative of the function .
Definition: EWSMTwoLoopQCD.cpp:309
EWSMcache::Cl3_4Phi
double Cl3_4Phi() const
The constant .
Definition: EWSMcache.h:474
EWSMTwoLoopQCD::deltaKappa_rem_f
gslpp::complex deltaKappa_rem_f(const Particle f, const double Mw_i) const
Remainder contribution of to the effective couplings , denoted as .
Definition: EWSMTwoLoopQCD.cpp:58
EWSMcache::Q_f
double Q_f(const Particle f) const
The charge of an SM fermion .
Definition: EWSMcache.h:268
EWSMcache::alsMt
double alsMt() const
The strong coupling at NNLO.
Definition: EWSMcache.h:378
FULLNNLO
Definition: OrderScheme.h:38
QCD::DOWN
Definition: QCD.h:325
EWSMcache::Xt_alpha
double Xt_alpha(const double Mw_i) const
The quantity with the coupling .
Definition: EWSMcache.h:355
EWSMTwoLoopQCD.h
StandardModel::getAle
double getAle() const
A get method to retrieve the fine-structure constant .
Definition: StandardModel.h:745
EWSMcache::v_f
double v_f(const Particle f, const double Mw_i) const
The tree-level vector coupling for , denoted as .
Definition: EWSMcache.h:290
EWSMTwoLoopQCD::DeltaAlpha_t
double DeltaAlpha_t(const double s) const
Top-quark contribution of to the electromagnetic coupling , denoted as .
Definition: EWSMTwoLoopQCD.cpp:25