a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
EWSMTwoLoopEW.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 <TMath.h>
10#include "EWSMTwoLoopEW.h"
11
12/* include O(alpha^2 Mt2/Mz2) in addition to O(alpha^2 M_t^4/M_Z^4) */
13#define EW_SUBLEADING_ALPHA2
14
16: cache(cache_i), myOneLoopEW(cache_i)
17{
18}
19
20
22
23double EWSMTwoLoopEW::DeltaAlpha_l(const double s) const
24{
28 double log_l[3];
29 if (s == cache.getSM().getMz() * cache.getSM().getMz()) {
30 log_l[0] = 2.0 * cache.logMZtoME();
31 log_l[1] = 2.0 * cache.logMZtoMMU();
32 log_l[2] = 2.0 * cache.logMZtoMTAU();
33 } else {
34 log_l[0] = log(xl[0]);
35 log_l[1] = log(xl[1]);
36 log_l[2] = log(xl[2]);
37 }
38
39 double twoLoop[3];
40 for (int i = 0; i < 3; i++) {
41 twoLoop[i] = -5.0 / 24.0 + cache.getZeta3() + log_l[i] / 4.0
42 + 3.0 / xl[i] * log_l[i];
43 }
44
45 return ( pow(cache.getSM().getAle() / M_PI, 2.0)
46 *(twoLoop[0] + twoLoop[1] + twoLoop[2]));
47}
48
49double EWSMTwoLoopEW::DeltaAlpha_t(const double s) const
50{
51 return (0.0);
52}
53
54double EWSMTwoLoopEW::DeltaRho(const double Mw_i) const
55{
56 double Mz = cache.getSM().getMz();
57 double Mw = Mw_i;
58 double sW2 = cache.getSM().sW2(Mw);
59 double cW2 = cache.getSM().cW2(Mw);
60
61 double DeltaRho = 0.0;
62
63#ifndef EW_SUBLEADING_ALPHA2
64 /* O(\alpha^2 Mt^4/Mz^4) */
65 DeltaRho += 3.0 * rho_2();
66 DeltaRho *= pow(cache.Xt_alpha(Mw), 2.0);
67#else
68 /* O(\alpha^2 Mt^4/Mz^4 + \alpha^2 Mt^2/Mz^2) */
69 double zt = Mz * Mz / cache.getSM().getMtpole() / cache.getSM().getMtpole();
70 DeltaRho += 3.0 * pow(cache.Xt_alpha(Mw), 2.0)
71 *(DeltaRho2(Mw) + 4.0 * zt * cW2 * DeltaRho2Add(Mw));
72#endif
73
74 /* add O(alpha^2) contribution from the Z-gamma mixing */
75 DeltaRho += -pow(cache.getSM().getAle() / 4.0 / M_PI, 2.0) * cW2 / sW2
76 * pow(myOneLoopEW.PibarZgamma_fer(Mz, Mz*Mz, Mw).real(), 2.0);
77
78 return DeltaRho;
79}
80
81double EWSMTwoLoopEW::DeltaR_rem(const double Mw_i) const
82{
83 double DeltaRrem = 0.0;
84
85#ifdef EW_SUBLEADING_ALPHA2
86 /* O(\alpha^2 Mt^4/Mz^4 + \alpha^2 Mt^2/Mz^2) */
87 double Mw = Mw_i;
88 double sW2 = cache.getSM().sW2(Mw);
89 DeltaRrem += 3.0 * pow(cache.getSM().getAle() * cache.getSM().getMtpole() / 4.0 / M_PI / sW2 / Mw, 2.0)
90 *(DeltaRw2(Mw) + sW2 * deltaEoverE2() + f2Add(Mw) / 4.0);
91#endif
92
93 return DeltaRrem;
94}
95
96gslpp::complex EWSMTwoLoopEW::deltaRho_rem_f(const Particle f, const double Mw_i) const
97{
98 if (f.is("TOP")) return ( gslpp::complex(0.0, 0.0, false));
99
100 gslpp::complex dRho = gslpp::complex(0.0, 0.0, false);
101
102#ifdef EW_SUBLEADING_ALPHA2
103 /* O(\alpha^2 Mt^4/Mz^4 + \alpha^2 Mt^2/Mz^2) */
104 double Mz = cache.getSM().getMz();
105 double Mw = Mw_i;
106 double cW2 = cache.getSM().cW2(Mw);
107 double zt = Mz * Mz / cache.getSM().getMtpole() / cache.getSM().getMtpole();
108 dRho += 3.0 * pow(cache.Xt_alpha(Mw), 2.0)
109 *(16.0 * zt * cW2 * DeltaEta2(Mw) + 4.0 * zt * cW2 * DeltaEta2Add_f(f, Mw));
110#endif
111
112 return dRho;
113}
114
115gslpp::complex EWSMTwoLoopEW::deltaKappa_rem_f(const Particle f, const double Mw_i) const
116{
117 if (f.is("TOP")) return ( gslpp::complex(0.0, 0.0, false));
118
119 gslpp::complex dKappa = gslpp::complex(0.0, 0.0, false);
120
121#ifdef EW_SUBLEADING_ALPHA2
122 /* O(\alpha^2 Mt^4/Mz^4 + \alpha^2 Mt^2/Mz^2) */
123 double Mz = cache.getSM().getMz();
124 double Mw = Mw_i;
125 double cW2 = cache.getSM().cW2(Mw);
126 double zt = Mz * Mz / cache.getSM().getMtpole() / cache.getSM().getMtpole();
127 dKappa += 3.0 * pow(cache.Xt_alpha(Mw), 2.0)
128 *(16.0 * zt * cW2 * DeltaKappa2(Mw) + 4.0 * zt * cW2 * DeltaKappa2Add_f(f, Mw));
129#endif
130
131 return dKappa;
132}
133
134
136
138{
139 double a = cache.getSM().getMHl() * cache.getSM().getMHl() / cache.getSM().getMtpole() / cache.getSM().getMtpole();
140 if (a <= 0.0) throw std::runtime_error("a is out of range in EWSMTwoLoopEW::rho_2");
141 double g_a = g(a);
142 double f_a_0 = f0(a); // f(a,0)
143 double f_a_1 = f1(a); // f(a,1)
144 double log_a = -2.0 * cache.logMTOPtoMH();
145 return ( 25.0 - 4.0 * a + 0.5 * (a * a - 12.0 * a - 12.0) * log_a
146 + (a - 2.0) / 2.0 / a * M_PI * M_PI + 0.5 * (a - 4.0) * sqrt(a) * g_a
147 - 3.0 / a * (a - 1.0)*(a - 1.0)*(a - 2.0) * f_a_0
148 + 3.0 * (a * a - 6.0 * a + 10.0) * f_a_1);
149}
150
152{
153 double a = cache.getSM().getMHl() * cache.getSM().getMHl() / cache.getSM().getMtpole() / cache.getSM().getMtpole();
154 if (a <= 0.0) throw std::runtime_error("a is out of range in EWSMTwoLoopEW::tau_2");
155 double g_a = g(a);
156 double f_a_0 = f0(a); // f(a,0)
157 double f_a_1 = f1(a); // f(a,1)
158 double log_a = -2.0 * cache.logMTOPtoMH();
159 return ( 9.0 - 13.0 / 4.0 * a - 2.0 * a * a - a / 4.0 * (19.0 + 6.0 * a) * log_a
160 - a * a / 4.0 * (7.0 - 6.0 * a) * log_a * log_a
161 - (1.0 / 4.0 + 7.0 / 2.0 * a * a - 3.0 * a * a * a) * M_PI * M_PI / 6.0
162 + (a / 2.0 - 2.0) * sqrt(a) * g_a
163 + (a - 1.0)*(a - 1.0)*(4.0 * a - 7.0 / 4.0) * f_a_0
164 - (a * a * a - 33.0 / 4.0 * a * a + 18.0 * a - 7.0) * f_a_1);
165}
166
167double EWSMTwoLoopEW::g(const double a) const
168{
169 if (a >= 0.0 && a <= 4.0) {
170 double phi = 2.0 * asin(sqrt(a / 4.0));
171 return ( sqrt(4.0 - a)*(M_PI - phi));
172 } else if (a > 4.0) {
173 double y = 4.0 / a;
174 double xi = (sqrt(1.0 - y) - 1.0) / (sqrt(1.0 - y) + 1.0);
175 return ( sqrt(a - 4.0) * log(-xi));
176 } else
177 throw std::runtime_error("Out of range in EWSMTwoLoopEW::g()");
178}
179
180double EWSMTwoLoopEW::f0(const double a) const
181{
182 if (a >= 0.0)
183 return ( cache.getPolyLog().Li2(1.0 - a).real()); // 1-a<1
184 else
185 throw std::runtime_error("Out of range in EWSMTwoLoopEW::f0()");
186}
187
188double EWSMTwoLoopEW::f1(const double a) const
189{
190 if (a >= 0.0 && a <= 4.0) {
191 double y = 4.0 / a;
192 double phi = 2.0 * asin(sqrt(a / 4.0));
193 return ( -2.0 / sqrt(y - 1.0) * cache.getClausen().Cl2(phi));
194 } else if (a > 4.0) {
195 double y = 4.0 / a; // 0<y<1
196 double xi = (sqrt(1.0 - y) - 1.0) / (sqrt(1.0 - y) + 1.0); // -1<xi<0
197 return ( -1.0 / sqrt(1.0 - y)*(cache.getPolyLog().Li2(xi).real()
198 - cache.getPolyLog().Li2(1.0 / xi).real()));
199 } else
200 throw std::runtime_error("Out of range in EWSMTwoLoopEW::f1()");
201}
202
203double EWSMTwoLoopEW::DeltaRho2(const double Mw_i) const
204{
205 double Mt = cache.getSM().getMtpole();
206 double mh = cache.getSM().getMHl();
207 double ht = mh * mh / Mt / Mt;
208
209 double rho2;
210 if (mh < 3.8 * Mt) {
211 rho2 = -15.642 + 0.036382 * Mt + pow(ht, 1.0 / 4.0)*(2.301 - 0.01343 * Mt)
212 + sqrt(ht)*(0.01809 * Mt - 9.953) + ht * (5.687 - 0.01568 * Mt)
213 + pow(ht, 3.0 / 2.0)*(0.005369 * Mt - 1.647)
214 + ht * ht * (0.1852 - 0.000646 * Mt);
215 } else {
216 double Mz = cache.getSM().getMz();
217 double Mw = Mw_i;
218 double Mz2 = Mz*Mz, Mw2 = Mw*Mw, Mt2 = Mt*Mt;
219 double cW2 = cache.getSM().cW2(Mw), sW2 = cache.getSM().sW2(Mw);
220 double cW4 = cW2*cW2, cW6 = cW4*cW2;
221 double zt = Mz * Mz / Mt / Mt;
222 double ht2 = ht*ht, ht3 = ht2*ht, ht4 = ht3*ht, ht5 = ht4*ht, ht6 = ht5*ht;
223 double mu = Mt; // renormalization scale
224
225 double B0_Mt2_Mz2_Mw2_Mw2 = cache.getPV().B0(Mt2, Mz2, Mw2, Mw2).real();
226 double B0_Mt2_Mw2_Mw2_Mz2 = cache.getPV().B0(Mt2, Mw2, Mw2, Mz2).real();
227 double Li21mht = cache.getPolyLog().Li2(1.0 - ht).real();
228
229 return ( 25.0 - 4.0 * ht + (1.0 / 2.0 - 1.0 / ht) * M_PI * M_PI
230 + (ht - 4.0) * sqrt(ht) * g(ht) / 2.0
231 + (-6.0 - 6.0 * ht + ht2 / 2.0) * log(ht)
232 + (6.0 / ht - 15.0 + 12.0 * ht - 3.0 * ht2) * Li21mht
233 + 3.0 / 2.0 * (-10.0 + 6.0 * ht - ht2) * phi(ht / 4.0)
234 + zt * (1.0 / (54.0 * cW2 * (ht - 4.0) * ht)
235 *(-1776.0 * cW4
236 + (72.0 - 6250.0 * cW2 - 3056.0 * cW4 + 3696.0 * cW6) * ht
237 + (-18.0 + 1283.0 * cW2 + 1371.0 * cW4 - 1436.0 * cW6) * ht2
238 + (68.0 * cW2 - 124.0 * cW4 + 128.0 * cW6) * ht3)
239 + (6.0 * cW2 * ht - 37.0 * cW2 - 119.0 * ht2 + 56.0 * cW2 * ht2)
240 * M_PI * M_PI / 27.0 / ht2
241 + (32.0 * cW4 / 3.0 - 2.0 / 3.0 - 12.0 * cW2) * B0_Mt2_Mz2_Mw2_Mw2
242 + (20.0 / 3.0 + 1.0 / 3.0 / cW2 - 8.0 * cW2) * B0_Mt2_Mw2_Mw2_Mz2
243 + (17.0 - 58.0 * cW2 + 32.0 * cW4)*(4.0 - ht) * sqrt(ht) * g(ht) / 27.0
244 - 40.0 * sW2 * (4.0 - ht) * Lambda(ht) / 3.0 / ht
245 + 2.0 * cW2 * (37.0 - 6.0 * ht - 12.0 * ht2 - 22.0 * ht3 + 9.0 * ht4)
246 * Li21mht / 9.0 / ht2
247 - (1.0 + 14.0 * cW2 + 16.0 * cW4) * log(cW2) / 3.0
248 + (11520.0 - 15072.0 * cW2
249 - (7170.0 - 8928.0 * cW2 - 768.0 * cW4) * ht
250 + (3411.0 - 7062.0 * cW2 + 3264.0 * cW4) * ht2
251 - (1259.0 - 3547.0 * cW2 + 2144.0 * cW4) * ht3
252 + (238.0 - 758.0 * cW2 + 448.0 * cW4) * ht4
253 - (17.0 - 58.0 * cW2 + 32.0 * cW4) * ht5)
254 * log(ht) / 27.0 / (ht - 4.0) / (ht - 4.0) / ht
255 + 8.0 / 9.0 * (4.0 - 26.0 * cW2 - 5.0 * cW4) * log(Mt * Mt / mu / mu)
256 + (3.0 + 5.0 * cW2 - 26.0 * cW4 - 48.0 * cW6) * log(zt) / 9.0 / cW2
257 + (3840.0 * sW2 - (4310.0 - 4224.0 * cW2 - 256.0 * cW4) * ht
258 + (1706.0 - 1312.0 * cW2 - 320.0 * cW4) * ht2
259 - (315.0 + 476.0 * cW2 - 64.0 * cW4) * ht3
260 + (24.0 + 454.0 * cW2) * ht4 - 112.0 * cW2 * ht5
261 + 9.0 * cW2 * ht6)
262 * phi(ht / 4.0) / 9.0 / (ht - 4.0) / (ht - 4.0) / ht2));
263 }
264 return rho2;
265}
266
267double EWSMTwoLoopEW::DeltaRho2Add(const double Mw_i) const
268{
269 double Mz = cache.getSM().getMz(), Mz2 = Mz*Mz;
270 double Mw = Mw_i, Mw2 = Mw*Mw;
271 double cW2 = cache.getSM().cW2(Mw);
272 double cW4 = cW2*cW2, cW6 = cW4*cW2;
273 double Mt = cache.getSM().getMtpole(), Mt2 = Mt*Mt;
274 double zt = Mz2 / Mt2;
275 double mu = Mt; // renormalization scale
276
277 double B0_Mt2_Mz2_Mw2_Mw2 = cache.getPV().B0(Mt2, Mz2, Mw2, Mw2).real();
278 double B0_Mt2_Mw2_0_Mw2 = cache.getPV().B0(Mt2, Mw2, 0.0, Mw2).real();
279 double B0_Mt2_Mw2_Mw2_Mz2 = cache.getPV().B0(Mt2, Mw2, Mw2, Mz2).real();
280
281 double dRho2add = 542.0 / 27.0 - 2.0 / 3.0 / cW2 - 800.0 * cW2 / 27.0
282 + 1.0 / 3.0 * (1.0 + 26.0 * cW2 + 24.0 * cW4) * B0_Mt2_Mz2_Mw2_Mw2
283 + 4.0 * cW2 * B0_Mt2_Mw2_0_Mw2
284 - (11.0 / 3.0 + 1.0 / 3.0 / cW2 + 4.0 * cW2) * B0_Mt2_Mw2_Mw2_Mz2
285 - (2.0 / 3.0 + 4.0 * cW2 / 3.0 - 8.0 * cW4) * log(cW2)
286 + (1.0 / cW2 - 38.0 / 3.0 + 34.0 * cW2 / 3.0) * log(Mt2 / mu / mu)
287 + 2.0 * (3.0 - 62.0 * cW2 + 74.0 * cW4 + 36.0 * cW6) * log(zt) / 9.0 / cW2;
288 return dRho2add;
289}
290
291double EWSMTwoLoopEW::DeltaRw2(const double Mw_i) const
292{
293 double Mz = cache.getSM().getMz(), Mz2 = Mz*Mz;
294 double Mw = Mw_i, Mw2 = Mw*Mw;
295 double cW2 = cache.getSM().cW2(Mw), sW2 = cache.getSM().sW2(Mw);
296 double cW4 = cW2*cW2, cW6 = cW4*cW2;
297 double Mt = cache.getSM().getMtpole(), Mt2 = Mt*Mt;
298 double mh = cache.getSM().getMHl(), mh2 = mh*mh;
299 double zt = Mz2 / Mt2;
300 double zt2 = zt*zt;
301 double ht = mh2 / Mt2;
302 double ht2 = ht*ht, ht3 = ht2*ht, ht4 = ht3*ht, ht5 = ht4*ht;
303 double mu = Mt; // renormalization scale
304
305 double dRw2;
306 if (mh < 0.3 * Mt) {
307 double B0_Mt2_Mw2_mh2_Mw2 = cache.getPV().B0(Mt2, Mw2, mh2, Mw2).real();
308 double B0_Mt2_Mw2_Mw2_Mz2 = cache.getPV().B0(Mt2, Mw2, Mw2, Mz2).real();
309 dRw2 = -13.0 / 144.0 - 1.0 / 48.0 / cW4 - 41.0 / 96.0 / cW2 + 61.0 * cW2 / 72.0
310 + (7.0 - 16.0 * cW2) / 27.0 * M_PI * sqrt(ht) - M_PI * M_PI / 36.0
311 - 5.0 * ht2 / 144.0 / cW4 / zt2 + 35.0 * ht / 288.0 / cW2 / zt
312 + 5.0 / 12.0 * (1.0 + ht2 / 12.0 / cW4 / zt2 - ht / 3.0 / cW2 / zt) * B0_Mt2_Mw2_mh2_Mw2
313 + (1.0 + 20.0 * cW2 - 24.0 * cW4) / 48.0 / cW4 * B0_Mt2_Mw2_Mw2_Mz2
314 - (5.0 * sW2 * ht2 + 3.0 * ht * zt + 48.0 * cW2 * ht * zt - 60.0 * cW4 * ht * zt
315 - 3.0 * cW2 * zt2 - 8.0 * cW4 * zt2 + 20.0 * cW6 * zt2) * log(cW2)
316 / 144.0 / cW2 / sW2 / zt / (ht - cW2 * zt)
317 + 5.0 * ht * (ht2 - 4.0 * cW2 * ht * zt + 12.0 * cW4 * zt2) * log(ht)
318 / 144.0 / cW4 / zt2 / (ht - cW2 * zt)
319 + (17.0 / 36.0 - 13.0 * cW2 / 18.0) * log(Mt2 / mu / mu)
320 - (5.0 * cW2 * ht2 - 3.0 * ht * zt - 60.0 * cW2 * ht * zt + 60.0 * cW4 * ht * zt
321 + (3.0 * cW2 + 60.0 * cW4 - 20.0 * cW6) * zt2) * log(zt)
322 / 144.0 / cW4 / zt / (ht - cW2 * zt);
323 } else {
324 double B0_Mt2_Mw2_Mw2_Mz2 = cache.getPV().B0(Mt2, Mw2, Mw2, Mz2).real();
325 dRw2 = -121.0 / 288.0 - 1.0 / 48.0 / cW4 - 41.0 / 96.0 / cW2 + 77.0 * cW2 / 12.0
326 + 19.0 / 72.0 / ht + (41.0 / 216.0 - 4.0 * cW2 / 27.0) * ht
327 - (19.0 + 21.0 * ht) * M_PI * M_PI / 432.0 / ht2
328 - (1.0 / 2.0 - 1.0 / 48.0 / cW4 - 5.0 / 12.0 / cW2) * B0_Mt2_Mw2_Mw2_Mz2
329 + (16.0 * cW2 - 7.0) / 216.0 * (ht - 4.0) * sqrt(ht) * g(ht)
330 - (1.0 / 12.0 - 1.0 / 3.0 / ht) * Lambda(ht)
331 + (19.0 + 21.0 * ht - 12.0 * ht2 - 31.0 * ht3 + 9.0 * ht4) / 72.0 / ht2
332 * cache.getPolyLog().Li2(1.0 - ht).real()
333 - (1.0 + 21.0 * cW2 - 25.0 * cW4) * log(cW2) / 48.0 / cW2 / sW2
334 + (17.0 / 36.0 - 13.0 * cW2 / 18.0) * log(Mt2 / mu / mu)
335 + (1.0 + 20.0 * cW2 - 25.0 * cW4) * log(zt) / 48.0 / cW4
336 + (372.0 + (96.0 * cW2 - 213.0) * ht + (432.0 * cW2 - 318.0) * ht2
337 + (97.0 - 160.0 * cW2) * ht3 - (7.0 - 16.0 * cW2) * ht4)
338 / 216.0 / (ht - 4.0) / ht * log(ht)
339 + (96.0 - (384.0 - 64.0 * cW2) * ht - (2.0 + 64.0 * cW2) * ht2 + 231.0 * ht3
340 - 85.0 * ht4 + 9.0 * ht5) / 144.0 / (ht - 4.0) / ht2 * phi(ht / 4.0);
341 }
342 return dRw2;
343}
344
346{
347 double Mt = cache.getSM().getMtpole(), Mt2 = Mt*Mt;
348 double mh = cache.getSM().getMHl(), mh2 = mh*mh;
349 double ht = mh2 / Mt2;
350 double ht2 = ht*ht, ht3 = ht2*ht;
351 double mu = Mt; // renormalization scale
352
353 double dEoE2;
354 if (mh < 0.3 * Mt) {
355 dEoE2 = 61.0 / 72.0 - 16.0 * sqrt(ht) * M_PI / 27.0 - 13.0 / 18.0 * log(Mt2 / mu / mu);
356 } else {
357 dEoE2 = (231.0 - 32.0 * ht) / 216.0 - 2.0 / 27.0 * (4.0 - ht) * sqrt(ht) * g(ht)
358 + 2.0 * (6.0 + 27.0 * ht - 10.0 * ht2 + ht3) / 27.0 / (ht - 4.0) * log(ht)
359 - 13.0 / 18.0 * log(Mt2 / mu / mu)
360 - 4.0 * (ht - 1.0) / 9.0 / (ht - 4.0) / ht * phi(ht / 4.0);
361 }
362 return dEoE2;
363}
364
365double EWSMTwoLoopEW::f2Add(const double Mw_i) const
366{
367 double Mz = cache.getSM().getMz(), Mz2 = Mz*Mz;
368 double Mw = Mw_i, Mw2 = Mw*Mw;
369 double cW2 = cache.getSM().cW2(Mw), sW2 = cache.getSM().sW2(Mw);
370 double Mt = cache.getSM().getMtpole(), Mt2 = Mt*Mt;
371 double zt = Mz2 / Mt2;
372
373 double B0_Mt2_Mw2_0_Mw2 = cache.getPV().B0(Mt2, Mw2, 0.0, Mw2).real();
374 double B0_Mt2_Mw2_Mw2_Mz2 = cache.getPV().B0(Mt2, Mw2, Mw2, Mz2).real();
375
376 double f2a = 10.0 / 3.0 + 1.0 / 3.0 / cW2 + 4.0 * cW2 * B0_Mt2_Mw2_0_Mw2
377 - (11.0 / 3.0 + 1.0 / 3.0 / cW2 + 4.0 * cW2) * B0_Mt2_Mw2_Mw2_Mz2
378 + (11.0 - 8.0 * cW2) * log(cW2) / 6.0 / sW2
379 - (11.0 / 3.0 + 1.0 / 3.0 / cW2) * log(zt);
380 return f2a;
381}
382
383double EWSMTwoLoopEW::DeltaEta2(const double Mw_i) const
384{
385 double Mz = cache.getSM().getMz(), Mz2 = Mz*Mz;
386 double Mw = Mw_i, Mw2 = Mw*Mw;
387 double cW2 = cache.getSM().cW2(Mw);
388 double cW4 = cW2*cW2, cW6 = cW4*cW2;
389 double Mt = cache.getSM().getMtpole(), Mt2 = Mt*Mt;
390 double mh = cache.getSM().getMHl(), mh2 = mh*mh;
391 double zt = Mz2 / Mt2;
392 double zt2 = zt*zt, zt3 = zt2*zt;
393 double ht = mh2 / Mt2;
394 double ht2 = ht*ht, ht3 = ht2*ht, ht4 = ht3*ht, ht5 = ht4*ht;
395 double mu = Mt; // renormalization scale
396
397 double dEta2;
398 if (mh < 0.57 * Mt) {
399 double B0_Mt2_Mz2_Mw2_Mw2 = cache.getPV().B0(Mt2, Mz2, Mw2, Mw2).real();
400 double B0_Mt2_Mz2_mh2_Mz2 = cache.getPV().B0(Mt2, Mz2, mh2, Mz2).real();
401 dEta2 = (ht3 - 6.0 * ht2 * zt + 11.0 * ht * zt2) / 9.0 / cW2 / (ht - 4.0 * zt) / zt2
402 + (49.0 - 289.0 * cW2 - 349.0 * cW4 + 292.0 * cW6) / 216.0 / cW2 / (1.0 - 4.0 * cW2)
403 + (1.0 + 18.0 * cW2 - 16.0 * cW4) / 12.0 / (1.0 - 4.0 * cW2) * log(cW2)
404 - (17.0 - 40.0 * cW2 + 32.0 * cW4) / 54.0 / cW2 * (sqrt(ht) * M_PI - 2.0)
405 + (11.0 * ht2 * zt - 2.0 * ht3 - 24.0 * ht * zt2 + 24.0 * zt3)
406 / 18.0 / cW2 / (ht - 4.0 * zt) / zt2 * log(ht)
407 + (1.0 - 4.0 * cW2 + 44.0 * cW4 - 32.0 * cW6) / 24.0 / cW2 / (1.0 - 4.0 * cW2)
408 * B0_Mt2_Mz2_Mw2_Mw2
409 + (13.0 * ht2 * zt - 2.0 * ht3 - 32.0 * ht * zt2 + 36.0 * zt3)
410 / 18.0 / cW2 / (ht - 4.0 * zt) / zt2 * B0_Mt2_Mz2_mh2_Mz2
411 - (17.0 - 34.0 * cW2 + 26.0 * cW4) / 36.0 / cW2 * log(Mt2 / mu / mu)
412 + (ht * (2.0 * ht - 5.0 * zt) / 18.0 / cW2 / zt / (ht - 4.0 * zt)
413 + (10.0 - 39.0 * cW2 - 70.0 * cW4 + 48.0 * cW6) / 36.0 / cW2 / (4.0 * cW2 - 1.0))
414 * log(zt);
415 } else {
416 double B0_Mt2_Mz2_Mw2_Mw2 = cache.getPV().B0(Mt2, Mz2, Mw2, Mw2).real();
417 dEta2 = (-17.0 + 40.0 * cW2 - 32.0 * cW4) * ht / 216.0 / cW2
418 + 5.0 / 144.0 / cW2 / (ht - 4.0)
419 + (707.0 - 4720.0 * cW2 + 5900.0 * cW4 - 3696.0 * cW6) / 864.0 / cW2 / (1.0 - 4.0 * cW2)
420 + (10.0 / 27.0 - 17.0 / 108.0 / cW2 - 8.0 * cW2 / 27.0)*(1.0 - ht / 4.0) * sqrt(ht) * g(ht)
421 + (1.0 + 18.0 * cW2 - 16.0 * cW4) / 12.0 / (1.0 - 4.0 * cW2) * log(cW2)
422 + (4.0 - ht) / 12.0 / cW2 / ht * Lambda(ht)
423 + (2.0 - 7.0 * cW2 - 70.0 * cW4 + 48.0 * cW6) / 36.0 / cW2 / (4.0 * cW2 - 1.0) * log(zt)
424 + (1.0 - 4.0 * cW2 + 44.0 * cW4 - 32.0 * cW6) / 24.0 / cW2 / (1.0 - 4.0 * cW2) * B0_Mt2_Mz2_Mw2_Mw2
425 - (17.0 - 34.0 * cW2 + 26.0 * cW4) / 36.0 / cW2 * log(Mt2 / mu / mu)
426 + ((4.0 * cW2 - 5.0)*(6.0 + 27.0 * ht - 10.0 * ht2 + ht3) / 54.0 / (ht - 4.0)
427 - (1152.0 + 606.0 * ht + 1467.0 * ht2 - 1097.0 * ht3 + 238.0 * ht4 - 17.0 * ht5)
428 / 432.0 / cW2 / (ht - 4.0) / (ht - 4.0) / ht) * log(ht)
429 + ((5.0 - 4.0 * cW2)*(ht - 1.0) / 9.0 / (ht - 4.0) / ht
430 - (384.0 + 10.0 * ht - 238.0 * ht2 + 63.0 * ht3 - 3.0 * ht4)
431 / 144.0 / cW2 / (ht - 4.0) / (ht - 4.0) / ht2) * phi(ht / 4.0);
432 }
433 return dEta2;
434}
435
436gslpp::complex EWSMTwoLoopEW::DeltaEta2Add_tmp(const double I3f, const double Qf,
437 const double Mw_i) const
438{
439 double Mz = cache.getSM().getMz(), Mz2 = Mz*Mz;
440 double Mw = Mw_i, Mw2 = Mw*Mw;
441 double cW2 = cache.getSM().cW2(Mw);
442 double cW4 = cW2*cW2, cW6 = cW4*cW2;
443 double Mt = cache.getSM().getMtpole(), Mt2 = Mt*Mt;
444 double zt = Mz2 / Mt2;
445 double mu = Mt; // renormalization scale
446
447 double B0_Mt2_Mz2_Mw2_Mw2 = cache.getPV().B0(Mt2, Mz2, Mw2, Mw2).real();
448
449 gslpp::complex dEta2add = 16.0 * M_PI * M_PI * DeltaEtaf1(I3f, Qf, Mw_i) + Vadd(I3f, Qf, Mw_i)
450 - (197.0 - 1378.0 * cW2 + 1064.0 * cW4) / 27.0 / (1.0 - 4.0 * cW2)
451 - (1.0 + 16.0 * cW2 - 20.0 * cW4 + 48.0 * cW6) / 3.0 / (1.0 - 4.0 * cW2)
452 * B0_Mt2_Mz2_Mw2_Mw2
453 - 2.0 * cW2 * (1.0 + 26.0 * cW2 + 24.0 * cW4) / 3.0 / (1.0 - 4.0 * cW2)
454 * log(cW2)
455 + (41.0 / 3.0 - 46.0 * cW2 / 3.0) * log(Mt2 / mu / mu)
456 + 2.0 * (50.0 - 283.0 * cW2 + 242.0 * cW4 - 72.0 * cW6)
457 / 9.0 / (1.0 - 4.0 * cW2) * log(zt);
458 return dEta2add;
459}
460
461gslpp::complex EWSMTwoLoopEW::DeltaEta2Add_f(const Particle f, const double Mw_i) const
462{
463 double I3f = cache.I3_f(f);
464 double Qf = cache.Q_f(f);
465 return DeltaEta2Add_tmp(I3f, Qf, Mw_i);
466}
467
468double EWSMTwoLoopEW::DeltaKappa2(const double Mw_i) const
469{
470 double Mz = cache.getSM().getMz(), Mz2 = Mz*Mz;
471 double Mw = Mw_i, Mw2 = Mw*Mw;
472 double cW2 = cache.getSM().cW2(Mw), sW2 = cache.getSM().sW2(Mw);
473 double Mt = cache.getSM().getMtpole(), Mt2 = Mt*Mt;
474 double mh = cache.getSM().getMHl(), mh2 = mh*mh;
475 double zt = Mz2 / Mt2;
476 double ht = mh2 / Mt2;
477 double ht2 = ht*ht, ht3 = ht2*ht;
478 double mu = Mt; // renormalization scale
479
480 double B0_Mt2_Mz2_Mw2_Mw2 = cache.getPV().B0(Mt2, Mz2, Mw2, Mw2).real();
481
482 double dKappa2;
483 if (mh < 0.57 * Mt) {
484 dKappa2 = (-175.0 + 366.0 * sW2) / 432.0
485 + (3.0 / 8.0 - sW2 / 3.0) * B0_Mt2_Mz2_Mw2_Mw2 - cW2 / 6.0 * log(cW2)
486 - 2.0 * M_PI / 27.0 * sqrt(ht)*(8.0 * sW2 - 3.0)
487 - (1.0 / 4.0 + 2.0 / 9.0 * sW2) * log(Mt2 / mu / mu)
488 + (3.0 * sW2 - 2.0) / 18.0 * log(zt);
489 } else {
490 dKappa2 = (-211.0 + 24.0 * ht + 462.0 * sW2 - 64.0 * ht * sW2) / 432.0
491 + (3.0 / 8.0 - sW2 / 3.0) * B0_Mt2_Mz2_Mw2_Mw2
492 - cW2 / 6.0 * log(cW2)
493 + (ht - 4.0) * sqrt(ht)*(8.0 * sW2 - 3.0) * g(ht) / 108.0
494 - (6.0 + 27.0 * ht - 10.0 * ht2 + ht3)*(3.0 - 8.0 * sW2)
495 / 108.0 / (ht - 4.0) * log(ht)
496 - (1.0 / 4.0 + 2.0 / 9.0 * sW2) * log(Mt2 / mu / mu)
497 + (3.0 * sW2 - 2.0) / 18.0 * log(zt)
498 + (ht - 1.0)*(8.0 * sW2 - 3.0) / 18.0 / (4.0 - ht) / ht * phi(ht / 4.0);
499 }
500 return dKappa2;
501}
502
503gslpp::complex EWSMTwoLoopEW::DeltaKappa2Add_tmp(const double I3f, const double Qf,
504 const double Mw_i) const
505{
506 double Mz = cache.getSM().getMz(), Mz2 = Mz*Mz;
507 double Mw = Mw_i;
508 double cW2 = cache.getSM().cW2(Mw);
509 double cW4 = cW2*cW2;
510 double Mt = cache.getSM().getMtpole(), Mt2 = Mt*Mt;
511 double zt = Mz2 / Mt2;
512 double mu = Mt; // renormalization scale
513
514 gslpp::complex i = gslpp::complex::i();
515
516 gslpp::complex dKappa2add = -238.0 * cW2 / 27.0 + 8.0 * cW4
517 - 2.0 * cW2 * sqrt(4.0 * cW2 - 1.0)*(3.0 + 4.0 * cW2)
518 * atan(1.0 / sqrt(4.0 * cW2 - 1.0))
519 - 16.0 / 9.0 * cW2 * log(zt)
520 + (1.0 - 12.0 * I3f * Qf + 8.0 * Qf * Qf - 8.0 * cW4 * Qf * Qf)
521 / 4.0 / cW2 * FV(1)
522 + 4.0 * cW2 * GV(1.0 / cW2) - 7.0 * cW2 * log(cW2)
523 - 17.0 / 3.0 * cW2 * log(mu * mu / Mz2)
524 + cW2 * (1.0 - 2.0 * Qf * I3f) * FV(1.0 / cW2)
525 - i * 80.0 / 9.0 * M_PI;
526 return dKappa2add;
527}
528
529gslpp::complex EWSMTwoLoopEW::DeltaKappa2Add_f(const Particle f, const double Mw_i) const
530{
531 double I3f = cache.I3_f(f);
532 double Qf = cache.Q_f(f);
533 return DeltaKappa2Add_tmp(I3f, Qf, Mw_i);
534}
535
536gslpp::complex EWSMTwoLoopEW::Vadd(const double I3f, const double Qf,
537 const double Mw_i) const
538{
539 double Mw = Mw_i, Mw2 = Mw*Mw;
540 double cW2 = cache.getSM().cW2(Mw), sW2 = cache.getSM().sW2(Mw);
541 double Mt = cache.getSM().getMtpole();
542 double mu = Mt; // renormalization scale
543
544 gslpp::complex V = 8.0 * cW2 * log(Mw2 / mu / mu)
545 + 3.0 * (2.0 * I3f * Qf - 4.0 * sW2 * Qf * Qf) * FV(1.0)
546 - 16.0 * cW2 * GV(1.0 / cW2)
547 + (1.0 - 4.0 * cW2 - 4.0 * (1.0 - 2.0 * cW2) * I3f * Qf) * FV(1.0 / cW2);
548 return V;
549}
550
551gslpp::complex EWSMTwoLoopEW::DeltaEtaf1(const double I3f, const double Qf,
552 const double Mw_i) const
553{
554 double Mz = cache.getSM().getMz(), Mz2 = Mz*Mz;
555 double Mw = Mw_i;
556 double cW2 = cache.getSM().cW2(Mw);
557
558 gslpp::complex SigmaPrime_ZZ = myOneLoopEW.SigmabarPrime_ZZ_bos_Mz2(Mz, Mw)
560
561 gslpp::complex dEtaf1 = 1.0 / 16.0 / M_PI / M_PI
562 * (-SigmaPrime_ZZ / cW2 - 4.0 * cW2 * log(cW2) + Vfi(I3f, Qf, Mz2, Mw));
563 return dEtaf1;
564}
565
566gslpp::complex EWSMTwoLoopEW::Vfi(const double I3f, const double Qf,
567 const double q2, const double Mw_i) const
568{
571 double I3aQaQa = I3i * Qi * Qi + I3f * Qf*Qf;
572 double I3aQa = I3i * Qi + I3f*Qf;
573 double QaQa = Qi * Qi + Qf*Qf;
574 double Mz = cache.getSM().getMz(), Mz2 = Mz*Mz;
575 double Mw = Mw_i, Mw2 = Mw*Mw;
576 double cW2 = cache.getSM().cW2(Mw), sW2 = cache.getSM().sW2(Mw);
577
578 gslpp::complex V = (1.0 - sW2 * (2.0 - 2.0 * I3aQaQa)) * FV(q2 / Mw2)
579 + 8.0 * cW2 * GV(q2 / Mw2)
580 - (1.0 - 6.0 * sW2 * I3aQa + 6.0 * sW2 * QaQa) / 2.0 / cW2 * FV(q2 / Mz2);
581 return V;
582}
583
584double EWSMTwoLoopEW::Lambda(const double x) const
585{
586 if (x >= 0.0 && x <= 4.0) {
587 return ( -1.0 / 2.0 / sqrt(x) * g(x) + M_PI / 2.0 * sqrt(4.0 / x - 1.0));
588 } else if (x > 4.0) {
589 return ( -1.0 / 2.0 / sqrt(x) * g(x));
590 } else
591 throw std::runtime_error("Out of range in EWSMTwoLoopEW::Lambda()");
592}
593
594double EWSMTwoLoopEW::phi(const double x) const
595{
596 if (x >= 0.0 && x <= 1.0) {
597 return ( 4.0 * sqrt(x / (1.0 - x)) * cache.getClausen().Cl2(2.0 * asin(sqrt(x))));
598 } else if (x > 1.0) {
599 double lambda = sqrt(1.0 - 1.0 / x);
600 return ( 1.0 / lambda * (-4.0 * cache.getPolyLog().Li2((1.0 - lambda) / 2.0).real()
601 + 2.0 * pow(log((1.0 - lambda) / 2.0), 2.0)
602 - pow(log(4.0 * x), 2.0) + M_PI * M_PI / 3.0));
603 } else
604 throw std::runtime_error("Out of range in EWSMTwoLoopEW::phi()");
605}
606
607gslpp::complex EWSMTwoLoopEW::FV(const double x) const
608{
609 if (x <= 0.0)
610 throw std::runtime_error("Out of range in EWSMTwoLoopEW::FV()");
611
612 gslpp::complex i = gslpp::complex::i();
613
614 return ( i * M_PI * (2.0 / x + 3.0 - 2.0 * (1.0 + 1.0 / x)*(1.0 + 1.0 / x) * log(1.0 + x))
615 + 2.0 / x + 7.0 / 2.0 - (3.0 + 2.0 / x) * log(x)
616 + (1.0 + 1.0 / x)*(1.0 + 1.0 / x)
617 *(2.0 * cache.getPolyLog().Li2(1.0 / (1.0 + x))
618 - M_PI * M_PI / 3.0 + pow(log(1.0 + x), 2.0)));
619}
620
621gslpp::complex EWSMTwoLoopEW::GV(const double x) const
622{
623 if (x <= 0.0 || x >= 4.0)
624 throw std::runtime_error("Out of range in EWSMTwoLoopEW::GV()");
625
626 double atanX = atan(sqrt(x / (4.0 - x)));
627
628 return ( (sqrt((4.0 - x) / x) * atanX - 1.0)*(1.0 / x + 1.0 / 2.0) + 9.0 / 8.0
629 + 1.0 / 2.0 / x - (1.0 + 1.0 / 2.0 / x)*4.0 / x * atanX * atanX);
630}
631
632
633
634
double Cl2(const double phi) const
The Clausen function of index 2, .
gslpp::complex SigmabarPrime_ZZ_bos_Mz2(const double mu, const double Mw_i) const
The derivative of the bosonic contribution to the self-energy of the boson for in the Unitary gauge...
gslpp::complex PibarZgamma_fer(const double mu, const double s, const double Mw_i) const
The fermionic contribution to the self-energy of the - mixing in the Unitary gauge,...
gslpp::complex SigmabarPrime_ZZ_fer_Mz2(const double mu, const double Mw_i) const
The derivative of the fermionic contribution to the self-energy of the boson for in the Unitary gau...
gslpp::complex FV(const double x) const
The auxiliary function .
gslpp::complex GV(const double x) const
The auxiliary function .
double deltaEoverE2() const
The auxiliary function .
double DeltaKappa2(const double Mw_i) const
The auxiliary function .
gslpp::complex deltaRho_rem_f(const Particle f, const double Mw_i) const
Remainder contribution of to the effective couplings , denoted as .
double DeltaRho2Add(const double Mw_i) const
The auxiliary function .
gslpp::complex DeltaEtaf1(const double I3f, const double Qf, const double Mw_i) const
The auxiliary function .
gslpp::complex DeltaEta2Add_tmp(const double I3f, const double Qf, const double Mw_i) const
The auxiliary function .
double f2Add(const double Mw_i) const
The auxiliary function .
gslpp::complex Vfi(const double I3f, const double Qf, const double q2, const double Mw_i) const
The auxiliary function .
double DeltaRho(const double Mw_i) const
Leading two-loop contribution of to , denoted as .
double f1(const double a) const
The auxiliary function .
double g(const double a) const
The auxiliary function .
double DeltaRw2(const double Mw_i) const
The auxiliary function .
double DeltaAlpha_t(const double s) const
Top-quark contribution of to the electromagnetic coupling , denoted as .
double f0(const double a) const
The auxiliary function .
gslpp::complex DeltaEta2Add_f(const Particle f, const double Mw_i) const
The auxiliary function for .
const EWSMOneLoopEW myOneLoopEW
An object of type EWSMOneLoopEW.
double tau_2() const
The function .
double phi(const double x) const
The auxiliary function .
gslpp::complex DeltaKappa2Add_tmp(const double I3f, const double Qf, const double Mw_i) const
The auxiliary function .
double Lambda(const double x) const
The auxiliary function .
EWSMTwoLoopEW(const EWSMcache &cache_i)
Constructor.
double DeltaAlpha_l(const double s) const
Leptonic contribution of to the electromagnetic coupling , denoted as .
double DeltaRho2(const double Mw_i) const
The auxiliary function .
gslpp::complex Vadd(const double I3f, const double Qf, const double Mw_i) const
The auxiliary function .
double rho_2() const
The function .
double DeltaR_rem(const double Mw_i) const
Remainder contribution of to , denoted as .
double DeltaEta2(const double Mw_i) const
The auxiliary function .
const EWSMcache & cache
A reference to an object of type EWSMcache.
gslpp::complex DeltaKappa2Add_f(const Particle f, const double Mw_i) const
The auxiliary function for .
gslpp::complex deltaKappa_rem_f(const Particle f, const double Mw_i) const
Remainder contribution of to the effective couplings , denoted as .
A class for cache variables used in computing radiative corrections to the EW precision observables.
Definition: EWSMcache.h:40
double Xt_alpha(const double Mw_i) const
The quantity with the coupling .
Definition: EWSMcache.h:355
double logMTOPtoMH() const
A cache method.
Definition: EWSMcache.cpp:134
const ClausenFunctions getClausen() const
A get method to access the member object of type ClausenFunctions.
Definition: EWSMcache.h:124
const Polylogarithms getPolyLog() const
A get method to access the member object of type Polylogarithms.
Definition: EWSMcache.h:115
double logMZtoMTAU() const
A cache method.
Definition: EWSMcache.cpp:106
double logMZtoME() const
A cache method.
Definition: EWSMcache.cpp:78
double logMZtoMMU() const
A cache method.
Definition: EWSMcache.cpp:92
const PVfunctions getPV() const
A get method to access the member reference to the object of type StandardModel passed to the constru...
Definition: EWSMcache.h:106
double getZeta3() const
A get method to access the value of the zeta function .
Definition: EWSMcache.h:146
double Q_f(const Particle f) const
The charge of an SM fermion .
Definition: EWSMcache.h:268
double mf(const Particle f, const double mu=0.0, const orders order=FULLNNLO) const
The mass of an SM fermion.
Definition: EWSMcache.cpp:49
double I3_f(const Particle f) const
The isospin of an SM fermion .
Definition: EWSMcache.h:278
const StandardModel & getSM() const
Definition: EWSMcache.h:56
An observable class for the -boson mass.
Definition: Mw.h:22
gslpp::complex B0(const double mu2, const double p2, const double m02, const double m12) const
.
Definition: PVfunctions.cpp:41
A class for particles.
Definition: Particle.h:26
bool is(std::string name_i) const
Definition: Particle.cpp:23
gslpp::complex Li2(const double x) const
The dilogarithm with a real argument, .
const double getMtpole() const
A get method to access the pole mass of the top quark.
Definition: QCD.h:600
@ MU
Definition: QCD.h:314
@ ELECTRON
Definition: QCD.h:312
@ TAU
Definition: QCD.h:316
const Particle & getLeptons(const QCD::lepton p) const
A get method to retrieve the member object of a lepton.
const double getMz() const
A get method to access the mass of the boson .
virtual const double getMHl() const
A get method to retrieve the Higgs mass .
virtual const double cW2(const double Mw_i) const
The square of the cosine of the weak mixing angle in the on-shell scheme, denoted as .
virtual const double sW2(const double Mw_i) const
The square of the sine of the weak mixing angle in the on-shell scheme, denoted as .
const double getAle() const
A get method to retrieve the fine-structure constant .
Test Observable.
A class for , relevant for mesons mixing in the Standard Model.
Definition: xi.h:23