a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
StandardModelMatching.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 "StandardModelMatching.h"
9#include "StandardModel.h"
10#include "QCD.h"
11#include <stdexcept>
12#include <sstream>
13
14#define ZEW_NUMERIC
15
16StandardModelMatching::StandardModelMatching(const StandardModel & SM_i)
17: SM(SM_i),
18mcdbd2(5, NDR, NLO),
19mcdbs2(5, NDR, NLO),
20mcdd2(5, NDR, NLO),
21mcdk2(5, NDR, NLO),
22mck(10, NDR, NLO, NLO_QED11),
23mckcc(10, NDR, NLO, NLO_QED11),
24mcbsg(8, NDR, NNLO),
25mcprimebsg(8, NDR, NNLO),
26mcBMll(13, NDR, NLO),
27mcprimeBMll(13, NDR, NLO),
28mcbnlep(10, NDR, NLO, NLO_QED11),
29mcbnlepCC(10, NDR, NLO),
30mcd1(10, NDR, NLO),
31mcd1Buras(10, NDR, NLO),
32mckpnn(2, NDR, NLO, NLO_QED11),
33mckmm(1, NDR, NLO),
34mcbsnn(2, NDR, NLO, NLO_QED11),
35mcbdnn(2, NDR, NLO, NLO_QED11),
36mcbsmm(8, NDR, NNLO, NLO_QED22),
37mcbdmm(8, NDR, NNLO, NLO_QED22),
38mculeptonnu(5, NDR, LO),
39mcDLij(2, NDR, LO),
40mcDLi3j(20, NDR, LO),
41mcmueconv(8, NDR, LO),
42mcgminus2mu(2, NDR, NLO),
43mcC(2, NDR, QCD2, QED2), // blocks have to have the same orders
44mcP(4, NDR, QCD2, QED2),
45mcM(2, NDR, QCD2, QED2),
46mcL(2, NDR, QCD2, QED2),
47mcQ(4, NDR, QCD2, QED2),
48mcB(1, NDR, QCD2, QED2),
49Vckm(3, 3, 0.) {
50 swa = 0.;
51 swb = 0.;
52 swc = 0.;
53 xcachea = 0.;
54 xcacheb = 0.;
55 xcachec = 0.;
56
57
58 for (int j = 0; j < 10; j++) {
59 CWD1ArrayLO[j] = 0.;
60 CWD1ArrayNLO[j] = 0.;
61 CWbnlepArrayLOqcd[j] = 0.;
62 CWbnlepArrayNLOqcd[j] = 0.;
63 CWbnlepArrayLOew[j] = 0.;
64 CWbnlepArrayNLOew[j] = 0.;
65 };
66
67
68 for (int j = 0; j < 19; j++) {
69 CWBMllArrayLO[j] = 0.;
70 CWBMllArrayNLO[j] = 0.;
71 }
72
73 for (int j = 0; j < 8; j++) {
74 CWbsgArrayLO[j] = 0.;
75 CWbsgArrayNLO[j] = 0.;
76 CWbsgArrayNNLO[j] = 0.;
77 CWprimebsgArrayLO[j] = 0.;
78 CWprimebsgArrayNLO[j] = 0.;
79 }
80
81 for (int j = 0; j < 8; j++) {
82 CWBsmmArrayNNLOqcd[j] = 0.;
83 CWBsmmArrayNLOqcd[j] = 0.;
84 CWBsmmArrayLOqcd[j] = 0.;
85 CWBsmmArrayNLOewt4[j] = 0.;
86 CWBsmmArrayNLOewt2[j] = 0.;
87 CWBsmmArrayNLOew[j] = 0.;
88 }
89
90 for (int j = 0; j < 8; j++) {
91 CWBdmmArrayNNLOqcd[j] = 0.;
92 CWBdmmArrayNLOqcd[j] = 0.;
93 CWBdmmArrayLOqcd[j] = 0.;
94 CWBdmmArrayNLOewt4[j] = 0.;
95 CWBdmmArrayNLOewt2[j] = 0.;
96 CWBdmmArrayNLOew[j] = 0.;
97 }
98
99
100 Nc = SM.getNc();
101 CF = SM.getCF();
102 gamma0 = 6. * (Nc - 1.) / Nc;
103 J5 = SM.Beta1(5) * gamma0 / 2. / SM.Beta0(5) / SM.Beta0(5) - ((Nc - 1.) / (2. * Nc) * (-21. + 57. / Nc - 19. / 3. * Nc + 20. / 3.)) / 2. / SM.Beta0(5);
104 BtNDR = 5. * (Nc - 1.) / 2. / Nc + 3. * CF;
105}
106
107StandardModelMatching::~StandardModelMatching() {
108}
109
110void StandardModelMatching::updateSMParameters() {
111 Mut = SM.getMut();
112 Muw = SM.getMuw();
113 Ale = SM.getAle();
114 Mt_muw = SM.Mrun(Muw, SM.getQuarks(QCD::TOP).getMass_scale(),
115 SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, FULLNNLO);
116 Mt_mut = SM.Mrun(Mut, SM.getQuarks(QCD::TOP).getMass_scale(),
117 SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, FULLNNLO);
118 alstilde = SM.Als(Muw, FULLNNNLO, true, true) / 4. / M_PI; // SM.Alstilde5(Muw) WHICH ONE TO USE?
119 aletilde = SM.Ale(Muw, FULLNLO) / 4. / M_PI; // Ale / 4. / M_PI; // WHERE IS ale(mu)?
120 GF = SM.getGF();
121 Mw_tree = SM.Mw_tree();
122 Mw = SM.Mw(); /* on-shell Mw */
123 sW2 = SM.sW2(); /* on-shell sW2 */
124 /* NP models should be added here after writing codes for Mw. */
125 /* The following is commented out till further review.*/
126 // if (SM.getModelName()=="StandardModel") {
127 // Mw = SM.Mw(); /* on-shell Mw */
128 // sW2 = SM.sW2(); /* on-shell sW2 */
129 // } else {
130 // Mw = Mw_tree;
131 // sW2 = 1.0 - Mw*Mw/SM.getMz()/SM.getMz();
132 // }
133 // sW2 = (M_PI * Ale ) / ( sqrt(2.) * GF * Mw * Mw ); // WARNING: only for
134 Vckm = SM.getVCKM();
135 lam_t = SM.getCKM().computelamt();
136 L = 2. * log(Muw / Mw);
137 Lz = 2. * log(Muw / SM.getMz());
138 mu_b = SM.getMub();
139}
140
141const double StandardModelMatching::x_c(const double mu, const orders order) const {
142 double mc = SM.Mrun(mu, SM.getQuarks(QCD::CHARM).getMass_scale(),
143 SM.getQuarks(QCD::CHARM).getMass(), QCD::CHARM, order);
144 return mc * mc / Mw / Mw;
145}
146
147const double StandardModelMatching::x_t(const double mu, const orders order) const {
148 double mt;
149
150 if (mu == Mut)
151 mt = Mt_mut;
152 else if (mu == Muw)
153 mt = Mt_muw;
154 else
155 mt = SM.Mrun(mu, SM.getQuarks(QCD::TOP).getMass_scale(),
156 SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, order);
157
158 //msbar mass here?
159 return mt * mt / Mw / Mw;
160}
161
162const double StandardModelMatching::mt2omh2(const double mu, const orders order) const {
163 double mt = SM.Mrun(mu, SM.getQuarks(QCD::TOP).getMass_scale(),
164 SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, order);
165 return (mt / SM.getMHl())*(mt / SM.getMHl());
166}
167
168const double StandardModelMatching::S0(double x) const {
169 return S0(x, x);
170}
171
172const double StandardModelMatching::S0(double x, double y) const
173{ // Buras 2000 (hep-ph/0007313v1) Appendix
174 if (fabs(1. - y / x) < LEPS){
175 return ((x * (-4. + 15. * x - 12. * x * x + x * x * x +
176 6. * x * x * log(x))) / (4. * pow(-1. + x, 3.)));
177 } else
178 return (x * y * ((1. / 4. + 3. / 2. / (1. - x) - 3. / 4. / pow(1. - x, 2.)) *
179 log(x) / (x - y) +
180 (1. / 4. + 3. / 2. / (1. - y) - 3. / 4. / pow(1. - y, 2.)) *
181 log(y) / (y - x) -
182 3. / 4. / (1. - x) / (1. - y)));
183}
184
185const double StandardModelMatching::S0p(double x) const {
186 double x2 = x * x;
187 return (x * (4. - 22. * x + 15. * x2 + 2. * x2 * x + x2 * x2 - 18. * x2 * log(x)) / 4. / pow(x - 1., 4.));
188}
189
190const double StandardModelMatching::S11(double x) const {
191 double x2 = x * x;
192 double x3 = x2 * x;
193 double x4 = x3 * x;
194 double xm3 = (x - 1.) * (x - 1.) * (x - 1.);
195 double xm4 = xm3 * (x - 1);
196
197 return (x * (4. - 39. * x + 168. * x2 + 11. * x3) / 4. / xm3
198 + 3. * x3 * gslpp_special_functions::dilog(1. - x) * (5. + x) / xm3
199 + 3. * x * log(x)*(-4. + 24. * x - 36. * x2 - 7. * x3 - x4) / 2.
200 / xm4 + 3. * x3 * pow(log(x), 2.) * (13. + 4. * x + x2) / 2.
201 / xm4);
202}
203
204const double StandardModelMatching::S18(double x) const {
205 double x2 = x * x;
206 double x3 = x2 * x;
207 double x4 = x3 * x;
208 double xm2 = (x - 1.) * (x - 1.);
209 double xm3 = xm2 * (x - 1.);
210 return ((-64. + 68. * x + 17. * x2 - 11. * x3) / 4. / xm2
211 + pow(M_PI, 2.) * 8. / 3. / x + 2. * gslpp_special_functions::dilog(1. - x) * (8. - 24. * x
212 + 20. * x2 - x3 + 7. * x4 - x4 * x) / (x * xm3)
213 + log(x)*(-32. + 68. * x - 32. * x2 + 28. * x3 - 3. * x4)
214 / (2. * xm3) + x2 * pow(log(x), 2.) * (4. - 7. * x + 7. * x2
215 - 2. * x3) / (2. * xm2 * xm2));
216}
217
218const double StandardModelMatching::S1(double x) const {
219 return (CF * S11(x) + (Nc - 1.) / 2. / Nc * S18(x));
220}
221
222/*******************************************************************************
223 * loop functions misiak base for b -> s gamma *
224 * operator basis: - current current *
225 * - qcd penguins *
226 * - magnetic and chromomagnetic penguins *
227 * - semileptonic *
228 * ****************************************************************************/
229
230const double StandardModelMatching::A0t(double x) const {
231 double x2 = x * x;
232 double x3 = x2 * x;
233
234 return ((-3. * x3 + 2. * x2) / (2. * pow(1. - x, 4.)) * log(x) +
235 (22. * x3 - 153. * x2 + 159. * x - 46.) / (36. * pow(1. - x, 3.)));
236}
237
238const double StandardModelMatching::B0t(double x) const {
239 return ( x / (4. * (1. - x) * (1. - x)) * log(x) + 1. / (4. * (1. - x)));
240}
241
242const double StandardModelMatching::C0t(double x) const {
243 return ( (3. * x * x + 2. * x) / (8. * (1. - x) * (1. - x)) * log(x) + (-x * x + 6. * x) / (8. * (1. - x)));
244}
245
246const double StandardModelMatching::D0t(double x) const {
247 double x2 = x * x;
248 double x3 = x2 * x;
249 double x4 = x3 * x;
250
251 return ( (-3. * x4 + 30. * x3 - 54. * x2 + 32. * x - 8.) / (18. * pow(1. - x, 4)) * log(x)
252 + (-47. * x3 + 237. * x2 - 312. * x + 104.) / (108. * pow(1. - x, 3.)));
253}
254
255const double StandardModelMatching::E0t(double x) const {
256 //***// CHECK THIS FUNCTION double x2 = x * x;
257
258 return (-9. * x * x + 16. * x - 4.) / (6. * pow((1. - x), 4)) * log(x) + (-7. * x * x * x - 21. * x * x + 42. * x + 4.) / (36 * pow((1. - x), 3));
259 //return (x * (18. - 11. * x - x * x) / (12. * pow(1. - x, 3.) + x * x * (15. - 16. * x + 4. * x * x) /(6. * pow(1. - x, 4.)) * log(x) - 2./3. * log(x)));
260}
261
262const double StandardModelMatching::F0t(double x) const {
263 double x2 = x * x;
264 double xm3 = (1. - x)*(1. - x)*(1. - x);
265
266 return ((3. * x2) / (2. * xm3 * (1. - x)) * log(x) + (5. * x2 * x - 9. * x2 + 30. * x - 8.) /
267 (12. * xm3));
268}
269
270const double StandardModelMatching::A1t(double x, double mu) const {
271 double x2 = x * x;
272 double x3 = x * x * x;
273 double x4 = x * x * x * x;
274 double xm2 = pow(1. - x, 2);
275 double xm3 = xm2 * (1. - x);
276 double xm4 = xm3 * (1. - x);
277 double xm5 = xm4 * (1. - x);
278 double mt = SM.Mrun(mu, SM.getQuarks(QCD::TOP).getMass_scale(),
279 SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, FULLNNLO);
280
281 return ((32. * x4 + 244. * x3 - 160. * x2 + 16. * x) / 9. / xm4 * gslpp_special_functions::dilog(1. - 1. / x) +
282 (-774. * x4 - 2826. * x3 + 1994. * x2 - 130. * x + 8.) / 81. / xm5 * log(x) +
283 (-94. * x4 - 18665. * x3 + 20682. * x2 - 9113. * x + 2006.) / 243. / xm4 +
284 ((-12. * x4 - 92. * x3 + 56. * x2) / 3. / (1. - x) / xm4 * log(x) +
285 (-68. * x4 - 202. * x3 - 804. * x2 + 794. * x - 152.) / 27. / xm4) * 2. * log(mu / mt));
286}
287
288const double StandardModelMatching::B1t(double x, double mu) const {
289 double x2 = x * x;
290 double xm2 = pow(1. - x, 2);
291 double xm3 = pow(1. - x, 3);
292 double mt = SM.Mrun(mu, SM.getQuarks(QCD::TOP).getMass_scale(),
293 SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, FULLNNLO);
294
295 return (-2. * x) / xm2 * gslpp_special_functions::dilog(1. - 1. / x) + (-x2 + 17. * x) / (3. * xm3) * log(x) + (13. * x + 3.) / (3. * xm2) +
296 ((2. * x2 + 2. * x) / xm3 * log(x) + (4. * x) / xm2) * 2. * log(mu / mt);
297 /*(2. * x)/xm2 * gsl_sf_dilog(1.-x) + (3. * x2 + x)/xm3 * log(x) * log(x) + (-11. * x2 - 5. * x)/(3. * xm3) * log(x) +
298 (-3. * x2 + 19. * x)/(3. * xm2) + 16. * x * (2. * (x - 1.) - (1. + x) * log(x))/(4. * xm3) * log(mu / Mw);*/
299}
300
301const double StandardModelMatching::C1t(double x, double mu) const {
302 double x2 = x * x;
303 double x3 = x * x2;
304 double xm2 = pow(1. - x, 2);
305 double xm3 = pow(1. - x, 3);
306 double mt = SM.Mrun(mu, SM.getQuarks(QCD::TOP).getMass_scale(),
307 SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, FULLNNLO);
308
309 return (-x3 - 4. * x) / xm2 * gslpp_special_functions::dilog(1. - 1. / x) + (3. * x3 + 14. * x2 + 23 * x) / (3. * xm3) * log(x) +
310 (4. * x3 + 7. * x2 + 29. * x) / (3. * xm2) + ((8. * x2 + 2. * x) / xm3 * log(x) + (x3 + x2 + 8. * x) / xm2) * 2. * log(mu / mt);
311 /*(x3 + 4. * x)/xm2 * gsl_sf_dilog(1.-x) + (x4 - x3 + 20. * x2)/(2. * xm3) * log(x) * log(x) +
312 (-3. * x4 - 3. * x3 - 35. * x2 + x)/(3. * xm3) * log(x) + (4. * x3 + 7. * x2 + 29. * x)/(3. * xm2) +
313 16. * x * (-8. + 7. * x + x3 - 2. * (1. + 4. * x) * log(x))/(8. * xm3) * log(mu / Mw);*/
314}
315
316const double StandardModelMatching::D1t(double x, double mu) const {
317 double x2 = x * x;
318 double x3 = x * x2;
319 double x4 = x * x3;
320 double xm4 = pow(1. - x, 4);
321 double xm5 = pow(1. - x, 5);
322 double mt = SM.Mrun(mu, SM.getQuarks(QCD::TOP).getMass_scale(),
323 SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, FULLNNLO);
324
325 return (380. * x4 - 1352. * x3 + 1656. * x2 - 784. * x + 256.) / (81. * xm4) * gslpp_special_functions::dilog(1. - 1. / x) +
326 (304. * x4 + 1716. * x3 - 4644. * x2 + 2768. * x - 720.) / (81. * xm5) * log(x) +
327 (-6175. * x4 + 41608. * x3 - 66723. * x2 + 33106. * x - 7000.) / (729. * xm4) +
328 ((648. * x4 - 720. * x3 - 232. * x2 - 160. * x + 32.) / (81. * xm5) * log(x) +
329 (-352. * x4 + 4912. * x3 - 8280. * x2 + 3304. * x - 880.) / (243. * xm4)) * 2. * log(mu / mt);
330}
331
332const double StandardModelMatching::F1t(double x, double mu) const {
333 double x2 = x * x;
334 double x3 = x * x2;
335 double x4 = x * x3;
336 double xm4 = pow(1. - x, 4);
337 double xm5 = pow(1. - x, 5);
338 double mt = SM.Mrun(mu, SM.getQuarks(QCD::TOP).getMass_scale(),
339 SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, FULLNNLO);
340
341 return ((4. * x4 - 40. * x3 - 41. * x2 - x) / 3. / xm4 * gslpp_special_functions::dilog(1. - 1. / x) +
342 (-144. * x4 + 3177. * x3 + 3661. * x2 + 250. * x - 32.) / 108. / xm5 * log(x)
343 + (-247. * x4 + 11890. * x3 + 31779. * x2 - 2966. * x + 1016.) / 648. / xm4
344 + ((17. * x3 + 31. * x2) / xm5 * log(x) + (-35. * x4 + 170. * x3 + 447. * x2
345 + 338. * x - 56.) / 18. / xm4)* 2. * log(mu / mt));
346}
347
348const double StandardModelMatching::E1t(double x, double mu) const
349 {
350 double x2 = x * x;
351 double x3 = x * x2;
352 double x4 = x * x3;
353 double xm4 = pow(1. - x, 4);
354 double xm5 = pow(1. - x, 5);
355 double mt = SM.Mrun(mu, SM.getQuarks(QCD::TOP).getMass_scale(),
356 SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, FULLNNLO);
357
358
359 return (515. * x4 - 614. * x3 - 81. * x2 - 190. * x + 40.) / (54. * xm4) * gslpp_special_functions::dilog(1. - 1. / x) +
360 (-1030. * x4 + 435. * x3 + 1373. * x2 + 1950. * x - 424.) / (108. * xm5) * log(x) +
361 (-29467. * x4 + 45604. * x3 - 30237. * x2 + 66532. * x - 10960.) / (1944. * xm4) +
362 ((-1125. * x3 + 1685. * x2 + 380. * x - 76.) / (54. * xm5) * log(x) +
363 (133. * x4 - 2758. * x3 - 2061. * x2 + 11522. * x - 1652.) / (324. * xm4)) * 2. * log(mu / mt);
364}
365
366const double StandardModelMatching::G1t(double x, double mu) const
367 {
368 double x2 = x * x;
369 double x3 = x * x2;
370 double x4 = x * x3;
371 double xm3 = pow(1. - x, 3);
372 double xm4 = pow(1. - x, 4);
373
374 double mt = SM.Mrun(mu, SM.getQuarks(QCD::TOP).getMass_scale(),
375 SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, FULLNNLO);
376
377 return (10. * x4 - 100. * x3 + 30. * x2 + 160. * x - 40.) / (27. * xm4) * gslpp_special_functions::dilog(1. - 1. / x) +
378 (30. * x3 - 42. * x2 - 332. * x + 68.) / (81. * xm4) * log(x) +
379 (-6. * x3 - 293. * x2 + 161. * x + 42.) / (81. * xm3) +
380 ((90. * x2 - 160. * x + 40.) / (27. * xm4) * log(x) +
381 (35. * x3 + 105. * x2 - 210. * x - 20.) / (81. * xm3)) * 2. * log(mu / mt);
382}
383
384const double StandardModelMatching::C7c_3L_at_mW(double x) const
385 {
386 double z = 1. / x;
387 return (1.525 - 0.1165 * z + 0.01975 * z * log(z) + 0.06283 * z * z + 0.005349 * z * z * log(z) + 0.01005 * z * z * log(z) * log(z)
388 - 0.04202 * z * z * z + 0.01535 * z * z * z * log(z) - 0.00329 * z * z * z * log(z) * log(z) + 0.002372 * z * z * z * z - 0.0007910 * z * z * z * z * log(z));
389}
390
391const double StandardModelMatching::C7t_3L_at_mt(double x) const
392 {
393 double z = 1. / x;
394 return (12.06 + 12.93 * z + 3.013 * z * log(z) + 96.71 * z * z + 52.73 * z * z * log(z)
395 + 147.9 * z * z * z + 187.7 * z * z * z * log(z) - 144.9 * z * z * z * z + 236.1 * z * z * z * z * log(z));
396}
397
398const double StandardModelMatching::C7t_3L_func(double x, double mu) const
399 {
400 double x2 = x * x;
401 double x3 = x * x2;
402 double x4 = x * x3;
403 double x5 = x * x4;
404 double xm1to5 = (x - 1.)*(x - 1.)*(x - 1.)*(x - 1.)*(x - 1.);
405 double xm1to6 = xm1to5 * (x - 1.);
406
407 double mt = SM.Mrun(mu, SM.getQuarks(QCD::TOP).getMass_scale(),
408 SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, FULLNNLO);
409
410 return ( 2. * log(mu / mt) * (gslpp_special_functions::dilog(1. - 1. / x) * (-592. * x5 - 22. * x4 + 12814. * x3 - 6376. * x2 + 512. * x) / 27. / xm1to5
411 + log(x) * (-26838. * x5 + 25938. * x4 + 627367. * x3 - 331956. * x2 + 16989. * x - 460.) / 729. / xm1to6
412 + (34400. * x5 + 276644. * x4 - 2668324. * x3 + 1694437. * x2 - 323354. * x + 53077.) / 2187. / xm1to5)
413 + 4. * log(mu / mt) * log(mu / mt) * (log(x)*(-63. * x5 + 532. * x4 + 2089. * x3 - 1118. * x2) / 9. / xm1to6
414 + (1186. * x5 - 2705. * x4 - 24791. * x3 - 16099. * x2 + 19229. * x - 2740.) / 162. / xm1to5));
415
416}
417
418const double StandardModelMatching::C8c_3L_at_mW(double x) const
419 {
420 double z = 1. / x;
421 return (-1.870 + 0.1010 * z - 0.1218 * z * log(z) + 0.1045 * z * z - 0.03748 * z * z * log(z)
422 + 0.01151 * z * z * log(z) * log(z) - 0.01023 * z * z * z + 0.004342 * z * z * z * log(z)
423 + 0.0003031 * z * z * z * log(z) * log(z) - 0.001537 * z * z * z * z + 0.0007532 * z * z * z * z * log(z));
424}
425
426const double StandardModelMatching::C8t_3L_at_mt(double x) const
427 {
428 double z = 1. / x;
429 return (-0.8954 - 7.043 * z - 98.34 * z * z - 46.21 * z * z * log(z) - 127.1 * z * z * z
430 - 181.6 * z * z * z * log(z) + 535.8 * z * z * z * z - 76.76 * z * z * z * z * log(z));
431}
432
433const double StandardModelMatching::C8t_3L_func(double x, double mu) const
434 {
435 double x2 = x * x;
436 double x3 = x * x2;
437 double x4 = x * x3;
438 double x5 = x * x4;
439 double xm1to5 = (x - 1.)*(x - 1.)*(x - 1.)*(x - 1.)*(x - 1.);
440 double xm1to6 = xm1to5 * (x - 1.);
441
442 double mt = SM.Mrun(mu, SM.getQuarks(QCD::TOP).getMass_scale(),
443 SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, FULLNNLO);
444
445 return ( 2. * log(mu / mt) * (gslpp_special_functions::dilog(1. - 1. / x) * (-148. * x5 + 1052. * x4 - 4811. * x3 - 3520. * x2 - 61. * x) / 18. / xm1to5
446 + log(x) * (-15984. * x5 + 152379. * x4 - 1358060. * x3 - 1201653. * x2 - 74190. * x + 9188.) / 1944. / xm1to6
447 + (109669. * x5 - 1112675. * x4 + 6239377. * x3 + 8967623. * x2 + 768722. * x - 42796.) / 11664. / xm1to5)
448 + 4. * log(mu / mt) * log(mu / mt) * (log(x) * (-139. * x4 - 2938. * x3 - 2683. * x2) / 12. / xm1to6
449 + (1295. * x5 - 7009. * x4 + 29495. * x3 + 64513. * x2 + 17458. * x - 2072.) / 216. / xm1to5));
450
451}
452
453const double StandardModelMatching::Tt(double x) const {
454 return ((-(16. * x + 8.) * sqrt(4. * x - 1.) * gslpp_special_functions::clausen(2. * asin(1. / 2. / sqrt(x)))) +((16. * x + 20. / 3.) * log(x)) + (32. * x) + (112. / 9.));
455}
456
457const double StandardModelMatching::Wt(double x) const {
458 double x2 = x * x;
459 double x3 = x * x * x;
460 double x4 = x * x * x * x;
461 double xm2 = pow(1. - x, 2);
462 double xm3 = xm2 * (1. - x);
463 double xm4 = xm3 * (1. - x);
464
465 return ((-32. * x4 + 38. * x3 + 15. * x2 - 18. * x) / 18. / xm4 * log(x) -
466 (-18. * x4 + 163. * x3 - 259. * x2 + 108. * x) / 36. / xm3);
467}
468
469const double StandardModelMatching::Eet(double x) const {
470 double x2 = x * x;
471 double xm2 = pow(1. - x, 2);
472 double xm3 = xm2 * (1. - x);
473 double xm4 = xm3 * (1. - x);
474
475 return ((x * (18. - 11. * x - x2)) / (12. * xm3) +
476 (log(x) * (x2 * (15. - 16. * x + 4. * x2)) / (6. * xm4)) - 2. * log(x) / 3.);
477}
478
479const double StandardModelMatching::Rest(double x, double mu) const
480 {
481 double mt = SM.Mrun(mu, SM.getQuarks(QCD::TOP).getMass_scale(), SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, FULLNNLO);
482
483 return (-37.01364013973161 + 7.870950908767437 * mt -
484 0.0015295355176062242 * mt * mt + 2.41071411865951 * Mw -
485 0.20348320015672194 * mt * Mw - 0.02858827491583899 * Mw * Mw +
486 0.001422822903303167 * mt * Mw * Mw + (4.257050684362808 + 0.17719711396626878 * mt -
487 0.8190947921716011 * Mw - 0.002315407459561656 * mt * Mw + 0.008797408866807221 * Mw * Mw) * log(
488 mu) + (0.49627858125619595 - pow(5.784745743815408, -8) * mt +
489 0.031869225004473686 * Mw - 0.00041193393986696286 * Mw * Mw) * log(
490 mu) * log(mu));
491}
492
493const double StandardModelMatching::Y0(double x) const {
494 return ( x / 8. * ((4 - 5 * x + x * x + 3 * x * log(x)) / pow(x - 1., 2.)));
495}
496
497const double StandardModelMatching::Y1(double x, double mu) const {
498 double x2 = x * x;
499 double x3 = x2 * x;
500 double x4 = x3 * x;
501 double logx = log(x);
502 double xm = 1. - x;
503 double xm2 = xm * xm;
504 double xm3 = xm2 * xm;
505
506 return ((10. * x + 10. * x2 + 4. * x3) / (3 * xm2) - (2. * x - 8. * x2 - x3 - x4) / xm3 * logx
507 + (2. * x - 14. * x2 + x3 - x4) / (2. * xm3) * pow(logx, 2.)
508 + (2. * x + x3) / xm2 * gslpp_special_functions::dilog(1. - x)
509 + 16. * x * (-4. + 3. * x + x3 - 6. * x * logx) / (8. * -xm3) * log(mu / Mw));
510}
511
512const double StandardModelMatching::C7LOeff(double x) const {
513 double x2 = x * x;
514 double x3 = x2 * x;
515
516 return ( (3. * x3 - 2. * x2) / (4. * pow(x - 1., 4.)) * log(x) + (-8. * x3 - 5. * x2 +
517 7. * x) / (24. * pow(x - 1., 3.)));
518}
519
520const double StandardModelMatching::C8LOeff(double x) const {
521 return ( -3. * x * x / (4. * pow(x - 1., 4.)) * log(x) + (-x * x * x + 5. * x * x + 2. * x) / (8. * pow(x - 1., 3)));
522}
523
524const double StandardModelMatching::C7NLOeff(double x) const {
525 double x2 = x * x;
526 double x3 = x2 * x;
527 double x4 = x3 * x;
528 double xm4 = pow(x - 1., 4.);
529 double xm5 = xm4 * (x - 1);
530 double logx = log(x);
531
532 double Li2 = gslpp_special_functions::dilog(1. - 1. / x);
533 return ( Li2 * (-16. * x4 - 122. * x3 + 80. * x2 - 8. * x) / (9. * xm4) +
534 (6. * x4 + 46. * x3 - 28. * x2) / (3. * xm5) * logx * logx +
535 (-102. * x4 * x - 588. * x4 - 2262. * x3 + 3244. * x2 - 1364. * x + 208.) / (81. * xm5) * logx +
536 (1646. * x4 + 12205. * x3 - 10740. * x2 + 2509. * x - 436.) / (486. * xm4));
537}
538
539const double StandardModelMatching::C8NLOeff(double x) const {
540 double x2 = x * x;
541 double x3 = x2 * x;
542 double x4 = x3 * x;
543 double xm4 = pow(x - 1., 4.);
544 double xm5 = xm4 * (x - 1);
545 double logx = log(x);
546
547 double Li2 = gslpp_special_functions::dilog(1. - 1. / x);
548 return (Li2 * (-4. * x4 + 40. * x3 + 41. * x2 + x) / (6. * xm4) +
549 (-17. * x3 - 31. * x2) / (2. * xm5) * logx * logx +
550 (-210. * x * x4 + 1086. * x4 + 4893. * x3 + 2857. * x2 - 1994. * x + 280.) / (216. * xm5) * logx +
551 (737. * x4 - 14102. * x3 - 28209. * x2 + 610. * x - 508.) / (1296. * xm4));
552}
553
554
555/******************************************************************************/
556
557/*******************************************************************************
558 * loop functions Buras base for nonlep. b decays *
559 * operator basis: - current current *
560 * - qcd penguins *
561 * - ew penguins *
562 * ****************************************************************************/
563
564const double StandardModelMatching::B0b(double x) const {
565 return ( 0.25 * (x / (1. - x) + x / (x * x - 2. * x + 1.) * log(x)));
566}
567
568const double StandardModelMatching::C0b(double x) const {
569 return ( x / 8. * ((x - 6.) / (x - 1.) + (3. * x + 2.) / (x * x - 2. * x + 1.) * log(x)));
570}
571
572const double StandardModelMatching::D0b(double x) const {
573 double x2 = x * x;
574 return ( -4. / 9. * log(x) + (-19. * x2 * x + 25. * x2) / (36. * (x2 * x - 3. * x2 + 3. * x - 1.))
575 + (x2 * (5. * x2 - 2. * x - 6.)) / (18. * pow(x - 1., 4.)) * log(x));
576}
577
578const double StandardModelMatching::D0b_tilde(double x) const {
579 return (D0b(x) - 4. / 9.);
580}
581
582const double StandardModelMatching::E0b(double x) const {
583 double x2 = x * x;
584
585 return ( -2. / 3. * log(x) + (18. * x - 11. * x2 - x2 * x) / (12. * (-x2 * x + 3. * x2 - 3. * x + 1.)) +
586 (x2 * (15. - 16. * x + 4. * x2)) / (6. * pow(1. - x, 4.)) * log(x));
587}
588
589/******************************************************************************
590 * Loop functions for QED corrections *
591 ******************************************************************************/
592
593const double StandardModelMatching::B1d(double x, double mu) const {
594 double xmo = x - 1.;
595 double dilog1mx = gslpp_special_functions::dilog(1. - x);
596 double mut = SM.getQuarks(QCD::TOP).getMass_scale();
597 double xmut = mut * mut / Mw / Mw;
598 double logxw;
599
600 if (mu == Muw)
601 logxw = L;
602 else
603 logxw = 2. * log(mu / Mw);
604
605 return (-(8. - 183. * x + 47. * x * x) / 24. / xmo / xmo - (8. + 27. * x + 93. * x * x) / 24. / xmo / xmo / xmo * log(x) +
606 (27. * x + 71. * x * x - 2. * x * x * x) / 24. / xmo / xmo / xmo * log(x) * log(x) - (2. - 3. * x - 9. * x * x + x * x * x) / 6. / x / xmo / xmo * dilog1mx +
607 (2. + x) / 36. / x * M_PI * M_PI + B0b(x)*(19. / 6. - logxw) + 4. * x / xmo / xmo * log(xmut) -
608 (2. * x + 2. * x * x) / xmo / xmo / xmo * log(x) * log(xmut));
609}
610
611const double StandardModelMatching::B1d_tilde(double x, double mu) const {
612 double xmo = x - 1.;
613 double dilog1mx = gslpp_special_functions::dilog(1. - x);
614 double logxw;
615
616 if (mu == Muw)
617 logxw = L;
618 else
619 logxw = 2. * log(mu / Mw);
620
621 return (-(8. - 23. * x) / 8. / xmo - (8. - 5. * x) / 8. / xmo / xmo * log(x) + (3. * x + 2. * x * x) / 8. / xmo / xmo * log(x) * log(x) +
622 (2. - 3. * x + 3. * x * x + x * x * x) / 2. / x / xmo / xmo * dilog1mx - (2. + x) / 12. / x * M_PI * M_PI + B0b(x)*(5. / 2. + 3. * logxw));
623}
624
625const double StandardModelMatching::B1u(double x, double mu) const {
626 double xmo = x - 1.;
627 double dilog1mx = gslpp_special_functions::dilog(1. - x);
628 double mut = SM.getQuarks(QCD::TOP).getMass_scale();
629 double xmut = mut * mut / Mw / Mw;
630 double logxw;
631
632 if (mu == Muw)
633 logxw = L;
634 else
635 logxw = 2. * log(mu / Mw);
636
637 return (-(46. * x + 18. * x * x) / 3. / xmo / xmo - (16. * x - 80. * x * x) / 3. / xmo / xmo / xmo * log(x) -
638 (9. * x + 23. * x * x) / 2. / xmo / xmo / xmo * log(x) * log(x) - 6. * x / xmo / xmo * dilog1mx +
639 B0b(x)*(-38. / 3. + 4. * logxw) - 16. * x / xmo / xmo * log(xmut) + (8. * x + 8. * x * x) / xmo / xmo / xmo * log(x) * log(xmut));
640}
641
642const double StandardModelMatching::B1u_tilde(double x, double mu) const {
643 double xmo = x - 1.;
644 double logx = log(x);
645 double dilog1mx = gslpp_special_functions::dilog(1. - x);
646 double logxw;
647
648 if (mu == Muw)
649 logxw = L;
650 else
651 logxw = 2. * log(mu / Mw);
652
653 return (-6. * x / xmo - 3. * x * logx * logx / 2. / xmo / xmo - 6. * x * dilog1mx / xmo / xmo - B0b(x)*(10. + 12. * logxw));
654}
655
656const double StandardModelMatching::C1ew(double x) const {
657 double xmo = x - 1.;
658 double dilog1mx = gslpp_special_functions::dilog(1. - x);
659 double mut = SM.getQuarks(QCD::TOP).getMass_scale();
660 double xmut = mut * mut / Mw / Mw;
661
662 return ((29. * x + 7. * x * x + 4. * x * x * x) / 3. / xmo / xmo + (x - 35. * x * x - 3. * x * x * x - 3. * x * x * x * x) / 3. / xmo / xmo / xmo * log(x) +
663 (20. * x * x - x * x * x + x * x * x * x) / 2. / xmo / xmo / xmo * log(x) * log(x) + (4. * x + x * x * x) / xmo / xmo * dilog1mx +
664 (8. * x + x * x + x * x * x) / xmo / xmo * log(xmut) - (2. * x + 8. * x * x) / xmo / xmo / xmo * log(x) * log(xmut));
665}
666
667const double StandardModelMatching::Zew(double xt, double xz) const {
668 double z0ew, z1ew;
669#ifdef ZEW_NUMERIC
670 z0ew = 5.1795 + 0.038 * (Mt_muw - 166.) + 0.015 * (Mw - 80.394);
671 z1ew = -2.1095 + 0.0067 * (Mt_muw - 166.) + 0.026 * (Mw - 80.394);
672#else
673 double xt2 = xt*xt;
674 double xt3 = xt2*xt;
675 double xz2 = xz*xz;
676 double xz3 = xz2*xz;
677 double xtmo = xt - 1.;
678 double xzmo = xz - 1.;
679 double dilog1mxt = gslpp_special_functions::dilog(1. - xt);
680 double dilog1mxz = gslpp_special_functions::dilog(1. - xz);
681 z0ew = -xt * (20. - 20. * xt2 - 457. * xz + 19. * xt * xz + 8. * xz2) / 32. / xtmo / xz +
682 xt * (10. * xt3 - 11. * xt2 * xz - xt * (30. - 16. * xz) + 4. * (5. - 17. * xz + xz2)) / 16. / xtmo / xtmo / xz * log(xt) +
683 xt * (10. - 10. * xt2 - 17. * xz - xt * xz - 4. * xz2) / 16. / xtmo / xz * log(xz) -
684 xz * (10. * xt2 - xt * (4. - xz) + 8. * xz) / 32. / xtmo / xtmo * log(xt) * log(xt) - xz2 / 4. * log(xz) * log(xz) -
685 ((8. + 12. * xt + xt2) / 4. / xz - 5. * xtmo * xtmo * (2. + xt) / 16. / xz2 -
686 (12. - 3. * xt3 - 3. * xt2 * (4. - xz) + 4. * xt * (3. - xz) + 4. * xz - xz2) / 8. / xtmo / xtmo) * log(xt) * log(xz) -
687 ((8. + 12. * xt + xt2) / 2. / xz - 5. * xtmo * xtmo * (2. + xt) / 8. / xz2 - 3. * (4. + 8. * xt + 2. * xt2 - xt3) / 4. / xtmo / xtmo) * dilog1mxt +
688 xzmo * xzmo * (5. - 6. * xz - 5. * xz2) / 4. / xz2 * dilog1mxz - (5. - 16. * xz + 12. * xz2 + 2 * xz3 * xz) / 24. / xz2 * M_PI * M_PI +
689 xt * (4. - xz)*(88. - 30. * xz - 25. * xz2 - 2. * xt * (44. - 5. * xz - 6. * xz2)) / 32. / xtmo / xtmo / xz * phi_z(xz / 4.) +
690 (16. * xt3 * xt - xt * (20. - xz) * xz2 + 8. * xz3 - 8. * xt3 * (14. + 5. * xz) + 8. * xt2 * (12. - 7. * xz + xz2)) / 32. / xtmo / xtmo / xz * phi_z(xz / 4. / xt) -
691 ((22. + 33. * xt - xt2) / 16. / xtmo / xz - 5. * xtmo * (2. + xt) / 16. / xz2 +
692 (2. + 5. * xt2 + 10. * xz + xt * (15. + xz)) / 16. / xtmo / xtmo) * phi_xy(xt, xz);
693
694 z1ew = xt * (20. - 20. * xt2 - 265. * xz + 67. * xt * xz + 8. * xz2) / 48. / xtmo / xz -
695 xt * (10. * xt3 - 15. * xt2 * xz + 4. * (5. - 7. * xz + 2. * xz2) - xt * (30. + 20. * xz + 4. * xz2)) / 24. / xtmo / xtmo / xz * log(xt) -
696 xt * (10. - 10. * xt2 - 33. * xz + 15. * xt * xz - 4. * xz2) / 24. / xtmo / xz * log(xz) +
697 xz * (8. - 16. * xt + 2. * xt2 + 10. * xz + 7. * xt * xz) / 48. / xtmo / xtmo * log(xt) * log(xt) + xz * (4. + 5. * xz) / 24. * log(xz) * log(xz) +
698 ((20. + 6. * xt + xt2) / 12. / xz - 5. * xtmo * xtmo * (2. + xt) / 24. / xz2 +
699 (3. * xt3 + 2. * xt2 * (12. - xz) - xt * (18. - 16. * xz + xz2) - 2. * (9. + 4. * xz - xz2)) / 12. / xtmo / xtmo) * log(xt) * log(xz) +
700 ((20. + 6. * xt + xt2) / 6. / xz - 5. * xtmo * xtmo * (2. + xt) / 12. / xz2 - (6. + 6. * xt - 8. * xt2 - xt3) / 2. / xtmo / xtmo) * dilog1mxt -
701 xzmo * xzmo * (5. - 10. * xz - 7. * xz2) / 6. / xz2 * dilog1mxz + (10. - 40. * xz + 36. * xz2 + 4. * xz3 + 5. * xz3 * xz) / 72. / xz2 * M_PI * M_PI +
702 xt * (xz - 4.)*(24. - 26. * xz - 13. * xz2 - 6. * xt * (4. - xz - xz2)) / 16. / xtmo / xtmo / xz * phi_z(xz / 4.) - (2. * xt2 * (2. + xt) / 3. / xtmo / xz -
703 (24. * xt3 + 12. * xt2 * (14. + xz) - 2. * xz * (4. + 5. * xz) - xt * (80. - 36. * xz + 7. * xz2)) / 48. / xtmo / xtmo) * phi_z(xz / 4. / xt) +
704 ((10. - xt - xt2) / 8. / xtmo / xz - 5 * xtmo * (2. + xt) / 24. / xz2 + (6. + 3. * xt2 + 14. * xz + 5. * xt * (7. - xz)) / 24. / xtmo / xtmo) * phi_xy(xt, xz);
705#endif
706 return (z0ew + sW2 * z1ew);
707}
708
709const double StandardModelMatching::Gew(double xt, double xz, double mu) const {
710 double xmuw = mu * mu / Mw / Mw;
711
712 return (Zew(xt, xz) + 5. * C0b(xt) + 6. * C0b(xt) * log(xmuw));
713}
714
715const double StandardModelMatching::Hew(double xt, double xz, double mu) const {
716 double xmuw = mu * mu / Mw / Mw;
717
718 return (Zew(xt, xz) - 7. * C0b(xt) + 6. * C0b(xt) * log(xmuw));
719}
720
721/******************************************************************************/
722/* loop functions for rare K and B decays, K-> pi nu nu & B-> Xs nu nu */
723
724/******************************************************************************/
725
726const double StandardModelMatching::X0t(double x) const {
727 return ((x / 8.) * ((x + 2.) / (x - 1.) + (3. * x - 6.) / ((x - 1.) * (x - 1.)) * log(x)));
728}
729
730const double StandardModelMatching::X1t(double x) const {
731
732 double x2 = x * x;
733 double x3 = x2 * x;
734 double x4 = x3 * x;
735 double xm3 = (1. - x) * (1. - x) * (1. - x);
736 double logx = log(x);
737
738 return (-(29. * x - x2 - 4. * x3) / (3. * (1. - x) * (1. - x))
739 - logx * (x + 9. * x2 - x3 - x4) / xm3
740 + logx * logx * (8. * x + 4. * x2 + x3 - x4) / (2. * xm3)
741 - gslpp_special_functions::dilog(1. - x) * (4. * x - x3) / ((1. - x) * (1. - x))
742 - 8. * x * log(Mut * Mut / Muw / Muw) * (8. - 9. * x + x3 + 6. * logx) / (8. * xm3));
743}
744
745const double StandardModelMatching::Xewt(double x, double a, double mu) const {
746 double b = 0.;
747 double swsq = SM.sW2_ND();
748
749 double A[17], C[17];
750
751 double x2 = x * x;
752 double x3 = x2 * x;
753 double x4 = x3 * x;
754 double x5 = x4 * x;
755 double x6 = x5 * x;
756 double x7 = x6 * x;
757 double x8 = x7 * x;
758 double M_PI2 = M_PI * M_PI;
759 double a2 = a * a;
760 double a3 = a2 * a;
761 double a4 = a3 * a;
762 double a5 = a4 * a;
763 double a6 = a5 * a;
764 double xm2 = (x - 1.) * (x - 1.);
765 double xm3 = xm2 * (x - 1.);
766 double axm = a * x - 1.;
767 double logx = log(x);
768 double loga = log(a);
769
770 A[0] = (16. - 48. * a) * M_PI2 + (288. * a - (32. - 88. * a) * M_PI2) * x
771 + (2003. * a + 4. * (4. - 6. * a - a2) * M_PI2) * x2
772 + (9. * a * (93. + 28. * a) - 4. * a * (3. - 2. * a + 8. * a2) * M_PI2) * x3
773 + (3. * a * (172. - 49. * a - 32. * a2) + 4. * a * (20. - a + 16. * a2) * M_PI2) * x4
774 - (3. * a * (168. + 11. * a - 24. * a2) + 4. * a * (45. + 8. * a2) * M_PI2) * x5
775 + 96. * a * M_PI2 * x6;
776
777 A[1] = -768. * x - (525. - 867. * a) * x2 + (303. + 318. * a) * x3 - 195. * a * x4;
778
779 A[2] = -8. * (95. - 67. * a + 11. * a2) * x2 + 2. * (662. - 78. * a - 177. * a2 + 40. * a3) * x3
780 - (608. + 476. * a - 595. * a2 + 114. * a3) * x4
781 + (44. + 188. * a - 321. * a2 + 103. * a3 - 8. * a4) * x5
782 - a * (28. - 72. * a + 33. * a2 - 4. * a3) * x6;
783
784 A[3] = +48. - 10. * (57. + 4. * a) * x + 51. * (29. + 10. * a) * x2 -
785 (841. + 1265. * a) * x3 + (308. + 347. * a) * x4
786 - (28. - 40. * a) * x5 + 12. * a * x6;
787
788 A[4] = +768. + (816. - 768. * a) * x + (1240. - 1232. * a) * x2
789 - 4. * (415. + 2. * a) * x3 + (311. + 722. * a) * x4
790 + (145. - 267. * a) * x5 - (36. + 51. * a) * x6 + 20. * a * x7;
791
792 A[5] = +328. * x - (536. + 900. * a) * x2 + (208. + 1584. * a + 670. * a2) * x3
793 - a * (668. + 1161. * a + 225. * a2) * x4
794 + a2 * (479. + 362. * a + 28. * a2) * x5
795 - a3 * (143. + 42. * a) * x6 + 16. * a4 * x7;
796
797 A[6] = +32. - 4. * (44. - 9. * a) * x + (384. - 322. * a - 400. * a2) * x2
798 - (400. - 869. * a - 1126. * a2 - 696. * a3) * x3
799 + 2. * (80. - 488. * a - 517. * a2 - 631. * a3 - 264. * a4) * x4
800 + (48. + 394. * a + 269. * a2 + 190. * a3 + 882. * a4 + 196. * a5) * x5
801 - (64. - 58. * a - 89. * a2 - 95. * a3 + 34. * a4 + 296. * a5 + 32. * a6) * x6
802 + (16. - 59. * a - 79. * a2 + 256. * a3 - 239. * a4
803 + 57. * a5 + 48. * a6) * x7
804 + (1. - a) * (1. - a) * (1. - a) * a2 * (29. + 16. * a) * x8;
805
806 A[7] = +28. * a2 * x2 - 32. * a3 * x3;
807
808 A[8] = -288. + 36. * (1. + 8. * a) * x + 6. * (647. + 87. * a) * x2 + 5. * (55. - 927. * a - 132. * a2) * x3
809 - (1233. + 98. * a - 879. * a2 - 192. * a3) * x4
810 + (360. + 1371. * a - 315. * a2 - 264. * a3) * x5
811 - 24. * a * (17. - 4. * a2) * x6;
812
813 A[9] = +32. + 4. * (-44. + 29. * a) * x - 12. * (-32. + 77. * a + 31. * a2) * x2
814 + 2. * (-200. + 837. * a + 767. * a2 + 182. * a3) * x3
815 - 2. * (-80. + 625. * a + 905. * a2 + 520. * a3 + 82. * a4) * x4
816 + (48. + 1079. * a + 590. * a2 + 1002. * a3 + 462. * a4 + 32. * a5) * x5
817 + (-64. - 1160. * a - 501. * a2 - 364. * a3 - 486. * a4 - 72. * a5) * x6
818 + (16. + 729. * a + 1038. * a2 + 38. * a3 + 238. * a4 + 52. * a5) * x7
819 - a * (192. + 743. * a + 50. * a3 + 12. * a4) * x8 + 192. * a2 * x8 * x;
820
821 A[10] = +16. * x + 324. * x2 - 36. * x4;
822
823 A[11] = +216. * x - 672. * x2 + 152. * x3;
824
825 A[12] = -16. * x + (16. - 42. * a) * x2 + (16. + 21. * a + 60. * a2) * x3
826 - (16 - 21. * a + 45. * a2 + 32. * a3) * x4 - a2 * (7. - 24. * a) * x5;
827
828 A[13] = -32. + (144. - 68. * a) * x + (-240. + 334. * a + 332. * a2) * x2
829 + (160. - 551. * a - 660. * a2 - 364. * a3) * x3
830 + a * (329. + 451. * a + 650. * a2 + 164. * a3) * x4
831 + (-48. - a - 59. * a2 - 523. * a3 - 316. * a4 - 32. * a5) * x5
832 + (16. - 43. * a - 93. * a2 + 255. * a3 + 287. * a4 + 32. * a5) * x6
833 - a2 * (-29. + 42. * a + 103. * a2 + 8. * a3) * x7;
834
835 A[14] = -144. * (1. - a)*(1. - a) * x2 + 144. * (1. - a) * (1. - a) * x3 - 36. * (1. - a) * (1. - a) * x4;
836
837 A[15] = -32. + 96. * a + (48. - 32. * a) * x - 176. * a * x2 - (16. - 74. * a) * x3 + 212. * a * x4;
838
839 A[16] = -32. + (64. - 100. * a) * x - 8. * (4. - 34. * a - 29. * a2) * x2
840 - 4. * a * (34. + 170. * a + 33. * a2) * x3
841 + 8. * a2 * (47. + 51. * a + 4. * a2) * x4 - 16. * a3 * (15. + 4. * a) * x5
842 + 32. * a4 * x6;
843
844 C[0] = 1. / (3. * a * xm2 * x);
845
846 C[1] = phi1(0.25) / (xm3 * axm);
847
848 C[2] = phi1(0.25 * a) / (2. * xm3 * axm);
849
850 C[3] = phi1(1. / 4. / x) / (2. * xm3 * axm);
851
852 C[4] = phi1(0.25 * x) / (2. * xm3 * axm);
853
854 C[5] = phi1(a * x * 0.25) / (xm3 * axm);
855
856 C[6] = phi2(1. / a / x, 1. / a) / (2. * a2 * x2 * xm3 * axm);
857
858 C[7] = loga * log(a) / axm;
859
860 C[8] = logx / (xm3 * axm * 3.);
861
862 C[9] = logx * logx / ((x - 1.) * xm3 * axm * 2. * a * x);
863
864 C[10] = 2. * log(mu / Mw) / xm2;
865
866 C[11] = logx * 2. * log(mu / Mw) / xm3;
867
868 C[12] = loga / (xm2 * axm);
869
870 C[13] = logx * loga / (2. * a * xm3 * x * axm);
871
872 C[14] = gslpp_special_functions::dilog(1. - a) / xm2;
873
874 C[15] = gslpp_special_functions::dilog(1. - x) / a / x;
875
876 C[16] = gslpp_special_functions::dilog(1. - a * x) / (a * x * xm2);
877
878 for (int i=0; i<17; i++){
879 b += C[i]*A[i];
880 }
881
882 return (b / 128. / swsq);
883}
884
885const double StandardModelMatching::phi1(double z) const {
886 if (z >= 0.) {
887 if (z < 1){
888 return(4. * sqrt(z / (1. - z)) * gslpp_special_functions::clausen(2. * asin(sqrt(z))));
889 }
890 else{
891 return((1. / sqrt(1. - 1. / z)) * (2. * log(0.5 * ( 1. - sqrt(1. - 1. / z) )) * log(0.5 * ( 1. - sqrt(1. - 1. / z) ))- 4. * gslpp_special_functions::dilog(0.5 * ( 1. - sqrt(1. - 1. / z) )) - log(4. * z) * log(4. * z)
892 + M_PI * M_PI / 3.));
893 }
894 } else {
895 std::stringstream out;
896 out << z;
897 throw std::runtime_error("StandardModelMatching::phi1(double z)" + out.str() + " <0");
898 }
899 return (0.);
900}
901
902const double StandardModelMatching::phi2(double x, double y) const{
903
904 double l2 = (1. - x - y) * (1. - x - y) - 4. * x * y;
905
906 if ( l2 >= 0. and (sqrt(x) + sqrt(y)) <= 1.){
907 double l = sqrt(l2);
908 return( 1. / l * (M_PI * M_PI/3. + 2. * log(0.5 * (1. + x - y - l)) * log(0.5 * (1. - x + y - l)) - log(x) * log(y) - 2. * gslpp_special_functions::dilog(0.5 * (1. + x - y - l)) - 2. * gslpp_special_functions::dilog(0.5 * (1. -x + y - l))));
909 }
910 else if(l2 < 0. and (sqrt(x) + sqrt(y)) > 1.){
911 return(2./sqrt( -l2) * (gslpp_special_functions::clausen(2. * acos((-1. + x + y)/(2. * sqrt(x * y))))
912 +gslpp_special_functions::clausen(2. * acos((1. + x - y) / (2. * sqrt(x))))
913 +gslpp_special_functions::clausen(2. * acos((1. - x + y) / (2. * sqrt(y))))));
914 }
915 else{
916 std::stringstream out;
917 out << x;
918 std::cout << "Error in StandardModelMatching::phi2(double x, double y): wrong input " << out.str() << std::endl;
919 SM.setSMSuccess(false);
920 }
921 return (0.);
922}
923
924const double StandardModelMatching::phi_z(double z) const {
925 double beta = sqrt(1. - 1. / z);
926 double clausen = gslpp_special_functions::clausen(2. * asin(sqrt(z)));
927 double dilog = gslpp_special_functions::dilog((1. - beta) / 2.);
928
929 if (z > 0.) {
930 if (z <= 1.) {
931 return (4. * sqrt(z / (1. - z)) * clausen);
932 } else {
933 return (1. / beta * (2. * log((1. - beta) / 2.) * log((1. - beta) / 2.) - 4. * dilog - log(4. * z) * log(4. * z) + M_PI * M_PI / 3.));
934 }
935 } else {
936 std::stringstream out;
937 out << z;
938 throw std::runtime_error("StandardModelMatching::phi_z(double z)" + out.str() + " <0");
939 }
940}
941
942const double StandardModelMatching::phi_xy(double x, double y) const {
943 double lambda = sqrt((1. - x - y)*(1. - x - y) - 4. * x * y);
944 double diloga = gslpp_special_functions::dilog((1. + x - y - lambda) / 2.);
945 double dilogb = gslpp_special_functions::dilog((1. - x + y - lambda) / 2.);
946 double clausenxy = gslpp_special_functions::clausen(2. * acos((-1. + x + y) / 2. / sqrt(x * y)));
947 double clausenx = gslpp_special_functions::clausen(2. * acos((1. + x - y) / 2. / sqrt(x)));
948 double clauseny = gslpp_special_functions::clausen(2. * acos((1. - x + y) / 2. / sqrt(y)));
949
950 if ((lambda * lambda) >= 0.) {
951 return (lambda * (2. * log((1. + x - y - lambda) / 2.) * log((1. - x + y - lambda) / 2.) - log(x) * log(y) -
952 2. * diloga - 2. * dilogb + M_PI * M_PI / 3.));
953 } else if ((lambda * lambda) < 0.) {
954 return (-2. * sqrt(-lambda * lambda)*(clausenxy + clausenx + clauseny));
955 } else {
956 std::stringstream out;
957 out << x;
958 throw std::runtime_error("StandardModelMatching::phi_xy(double x, double y) wrong" + out.str());
959 }
960}
961
962/*******************************************************************************
963 * Wilson coefficients Buras base for Delta B = 2 observables *
964 * ****************************************************************************/
965
966std::vector<WilsonCoefficient>& StandardModelMatching::CMdbd2() {
967 double gammam = 6. * CF;
968 double Bt;
969
970
971 double xt = x_t(Mut);
972 gslpp::complex co = GF / 4. / M_PI * Mw * SM.getCKM().computelamt_d();
973
974 vmcdb.clear();
975
976 switch (mcdbd2.getScheme()) {
977 case NDR:
978 Bt = BtNDR;
979 break;
980 case HV:
981 case LRI:
982 default:
983 std::stringstream out;
984 out << mcdbd2.getScheme();
985 throw std::runtime_error("StandardModelMatching::CMdb2(): scheme " + out.str() + "not implemented");
986 }
987
988 mcdbd2.setMu(Mut);
989
990 switch (mcdbd2.getOrder()) {
991 case NNLO:
992 case NLO:
993 mcdbd2.setCoeff(0, co * co * 4. * (SM.Als(Mut, FULLNLO) / 4. / M_PI * (S1(xt) + //* CHECK ORDER *//
994 (Bt + gamma0 * log(Mut / Mw)) * S0(xt, xt) + 2. * gammam * S0p(xt) * log(Mut / Mw))), NLO);
995#if SUSYFIT_DEBUG & 2
996 std::cout << "Mw = " << Mw << " mt(mut=" << Mut << ")= " << Mt_mut << " xt(mut=" << Mut << ")= " << xt << "matching of DB=2: S0(xt) = " << S0(xt) <<
997 ", S1(xt) = " << S1(xt) +
998 (Bt + gamma0 * log(Muw / Mw)) * S0(xt, xt) + 2. * gammam * S0p(xt) * log(Muw / Mw)
999 << ", lambdat_d^2 = " << SM.getCKM().computelamt_d() * SM.getCKM().computelamt_d()
1000 << std::endl;
1001#endif
1002 case LO:
1003 mcdbd2.setCoeff(0, co * co * 4. * (S0(xt, xt)), LO);
1004 break;
1005 default:
1006 std::stringstream out;
1007 out << mcdbd2.getOrder();
1008 throw std::runtime_error("StandardModelMatching::CMdbd2(): order " + out.str() + "not implemented");
1009 }
1010
1011
1012 vmcdb.push_back(mcdbd2);
1013 return (vmcdb);
1014}
1015
1016std::vector<WilsonCoefficient>& StandardModelMatching::CMdbs2() {
1017 double gammam = 6. * CF;
1018 double Bt;
1019 double xt = x_t(Mut);
1020 gslpp::complex co = GF / 4. / M_PI * Mw * SM.getCKM().computelamt_s();
1021
1022 vmcds.clear();
1023
1024 switch (mcdbs2.getScheme()) {
1025 case NDR:
1026 Bt = BtNDR;
1027 break;
1028 case HV:
1029 case LRI:
1030 default:
1031 std::stringstream out;
1032 out << mcdbs2.getScheme();
1033 throw std::runtime_error("StandardModelMatching::CMdbs2(): scheme " + out.str() + "not implemented");
1034 }
1035
1036 mcdbs2.setMu(Mut);
1037
1038 switch (mcdbs2.getOrder()) {
1039 case NNLO:
1040 case NLO:
1041 mcdbs2.setCoeff(0, co * co * 4. * (SM.Als(Mut, FULLNLO) / 4. / M_PI * (S1(xt) + //* CHECK ORDER *//
1042 (Bt + gamma0 * log(Mut / Mw)) * S0(xt, xt) + 2. * gammam * S0p(xt) * log(Mut / Mw))), NLO);
1043 case LO:
1044 mcdbs2.setCoeff(0, co * co * 4. * (S0(xt, xt)), LO);
1045 break;
1046 default:
1047 std::stringstream out;
1048 out << mcdbs2.getOrder();
1049 throw std::runtime_error("StandardModelMatching::CMdbs2(): order " + out.str() + "not implemented");
1050 }
1051
1052 vmcds.push_back(mcdbs2);
1053 return (vmcds);
1054}
1055
1056/*******************************************************************************
1057 * Wilson coefficients Buras base for Delta S = 2 observables *
1058 * Comment: they look empty because they are computed in Flavour/EvolDF2.cpp *
1059 * due to historical reasons. *
1060 * ****************************************************************************/
1061
1062std::vector<WilsonCoefficient>& StandardModelMatching::CMdk2() {
1063 vmck2.clear();
1064
1065 mcdk2.setMu(Mut);
1066
1067 switch (mcdk2.getOrder()) {
1068 case NNLO:
1069 case NLO:
1070 mcdk2.setCoeff(0, 0., NLO);
1071 case LO:
1072 mcdk2.setCoeff(0, 0., LO);
1073 break;
1074 default:
1075 std::stringstream out;
1076 out << mcdk2.getOrder();
1077 throw std::runtime_error("StandardModelMatching::CMdk2(): order " + out.str() + "not implemented");
1078 }
1079
1080 vmck2.push_back(mcdk2);
1081 return (vmck2);
1082}
1083
1084std::vector<WilsonCoefficient>& StandardModelMatching::CMd1Buras() {
1085 vmcd1Buras.clear();
1086
1087 switch (mcd1Buras.getScheme()) {
1088 case NDR:
1089 //case HV:
1090 //case LRI:
1091 break;
1092 default:
1093 std::stringstream out;
1094 out << mcd1Buras.getScheme();
1095 throw std::runtime_error("StandardModelMatching::CMd1Buras(): scheme " + out.str() + "not implemented");
1096 }
1097
1098 mcd1Buras.setMu(Muw);
1099
1100 switch (mcd1Buras.getOrder()) {
1101 case NNLO:
1102 case NLO:
1103 mcd1Buras.setCoeff(0, SM.Als(Muw, FULLNLO) / 4. / M_PI * 11. / 2., NLO); //* CHECK ORDER *//
1104 mcd1Buras.setCoeff(1, SM.Als(Muw, FULLNLO) / 4. / M_PI * (-11. / 6.), NLO);
1105 for (int j = 2; j < 10; j++) {
1106 mcd1Buras.setCoeff(j, 0., NLO);
1107 }
1108 case LO:
1109 mcd1Buras.setCoeff(0, 0., LO);
1110 mcd1Buras.setCoeff(1, 1., LO);
1111 for (int j = 2; j < 10; j++) {
1112 mcd1Buras.setCoeff(j, 0., LO);
1113 }
1114 break;
1115 default:
1116 std::stringstream out;
1117 out << mcd1Buras.getOrder();
1118 throw std::runtime_error("StandardModelMatching::CMd1Buras(): order " + out.str() + "not implemented");
1119 }
1120
1121 vmcd1Buras.push_back(mcd1Buras);
1122
1123 return (vmcd1Buras);
1124
1125}
1126
1127std::vector<WilsonCoefficient>& StandardModelMatching::CMd1() {
1128 vmcd1.clear();
1129
1130 switch (mcd1.getScheme()) {
1131 case NDR:
1132 //case HV:
1133 //case LRI:
1134 break;
1135 default:
1136 std::stringstream out;
1137 out << mcd1.getScheme();
1138 throw std::runtime_error("StandardModelMatching::CMd1(): scheme " + out.str() + "not implemented");
1139 }
1140
1141 mcd1.setMu(Muw);
1142
1143 switch (mcd1.getOrder()) {
1144 case NNLO:
1145 case NLO:
1146 mcd1.setCoeff(0, SM.Als(Muw, FULLNLO) / 4. / M_PI * 15., NLO); //* CHECK ORDER *//
1147 for (int j = 1; j < 10; j++) {
1148 mcd1.setCoeff(j, 0., NLO);
1149 }
1150 case LO:
1151 mcd1.setCoeff(0, 0., LO);
1152 mcd1.setCoeff(1, 1., LO);
1153 for (int j = 2; j < 10; j++) {
1154 mcd1.setCoeff(j, 0., LO);
1155 }
1156 break;
1157 default:
1158 std::stringstream out;
1159 out << mcd1.getOrder();
1160 throw std::runtime_error("StandardModelMatching::CMd1(): order " + out.str() + "not implemented");
1161 }
1162
1163 vmcd1.push_back(mcd1);
1164
1165 return (vmcd1);
1166
1167}
1168
1169std::vector<WilsonCoefficient>& StandardModelMatching::CMdd2() {
1170 vmcd2.clear();
1171
1172 switch (mcdd2.getScheme()) {
1173 case NDR:
1174 case HV:
1175 case LRI:
1176 break;
1177 default:
1178 std::stringstream out;
1179 out << mcdd2.getScheme();
1180 throw std::runtime_error("StandardModelMatching::CMdd2(): scheme " + out.str() + "not implemented");
1181 }
1182
1183 mcdd2.setMu(Muw);
1184
1185 switch (mcdd2.getOrder()) {
1186 case NNLO:
1187 case NLO:
1188 for (int i = 0; i < 5; i++)
1189 mcdd2.setCoeff(i, 0., NLO);
1190 case LO:
1191 for (int j = 0; j < 5; j++)
1192 mcdd2.setCoeff(j, 0., LO);
1193 break;
1194 default:
1195 std::stringstream out;
1196 out << mcdd2.getOrder();
1197 throw std::runtime_error("StandardModelMatching::CMdd2(): order " + out.str() + "not implemented");
1198 }
1199
1200 vmcd2.push_back(mcdd2);
1201 return (vmcd2);
1202}
1203
1204std::vector<WilsonCoefficient>& StandardModelMatching::CMK() {
1205
1206 double xt = x_t(Muw);
1207
1208 vmck.clear();
1209
1210 switch (mck.getScheme()) {
1211 case NDR:
1212 //case HV:
1213 //case LRI:
1214 break;
1215 default:
1216 std::stringstream out;
1217 out << mck.getScheme();
1218 throw "StandardModel::CMK(): scheme " + out.str() + "not implemented";
1219 }
1220
1221 mck.setMu(Muw);
1222
1223 switch (mck.getOrder()) {
1224 case NNLO:
1225 case NLO:
1226 for (int j = 0; j < 10; j++) {
1227 mck.setCoeff(j, SM.Als(Muw, FULLNLO) / 4. / M_PI * setWCbnlep(j, xt, NLO), NLO);
1228 mck.setCoeff(j, Ale / 4. / M_PI * setWCbnlepEW(j, xt), NLO_QED11);
1229 }
1230 case LO:
1231 for (int j = 0; j < 10; j++) {
1232 mck.setCoeff(j, setWCbnlep(j, xt, LO), LO);
1233 mck.setCoeff(j, 0., LO_QED);
1234 }
1235 break;
1236 default:
1237 std::stringstream out;
1238 out << mck.getOrder();
1239 throw "StandardModelMatching::CMK(): order " + out.str() + "not implemented";
1240 }
1241
1242 vmck.push_back(mck);
1243 return (vmck);
1244}
1245
1246/*******************************************************************************
1247 * Wilson coefficients Buras base for K -> pi pi decays *
1248 * operator basis: - current current *
1249 * ****************************************************************************/
1250std::vector<WilsonCoefficient>& StandardModelMatching::CMKCC() {
1251
1252 double xt = x_t(Muw);
1253
1254 vmckcc.clear();
1255
1256 switch (mckcc.getScheme()) {
1257 case NDR:
1258 //case HV:
1259 //case LRI:
1260 break;
1261 default:
1262 std::stringstream out;
1263 out << mckcc.getScheme();
1264 throw "StandardModelMatching::CMKCC(): scheme " + out.str() + "not implemented";
1265 }
1266
1267 mckcc.setMu(Muw);
1268
1269 switch (mckcc.getOrder()) {
1270 case NNLO:
1271 case NLO:
1272 for (int j = 0; j < 2; j++) {
1273 mckcc.setCoeff(j, SM.Als(Muw, FULLNLO) / 4. / M_PI * setWCbnlep(j, xt, NLO), NLO);
1274 mckcc.setCoeff(j, Ale / 4. / M_PI * setWCbnlepEW(j, xt), NLO_QED11);
1275 }
1276 for (int j = 2; j < 10; j++) {
1277 mckcc.setCoeff(j, 0., NLO);
1278 mckcc.setCoeff(j, 0., NLO_QED11);
1279 }
1280 case LO:
1281 for (int j = 0; j < 2; j++) {
1282 mckcc.setCoeff(j, setWCbnlep(j, xt, LO), LO);
1283 mckcc.setCoeff(j, 0., LO_QED);
1284 }
1285 for (int j = 2; j < 10; j++) {
1286 mckcc.setCoeff(j, 0., LO);
1287 mckcc.setCoeff(j, 0., LO_QED);
1288 }
1289 break;
1290 default:
1291 std::stringstream out;
1292 out << mckcc.getOrder();
1293 throw "StandardModelMatching::CMKCC(): order " + out.str() + "not implemented";
1294 }
1295
1296 vmckcc.push_back(mckcc);
1297 return (vmckcc);
1298}
1299
1300/*******************************************************************************
1301 * Wilson coefficients misiak base for b -> s g *
1302 * operator basis: - current current *
1303 * - qcd penguins *
1304 * - magnetic and chromomagnetic penguins *
1305 * ****************************************************************************/
1306std::vector<WilsonCoefficient>& StandardModelMatching::CMbsg() {
1307 double xt = x_t(Muw);
1308
1309 vmcbsg.clear();
1310
1311 switch (mcbsg.getScheme()) {
1312 case NDR:
1313 //case HV:
1314 //case LRI:
1315 break;
1316 default:
1317 std::stringstream out;
1318 out << mcbsg.getScheme();
1319 throw std::runtime_error("StandardModelMatching::CMbsg(): scheme " + out.str() + "not implemented");
1320 }
1321
1322 mcbsg.setMu(Muw);
1323 double alSo4pi = SM.Alstilde5(Muw);
1324
1325 switch (mcbsg.getOrder()) {
1326 case NNLO:
1327 for (int j = 0; j < 8; j++) {
1328 mcbsg.setCoeff(j, alSo4pi * alSo4pi * setWCbsg(j, xt, NNLO), NNLO);
1329 }
1330 case NLO:
1331 for (int j = 0; j < 8; j++) {
1332 mcbsg.setCoeff(j, alSo4pi * setWCbsg(j, xt, NLO), NLO);
1333 }
1334 case LO:
1335 for (int j = 0; j < 8; j++) {
1336 mcbsg.setCoeff(j, setWCbsg(j, xt, LO), LO);
1337 }
1338 break;
1339 default:
1340 std::stringstream out;
1341 out << mcbsg.getOrder();
1342 throw std::runtime_error("StandardModelMatching::CMbsg(): order " + out.str() + "not implemented");
1343 }
1344
1345 vmcbsg.push_back(mcbsg);
1346 return (vmcbsg);
1347}
1348
1349std::vector<WilsonCoefficient>& StandardModelMatching::CMprimebsg() {
1350 vmcprimebsg.clear();
1351
1352 switch (mcprimebsg.getScheme()) {
1353 case NDR:
1354 //case HV:
1355 //case LRI:
1356 break;
1357 default:
1358 std::stringstream out;
1359 out << mcprimebsg.getScheme();
1360 throw std::runtime_error("StandardModelMatching::CMprimebsg(): scheme " + out.str() + "not implemented");
1361 }
1362
1363 mcprimebsg.setMu(Muw);
1364
1365 switch (mcprimebsg.getOrder()) {
1366 case NNLO:
1367 for (int j = 0; j < 8; j++) {
1368 mcprimebsg.setCoeff(j, 0., NNLO); //* CHECK ORDER *//
1369 }
1370 case NLO:
1371 for (int j = 0; j < 8; j++) {
1372 mcprimebsg.setCoeff(j, 0., NLO); //* CHECK ORDER *//
1373 }
1374 case LO:
1375 for (int j = 0; j < 8; j++) {
1376 mcprimebsg.setCoeff(j, 0., LO);
1377 }
1378 break;
1379 default:
1380 std::stringstream out;
1381 out << mcprimebsg.getOrder();
1382 throw std::runtime_error("StandardModelMatching::CMprimebsg(): order " + out.str() + "not implemented");
1383 }
1384
1385 vmcprimebsg.push_back(mcprimebsg);
1386 return (vmcprimebsg);
1387}
1388
1389/*******************************************************************************
1390 * Wilson coefficients calcoulus, misiak base for b -> s gamma *
1391 * ****************************************************************************/
1392
1393double StandardModelMatching::setWCbsg(int i, double x, orders order) {
1394 sw = sqrt(sW2);
1395
1396 if (swf == sw && xcachef == x) {
1397 switch (order) {
1398 case NNLO:
1399 return ( CWbsgArrayNNLO[i]);
1400 break;
1401 case NLO:
1402 return ( CWbsgArrayNLO[i]);
1403 break;
1404 case LO:
1405 return ( CWbsgArrayLO[i]);
1406 break;
1407 default:
1408 std::stringstream out;
1409 out << order;
1410 throw std::runtime_error("order" + out.str() + "not implemeted");
1411 }
1412 }
1413
1414 swf = sw;
1415 xcachef = x;
1416
1417 switch (order) {
1418 case NNLO:
1419 // One and two loop matching from arXiv:9910220.
1420 CWbsgArrayNNLO[0] = -Tt(x) + 7987. / 72. + 17. / 3. * M_PI * M_PI + 475. / 6. * L + 17. * L*L;
1421 CWbsgArrayNNLO[1] = 127. / 18. + 4. / 3. * M_PI * M_PI + 46. / 3. * L + 4. * L*L;
1422 CWbsgArrayNNLO[2] = G1t(x, Muw) - 680. / 243. - 20. / 81. * M_PI * M_PI - 68. / 81. * L - 20. / 27 * L*L;
1423 CWbsgArrayNNLO[3] = E1t(x, Muw) + 950. / 243. + 10. / 81. * M_PI * M_PI + 124. / 27. * L + 10. / 27. * L*L;
1424 CWbsgArrayNNLO[4] = -0.1 * G1t(x, Muw) + 2. / 15. * E0t(x) + 68. / 243. + 2. / 81. * M_PI * M_PI + 14. / 81. * L + 2. / 27. * L*L;
1425 CWbsgArrayNNLO[5] = -3. / 16. * G1t(x, Muw) + 0.25 * E0t(x) + 85. / 162. + 5. / 108. * M_PI * M_PI + 35. / 108. * L + 5. / 36 * L*L;
1426 // 3-loop matching from arXiv:0401041. Expansion around x = 1, on the basis of Fig.3 of the paper.
1427 CWbsgArrayNNLO[6] = (C7t_3L_at_mt(x) + C7t_3L_func(x, Muw)-(C7c_3L_at_mW(x) + 13763. / 2187. * L + 814. / 729. * L * L)) - 1. / 3. * CWbsgArrayNNLO[2] - 4. / 9. * CWbsgArrayNNLO[3] - 20. / 3. * CWbsgArrayNNLO[4] - 80. / 9. * CWbsgArrayNNLO[5]; //effective C7
1428 CWbsgArrayNNLO[7] = (C8t_3L_at_mt(x) + C8t_3L_func(x, Muw)-(C8c_3L_at_mW(x) + 16607. / 5832. * L + 397. / 486. * L * L)) + CWbsgArrayNNLO[2] - 1. / 6. * CWbsgArrayNNLO[3] - 20. * CWbsgArrayNNLO[4] - 10. / 3. * CWbsgArrayNNLO[5]; //effective C8
1429 case NLO:
1430 CWbsgArrayNLO[0] = 15. + 6 * L;
1431 CWbsgArrayNLO[3] = E0t(x) - 7. / 9. + 2. / 3. * L;
1432 CWbsgArrayNLO[6] = -0.5 * A1t(x, Muw) + 713. / 243. + 4. / 81. * L - 4. / 9. * CWbsgArrayNLO[3]; //effective C7
1433 CWbsgArrayNLO[7] = -0.5 * F1t(x, Muw) + 91. / 324. - 4. / 27. * L - 1. / 6. * CWbsgArrayNLO[3]; //effective C8
1434 case LO:
1435 CWbsgArrayLO[1] = 1.;
1436 CWbsgArrayLO[6] = -0.5 * A0t(x) - 23. / 36.; //effective C7
1437 CWbsgArrayLO[7] = -0.5 * F0t(x) - 1. / 3.; //effective C8
1438 break;
1439 default:
1440 std::stringstream out;
1441 out << order;
1442 throw std::runtime_error("order" + out.str() + "not implemeted");
1443 }
1444
1445 switch (order) {
1446 case NNLO:
1447 return ( CWbsgArrayNNLO[i]);
1448 break;
1449 case NLO:
1450 return ( CWbsgArrayNLO[i]);
1451 break;
1452 case LO:
1453 return ( CWbsgArrayLO[i]);
1454 break;
1455 default:
1456 std::stringstream out;
1457 out << order;
1458 throw std::runtime_error("order" + out.str() + "not implemeted");
1459 }
1460}
1461
1462/*******************************************************************************
1463 * Wilson coefficients misiak base for B -> K^*ll *
1464 * operator basis: - current current *
1465 * - qcd penguins *
1466 * - magnetic and chromomagnetic penguins *
1467 * - semileptonic *
1468 * ****************************************************************************/
1469std::vector<WilsonCoefficient>& StandardModelMatching::CMBMll(QCD::lepton lepton) {
1470 double xt = x_t(Muw); //* ORDER FULLNNLO*//
1471
1472 vmcBMll.clear();
1473
1474 switch (mcBMll.getScheme()) {
1475 case NDR:
1476 //case HV:
1477 //case LRI:
1478 break;
1479 default:
1480 std::stringstream out;
1481 out << mcBMll.getScheme();
1482 throw std::runtime_error("StandardModelMatching::CMBMll(): scheme " + out.str() + "not implemented");
1483 }
1484
1485 mcBMll.setMu(Muw);
1486
1487 switch (mcBMll.getOrder()) {
1488 case NNLO:
1489 case NLO:
1490 for (int j = 0; j < 13; j++) {
1491 mcBMll.setCoeff(j, SM.Als(Muw, FULLNNLO) / 4. / M_PI * setWCBMll(j, xt, NLO), NLO);
1492 }
1493 case LO:
1494 for (int j = 0; j < 13; j++) {
1495 mcBMll.setCoeff(j, setWCBMll(j, xt, LO), LO);
1496 }
1497 break;
1498 default:
1499 std::stringstream out;
1500 out << mcBMll.getOrder();
1501 throw std::runtime_error("StandardModelMatching::CMBMll(): order " + out.str() + "not implemented");
1502 }
1503
1504 vmcBMll.push_back(mcBMll);
1505 return (vmcBMll);
1506}
1507
1508/*******************************************************************************
1509 * Wilson coefficients calcoulus, misiak base for B -> K^*ll *
1510 * *****************************************************************************/
1511
1512double StandardModelMatching::setWCBMll(int i, double x, orders order) {
1513 sw = sqrt(sW2);
1514
1515 if (swa == sw && xcachea == x) {
1516 switch (order) {
1517 case NNLO:
1518 case NLO:
1519 return ( CWBMllArrayNLO[i]);
1520 break;
1521 case LO:
1522 return ( CWBMllArrayLO[i]);
1523 break;
1524 default:
1525 std::stringstream out;
1526 out << order;
1527 throw std::runtime_error("order" + out.str() + "not implemeted");
1528 }
1529 }
1530
1531 swa = sw;
1532 xcachea = x;
1533
1534 switch (order) {
1535 case NNLO:
1536 case NLO:
1537 CWBMllArrayNLO[0] = 15. + 6 * L;
1538 CWBMllArrayNLO[3] = E0t(x) - (7. / 9.) + (2. / 3. * L);
1539 CWBMllArrayNLO[6] = -0.5 * A1t(x, Muw) + 713. / 243. + 4. / 81. * L - 4. / 9. * CWBMllArrayNLO[3];
1540 CWBMllArrayNLO[7] = -0.5 * F1t(x, Muw) + 91. / 324. - 4. / 27. * L - 1. / 6. * CWBMllArrayNLO[3];
1541 CWBMllArrayNLO[8] = (1 - 4. * sW2) / (sW2) * C1t(x, Muw) - 1. / (sW2) * B1t(x, Muw) - D1t(x, Muw) + 1. / sW2 + 524. / 729. -
1542 128. * M_PI * M_PI / 243. - 16. * L / 3. - 128. * L * L / 81.;
1543 CWBMllArrayNLO[9] = (B1t(x, Muw) - C1t(x, Muw)) / sW2 - 1. / sW2;
1544 case LO:
1545 CWBMllArrayLO[1] = 1.;
1546 CWBMllArrayLO[6] = -0.5 * A0t(x) - 23. / 36.;
1547 CWBMllArrayLO[7] = -0.5 * F0t(x) - 1. / 3.;
1548 CWBMllArrayLO[8] = (1 - 4. * sW2) / (sW2) * C0t(x) - 1. / (sW2) * B0t(x) - D0t(x) + 38. / 27. + 1. / (4. * sW2) - (4. / 9.) * L + 8. / 9. * log(SM.getMuw() / mu_b);
1549 CWBMllArrayLO[9] = 1. / (sW2) * (B0t(x) - C0t(x)) - 1. / (4. * sW2);
1550 break;
1551 default:
1552 std::stringstream out;
1553 out << order;
1554 throw std::runtime_error("order" + out.str() + "not implemeted");
1555 }
1556
1557 switch (order) {
1558 case NNLO:
1559 case NLO:
1560 return ( CWBMllArrayNLO[i]);
1561 break;
1562 case LO:
1563 return ( CWBMllArrayLO[i]);
1564 break;
1565 default:
1566 std::stringstream out;
1567 out << order;
1568 throw std::runtime_error("order" + out.str() + "not implemeted");
1569 }
1570}
1571
1572/*******************************************************************************
1573 * Wilson coefficients misiak primed base for B -> K^*ll *
1574 * operator basis: - current current *
1575 * - qcd penguins *
1576 * - magnetic and chromomagnetic penguins *
1577 * - semileptonic *
1578 * ****************************************************************************/
1579std::vector<WilsonCoefficient>& StandardModelMatching::CMprimeBMll(QCD::lepton lepton) {
1580 vmcprimeBMll.clear();
1581 mcprimeBMll.setMu(Muw);
1582 switch (mcprimeBMll.getOrder()) {
1583 case NNLO:
1584 case NLO:
1585 for (int j = 0; j < 13; j++) {
1586 mcprimeBMll.setCoeff(j, 0., NLO);
1587 }
1588 case LO:
1589 for (int j = 0; j < 13; j++) {
1590 mcprimeBMll.setCoeff(j, 0., LO);
1591 }
1592 break;
1593 default:
1594 std::stringstream out;
1595 out << mcprimeBMll.getOrder();
1596 throw std::runtime_error("StandardModelMatching::CMprimeBMll(): order " + out.str() + "not implemented");
1597 }
1598 vmcprimeBMll.push_back(mcprimeBMll);
1599 return (vmcprimeBMll);
1600}
1601
1602
1603/******************************************************************************/
1604
1605/*******************************************************************************
1606 * Wilson coefficients Misiak base for bs -> mu mu. decays *
1607 * operator basis: - current current *
1608 * - qcd penguins *
1609 * - semileptonic *
1610 * ****************************************************************************/
1611
1612std::vector<WilsonCoefficient>& StandardModelMatching::CMbsmm() {
1613
1614 // The couplings are not used here, but in the Bsmumu class.
1615
1616 double xt = x_t(Muw);
1617
1618 vmcbsmm.clear();
1619
1620
1621
1622 switch (mcbsmm.getScheme()) {
1623 case NDR:
1624 //case HV:
1625 //case LRI:
1626 break;
1627 default:
1628 std::stringstream out;
1629 out << mcbsmm.getScheme();
1630 throw std::runtime_error("StandardModelMatching::CMbsmm(): scheme " + out.str() + "not implemented");
1631 }
1632
1633 mcbsmm.setMu(Muw);
1634
1635 switch (mcbsmm.getOrder()) {
1636 case NNLO:
1637
1638 for (int j = 0; j < 8; j++) {
1639
1640 mcbsmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 1)) *
1641 setWCBsmm(j, xt, NNLO), NNLO);
1642
1643 }
1644 case NLO:
1645
1646 for (int j = 0; j < 8; j++) {
1647
1648 mcbsmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 1)) *
1649 setWCBsmm(j, xt, NLO), NLO);
1650 }
1651 case LO:
1652
1653 for (int j = 0; j < 8; j++) {
1654
1655 mcbsmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 1)) * setWCBsmm(j, xt, LO), LO);
1656 }
1657 break;
1658 default:
1659 std::stringstream out;
1660 out << mcbsmm.getOrder();
1661 throw std::runtime_error("StandardModelMatching::CMbsmm(): order " + out.str() + "not implemented");
1662 }
1663
1664 switch (mcbsmm.getOrder_qed()) {
1665 case NLO_QED22:
1666 ;
1667 for (int j = 0; j < 8; j++) {
1668
1669 mcbsmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 1)) * setWCBsmmEW(j, xt, NLO_QED22), NLO_QED22);
1670
1671 }
1672 case NLO_QED12: /*absent at high energy */
1673 for (int j = 0; j < 8; j++) {
1674
1675 mcbsmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 1)) * 0., NLO_QED12);
1676 }
1677 case NLO_QED21:
1678
1679 for (int j = 0; j < 8; j++) {
1680
1681 mcbsmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 1)) * setWCBsmmEW(j, xt, NLO_QED21), NLO_QED21);
1682 }
1683 case NLO_QED02: /*absent at high energy */
1684 for (int j = 0; j < 8; j++) {
1685
1686 mcbsmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 1)) * 0., NLO_QED02);
1687 }
1688 case NLO_QED11:
1689 for (int j = 0; j < 8; j++) {
1690
1691 mcbsmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 1)) * setWCBsmmEW(j, xt, NLO_QED11), NLO_QED11);
1692
1693 }
1694 case LO_QED: /*absent at high energy */
1695 for (int j = 0; j < 8; j++) {
1696
1697 mcbsmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 1)) * 0., LO_QED);
1698 }
1699 break;
1700 default:
1701 std::stringstream out;
1702 out << mcbsmm.getOrder_qed();
1703 throw std::runtime_error("StandardModelMatching::CMbsmm(): order_qed " + out.str() + "not implemented");
1704 }
1705 vmcbsmm.push_back(mcbsmm);
1706 return (vmcbsmm);
1707
1708}
1709
1710/*******************************************************************************
1711 * Wilson coefficients Misiak base for bd -> mu mu. decays *
1712 * operator basis: - current current *
1713 * - qcd penguins *
1714 * - semileptonic *
1715 * ****************************************************************************/
1716
1717std::vector<WilsonCoefficient>& StandardModelMatching::CMbdmm() {
1718
1719 // The couplings are not used here, but in the Bdmumu class.
1720
1721 double xt = x_t(Muw);
1722
1723 vmcbdmm.clear();
1724
1725
1726
1727 switch (mcbdmm.getScheme()) {
1728 case NDR:
1729 //case HV:
1730 //case LRI:
1731 break;
1732 default:
1733 std::stringstream out;
1734 out << mcbdmm.getScheme();
1735 throw std::runtime_error("StandardModelMatching::CMbdmm(): scheme " + out.str() + "not implemented");
1736 }
1737
1738 mcbdmm.setMu(Muw);
1739
1740 switch (mcbdmm.getOrder()) {
1741 case NNLO:
1742
1743 for (int j = 0; j < 8; j++) {
1744
1745 mcbdmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 0)) *
1746 setWCBdmm(j, xt, NNLO), NNLO);
1747 }
1748 case NLO:
1749
1750 for (int j = 0; j < 8; j++) {
1751
1752 mcbdmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 0)) *
1753 setWCBdmm(j, xt, NLO), NLO);
1754 }
1755 case LO:
1756
1757 for (int j = 0; j < 8; j++) {
1758
1759 mcbdmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 0)) * setWCBdmm(j, xt, LO), LO);
1760 }
1761 break;
1762 default:
1763 std::stringstream out;
1764 out << mcbdmm.getOrder();
1765 throw std::runtime_error("StandardModelMatching::CMbdmm(): order " + out.str() + "not implemented");
1766 }
1767
1768 switch (mcbdmm.getOrder_qed()) {
1769 case NLO_QED22:
1770 ;
1771 for (int j = 0; j < 8; j++) {
1772
1773 mcbdmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 0)) * setWCBdmmEW(j, xt, NLO_QED22), NLO_QED22);
1774
1775 }
1776 case NLO_QED12: /*absent at high energy */
1777 for (int j = 0; j < 8; j++) {
1778
1779 mcbdmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 0)) * 0., NLO_QED12);
1780 }
1781 case NLO_QED21:
1782
1783 for (int j = 0; j < 8; j++) {
1784
1785 mcbdmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 0)) * setWCBdmmEW(j, xt, NLO_QED21), NLO_QED21);
1786 }
1787 case NLO_QED02: /*absent at high energy */
1788 for (int j = 0; j < 8; j++) {
1789
1790 mcbdmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 0)) * 0., NLO_QED02);
1791 }
1792 case NLO_QED11:
1793 for (int j = 0; j < 8; j++) {
1794
1795 mcbdmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 0)) * setWCBdmmEW(j, xt, NLO_QED11), NLO_QED11);
1796
1797 }
1798 case LO_QED: /*absent at high energy */
1799 for (int j = 0; j < 8; j++) {
1800
1801 mcbdmm.setCoeff(j, (Vckm(2, 2).conjugate() * Vckm(2, 0)) * 0., LO_QED);
1802 }
1803 break;
1804 default:
1805 std::stringstream out;
1806 out << mcbdmm.getOrder_qed();
1807 throw std::runtime_error("StandardModelMatching::CMbdmm(): order_qed " + out.str() + "not implemented");
1808 }
1809 vmcbdmm.push_back(mcbdmm);
1810 return (vmcbdmm);
1811
1812}
1813
1814/*******************************************************************************
1815 * Wilson coefficients matching, LEFT basis [1709.04486] ordered as CnueduVLLkkij^+, CnueduVLRkkij^+, CnueduSRRkkij^+, CnueduSRLkkij^+, CnueduTRRkkij^+ *
1816 * ****************************************************************************/
1817std::vector<WilsonCoefficient>& StandardModelMatching::CMdiujleptonknu(int i, int j, int k) {
1818
1819 vmculeptonnu.clear();
1820
1821 mculeptonnu.setMu(Muw);
1822
1823 switch (mculeptonnu.getOrder()) {
1824 case NNLO:
1825 case NLO:
1826 case LO:
1827 mculeptonnu.setCoeff(0, 4. * GF * Vckm(j, i) / sqrt(2.), LO);
1828 break;
1829 default:
1830 std::stringstream out;
1831 out << mculeptonnu.getOrder();
1832 throw std::runtime_error("StandardModelMatching::CMuleptonnu(): order " + out.str() + "not implemented");
1833 }
1834
1835 vmculeptonnu.push_back(mculeptonnu);
1836 return (vmculeptonnu);
1837
1838}
1839
1840/******************************************************************************/
1841
1842/*******************************************************************************
1843 * Wilson coefficients Buras base for b -> nonlep. decays *
1844 * operator basis: - current current *
1845 * - qcd penguins *
1846 * - magnetic and chromomagnetic penguins *
1847 * - semileptonic *
1848 * i=0 deltaS=0 deltaC=0; i=1 1,0 ; *
1849 * ****************************************************************************/
1850std::vector<WilsonCoefficient>& StandardModelMatching::CMbnlep(const int a) {
1851 gslpp::complex lambda;
1852
1853 switch (a) {
1854 case 0: lambda = SM.getCKM().computelamt_d();
1855 break;
1856 case 1: lambda = SM.getCKM().computelamt_s();
1857 break;
1858 default:
1859 std::stringstream out;
1860 out << a;
1861 throw std::runtime_error("case" + out.str() + "not implemented; implemented i=0,1,2,3");
1862 }
1863
1864 double xt = x_t(Muw);
1865 double co = (GF / sqrt(2));
1866
1867 vmcbnlep.clear();
1868
1869 switch (mcbnlep.getScheme()) {
1870 case NDR:
1871 //case HV:
1872 //case LRI:
1873 break;
1874 default:
1875 std::stringstream out;
1876 out << mcbnlep.getScheme();
1877 throw std::runtime_error("StandardModelMatching::CMbnlep(): scheme " + out.str() + "not implemented");
1878 }
1879
1880 mcbnlep.setMu(Muw);
1881
1882 switch (mcbnlep.getOrder()) {
1883 case NNLO:
1884 case NLO:
1885 for (int j = 0; j < 10; j++) {
1886 mcbnlep.setCoeff(j, co * lambda * SM.Als(Muw, FULLNLO) / 4. / M_PI * //* CHECK ORDER *//
1887 setWCbnlep(j, xt, NLO), NLO);
1888 mcbnlep.setCoeff(j, co * lambda * Ale / 4. / M_PI *
1889 setWCbnlepEW(j, xt), NLO_QED11);
1890 }
1891 case LO:
1892 for (int j = 0; j < 10; j++) {
1893 mcbnlep.setCoeff(j, co * lambda * setWCbnlep(j, xt, LO), LO);
1894 mcbnlep.setCoeff(j, 0., LO_QED);
1895 }
1896 break;
1897 default:
1898 std::stringstream out;
1899 out << mcbnlep.getOrder();
1900 throw std::runtime_error("StandardModelMatching::CMbnlep(): order " + out.str() + "not implemented");
1901 }
1902
1903 vmcbnlep.push_back(mcbnlep);
1904 return (vmcbnlep);
1905}
1906
1907/*******************************************************************************
1908 * Wilson coefficients Buras base for b -> nonlep. decays *
1909 * operator basis: - current current opertors *
1910 * i=0 deltaS=0 deltaC=0; i=1 1,0 ; i=2 0,1 ; i=3 1,1 *
1911 * ****************************************************************************/
1912std::vector<WilsonCoefficient>& StandardModelMatching::CMbnlepCC(const int a) {
1913 gslpp::complex lambda1 = 0.;
1914 //gslpp::complex lambda2 = 0.;
1915 //gslpp::matrix<gslpp::complex> ckm = SM.getVCKM();
1916
1917 switch (a) {
1918 case 0: lambda1 = SM.getCKM().computelamu_d();
1919 break;
1920 case 1: lambda1 = SM.getCKM().computelamu_s();
1921 break;
1922 case 2: lambda1 = Vckm(0, 2).conjugate() * Vckm(1, 0);
1923 break;
1924 case 3: lambda1 = Vckm(1, 2).conjugate() * Vckm(0, 0);
1925 break;
1926 case 4: lambda1 = Vckm(0, 2).conjugate() * Vckm(1, 1);
1927 break;
1928 case 5: lambda1 = Vckm(1, 2).conjugate() * Vckm(2, 1);
1929 break;
1930 default:
1931 std::stringstream out;
1932 out << a;
1933 throw std::runtime_error("case" + out.str() + " not existing; implemented i=0,1,2,3");
1934 }
1935
1936 double xt = x_t(Muw);
1937 double co = (GF / sqrt(2));
1938
1939 vmcbnlepCC.clear();
1940
1941 switch (mcbnlepCC.getScheme()) {
1942 case NDR:
1943 //case HV:
1944 //case LRI:
1945 break;
1946 default:
1947 std::stringstream out;
1948 out << mcbnlepCC.getScheme();
1949 throw std::runtime_error("StandardModelMatching::CMbnlepCC(): scheme " + out.str() + "not implemented");
1950 }
1951
1952 mcbnlepCC.setMu(Muw);
1953
1954 switch (mcbnlepCC.getOrder()) {
1955 case NNLO:
1956 case NLO:
1957 for (int j = 0; j < 2; j++) {
1958 mcbnlepCC.setCoeff(j, co * lambda1 * setWCbnlep(j, xt, NLO), NLO);
1959 }
1960 for (int j = 2; j < 10; j++) {
1961 mcbnlepCC.setCoeff(j, 0., NLO);
1962 }
1963 case LO:
1964 for (int j = 0; j < 2; j++) {
1965 mcbnlepCC.setCoeff(j, co * lambda1 * setWCbnlep(j, xt, LO), LO);
1966 }
1967 for (int j = 2; j < 10; j++) {
1968 mcbnlepCC.setCoeff(j, 0., LO);
1969 }
1970 break;
1971 default:
1972 std::stringstream out;
1973 out << mcbnlepCC.getOrder();
1974 throw std::runtime_error("StandardModelMatching::CMbnlepCC(): order " + out.str() + "not implemented");
1975 }
1976
1977 vmcbnlepCC.push_back(mcbnlepCC);
1978 return (vmcbnlepCC);
1979}
1980
1981std::vector<WilsonCoefficient>& StandardModelMatching::CMkpnn() {
1982
1983 double co = 4. * GF / sqrt(2.) * SM.alphaMz() / 2. / M_PI / SM.sW2_ND() ; //SM prefactor as in eq. (1.1) of arXiv:1009.0947
1984
1985 vmckpnn.clear();
1986
1987 mckpnn.setMu(Mut);
1988 double xt = x_t(Mut,FULLNNLO);
1989
1990 switch (mckpnn.getOrder()) {
1991 case NNLO:
1992 case NLO:
1993 mckpnn.setCoeff(0, co * lam_t * SM.Als(Mut, FULLNLO) / 4. / M_PI * X1t(xt), NLO);
1994 case LO:
1995 mckpnn.setCoeff(0, co * lam_t * X0t(xt), LO);
1996 break;
1997 default:
1998 std::stringstream out;
1999 out << mckpnn.getOrder();
2000 throw std::runtime_error("StandardModelMatching::CMkpnn(): order " + out.str() + " not implemented");
2001 }
2002
2003 switch (mckpnn.getOrder_qed()) {
2004 case NLO_QED11:
2005 mckpnn.setCoeff(0, co * lam_t * SM.Ale(Mut,FULLNLO) / 4. / M_PI * Xewt(xt, SM.getMHl() * SM.getMHl() / getMt_mut() / getMt_mut(), Mut), NLO_QED11);
2006 case LO_QED:
2007 mckpnn.setCoeff(0, 0. , LO_QED);
2008 break;
2009 default:
2010 std::stringstream out;
2011 out << mckpnn.getOrder_qed();
2012 throw std::runtime_error("StandardModelMatching::CMkpnn(): qed order " + out.str() + " not implemented");
2013 }
2014
2015 vmckpnn.push_back(mckpnn);
2016 return (vmckpnn);
2017
2018}
2019
2020std::vector<WilsonCoefficient>& StandardModelMatching::CMkmm() {
2021
2022 //PROBLEMI: mu e sin(theta_weak) e la scala di als
2023
2024 double xt = x_t(Muw);
2025
2026 vmckmm.clear();
2027
2028 mckmm.setMu(Mut);
2029
2030 switch (mckmm.getOrder()) {
2031 case NNLO:
2032 case NLO:
2033 mckmm.setCoeff(0, SM.Als(Muw, FULLNLO) / 4. / M_PI * lam_t.real() * Y1(xt, Muw) / SM.getCKM().getLambda(), NLO); //* CHECK ORDER *//
2034 case LO:
2035 mckmm.setCoeff(0, lam_t.real() * Y0(xt) / SM.getCKM().getLambda(), LO);
2036 break;
2037 default:
2038 std::stringstream out;
2039 out << mckmm.getOrder();
2040 throw std::runtime_error("StandardModelMatching::CMkmm(): order " + out.str() + "not implemented");
2041 }
2042
2043 vmckmm.push_back(mckmm);
2044 return (vmckmm);
2045
2046}
2047
2048std::vector<WilsonCoefficient>& StandardModelMatching::CMBXsnn(QCD::lepton lepton) {
2049
2050 double xt = x_t(Mut,FULLNNLO);
2051 sw2 = SM.sW2_ND();
2052
2053 vmcbsnn.clear();
2054
2055 mcbsnn.setMu(Mut);
2056
2057 switch (mcbsnn.getOrder()) {
2058 case NNLO:
2059 case NLO:
2060 mcbsnn.setCoeff(0, -1/sw2 * SM.Als(Mut, FULLNLO) / 4. / M_PI * X1t(xt), NLO); //* CHECK ORDER *// //* CKM ELEMENTS IN CLASS AND NOT HERE *//
2061 mcbsnn.setCoeff(1, 0., NLO);
2062 case LO:
2063 mcbsnn.setCoeff(0, -1/sw2 * X0t(xt), LO);
2064 mcbsnn.setCoeff(1, 0., LO);
2065 break;
2066 default:
2067 std::stringstream out;
2068 out << mcbsnn.getOrder();
2069 throw std::runtime_error("StandardModelMatching::CMXsnn(): order " + out.str() + "not implemented");
2070 }
2071
2072 switch (mcbsnn.getOrder_qed()) {
2073 case NLO_QED11:
2074 mcbsnn.setCoeff(0, -1/sw2 * SM.Ale(Mut,FULLNLO) / 4. / M_PI * Xewt(xt, SM.getMHl() * SM.getMHl() / getMt_mut() / getMt_mut(), Mut), NLO_QED11);
2075 mcbsnn.setCoeff(1, 0., NLO_QED11);
2076 case LO_QED:
2077 mcbsnn.setCoeff(0, 0., LO_QED);
2078 mcbsnn.setCoeff(1, 0., LO_QED);
2079 break;
2080 default:
2081 std::stringstream out;
2082 out << mcbsnn.getOrder_qed();
2083 throw std::runtime_error("StandardModelMatching::CMXsnn(): qed order " + out.str() + " not implemented");
2084 }
2085
2086 vmcbsnn.push_back(mcbsnn);
2087 return (vmcbsnn);
2088
2089}
2090
2091std::vector<WilsonCoefficient>& StandardModelMatching::CMBXdnn() {
2092
2093 double xt = x_t(Mut,FULLNNLO);
2094 sw2 = SM.sW2_ND();
2095
2096 vmcbdnn.clear();
2097
2098 mcbdnn.setMu(Mut);
2099
2100 switch (mcbdnn.getOrder()) {
2101 case NNLO:
2102 case NLO:
2103 mcbdnn.setCoeff(0, -1/sw2 * SM.Als(Mut, FULLNLO) / 4. / M_PI * X1t(xt), NLO); //* CHECK ORDER *// //* CKM ELEMENTS IN CLASS AND NOT HERE *//
2104 case LO:
2105 mcbdnn.setCoeff(0, -1/sw2 * X0t(xt), LO);
2106 break;
2107 default:
2108 std::stringstream out;
2109 out << mcbdnn.getOrder();
2110 throw std::runtime_error("StandardModelMatching::CXdnn(): order " + out.str() + "not implemented");
2111 }
2112
2113 switch (mcbdnn.getOrder_qed()) {
2114 case NLO_QED11:
2115 mcbdnn.setCoeff(0, -1/sw2 * SM.Ale(Mut,FULLNLO) / 4. / M_PI * Xewt(xt, SM.getMHl() * SM.getMHl() / getMt_mut() / getMt_mut(), Mut), NLO_QED11);
2116 case LO_QED:
2117 mcbdnn.setCoeff(0, 0., LO_QED);
2118 break;
2119 default:
2120 std::stringstream out;
2121 out << mcbdnn.getOrder_qed();
2122 throw std::runtime_error("StandardModelMatching::CXdnn(): qed order " + out.str() + " not implemented");
2123 }
2124
2125 vmcbdnn.push_back(mcbdnn);
2126 return (vmcbdnn);
2127
2128}
2129
2130/*******************************************************************************
2131 * Wilson coefficients calculus, MISIAK base for Bs to mu mu decay *
2132 * ****************************************************************************/
2133double StandardModelMatching::setWCBsmm(int i, double x, orders order) {
2134
2135 sw2 = (M_PI * Ale) / (sqrt(2.) * GF * Mw * Mw);
2136
2137 if (sw2d == sw2 && xcached == x) {
2138 switch (order) {
2139 case NNLO:
2140 return (CWBsmmArrayNNLOqcd[i]);
2141 break;
2142 case NLO:
2143 return (CWBsmmArrayNLOqcd[i]);
2144 break;
2145 case LO:
2146 return (CWBsmmArrayLOqcd[i]);
2147 break;
2148 default:
2149 std::stringstream out;
2150 out << order;
2151 throw std::runtime_error("order" + out.str() + "not implemeted");
2152 }
2153 }
2154
2155 sw2d = sw2;
2156 xcached = x;
2157
2158 switch (order) {
2159 case NNLO:
2160 CWBsmmArrayNNLOqcd[0] = sw2 * (-Tt(x) + 7987. / 72. + 17. * M_PI * M_PI / 3. + 475. * L / 6. + 17. * L * L);
2161 CWBsmmArrayNNLOqcd[1] = sw2 * (127. / 18. + 4. * M_PI * M_PI / 3. + 46. * L / 3. + 4. * L * L);
2162 CWBsmmArrayNNLOqcd[2] = sw2 * (G1t(x, Muw) - 680. / 243. - 20. * M_PI * M_PI / 81. - 68. * L / 81. - 20. * L * L / 27.);
2163 CWBsmmArrayNNLOqcd[3] = sw2 * (E1t(x, Muw) + 950. / 243. + 10. * M_PI * M_PI / 81. + 124. * L / 27. + 10. * L * L / 27.);
2164 CWBsmmArrayNNLOqcd[4] = sw2 * (-G1t(x, Muw) / 10. + 2. * E0t(x) / 15. + 68. / 243. + 2. * M_PI * M_PI / 81. + 14. * L / 81. + 2. * L * L / 27.);
2165 CWBsmmArrayNNLOqcd[5] = sw2 * (-3. * G1t(x, Muw) / 16. + E0t(x) / 4. + 85. / 162. + 5. * M_PI * M_PI / 108. + 35. * L / 108. + 5. * L * L / 36.);
2166
2167 case NLO:
2168 CWBsmmArrayNLOqcd[0] = sw2 * (15. + 6. * L);
2169 CWBsmmArrayNLOqcd[3] = sw2 * (Eet(x) - 2. / 3. + 2. * L / 3.);
2170
2171 case LO:
2172 CWBsmmArrayLOqcd[1] = sw2 * 1.;
2173
2174 break;
2175 default:
2176 std::stringstream out;
2177 out << order;
2178 throw std::runtime_error("order" + out.str() + "not implemeted");
2179 }
2180 switch (order) {
2181 case NNLO:
2182 return (CWBsmmArrayNNLOqcd[i]);
2183
2184 break;
2185 case NLO:
2186 return (CWBsmmArrayNLOqcd[i]);
2187
2188 break;
2189 case LO:
2190 return (CWBsmmArrayLOqcd[i]);
2191
2192 break;
2193 default:
2194 std::stringstream out;
2195 out << order;
2196 throw std::runtime_error("order" + out.str() + "not implemeted");
2197 }
2198}
2199
2200double StandardModelMatching::setWCBsmmEW(int i, double x, orders_qed order_qed) {
2201 sw2 = (M_PI * Ale) / (sqrt(2.) * GF * Mw * Mw);
2202
2203 double mt = SM.Mrun(Muw, SM.getQuarks(QCD::TOP).getMass_scale(),
2204 SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, FULLNNLO);
2205
2206 if (sw2e == sw2 && xcachee == x) {
2207 switch (order_qed) {
2208 case NLO_QED22:
2209 return (CWBsmmArrayNLOewt4[i]);
2210 break;
2211 case NLO_QED21:
2212 return (CWBsmmArrayNLOewt2[i]);
2213 break;
2214 case NLO_QED11:
2215 return (CWBsmmArrayNLOew[i]);
2216 break;
2217 default:
2218 std::stringstream out;
2219 out << order_qed;
2220 throw std::runtime_error("order_qed" + out.str() + "not implemeted");
2221 }
2222
2223 }
2224
2225
2226 sw2e = sw2;
2227 xcachee = x;
2228
2231
2232 switch (order_qed) {
2233 case NLO_QED22:
2236 CWBsmmArrayNLOewt4[7] = sw2 * (1. / (sw2)) * Rest(x, Muw);
2237 // std::cout << ">>>>>>> " << Rest(163.5*163.5/80.358/80.358, 80.358) << std::endl;
2238
2239 case NLO_QED21:
2241 // 2. / 27. * B1u_tilde(x, Muw) - 2. / 9. * C1ew(x) + 320. / 27. * B0b(x) + 160. / 27. * C0b(x)));
2243 // 2. / 9. * Gew(x, xz, Muw) - 88. / 9. * B0b(x) - 184. / 27. * C0b(x)));
2245 // 1. / 54. * B1u_tilde(x, Muw) + 1. / 18. * C1ew(x) - 32. / 27. * B0b(x) - 16. / 27. * C0b(x)));
2247 // 4. / 3. * B0b(x) + 2. / 3. * C0b(x)));
2248 CWBsmmArrayNLOewt2[6] = sw2 * ((1. - 4. * sw2) * C1t(x, Muw) / (sw2) - B1t(x, Muw) / (sw2)
2249 - D1t(x, Muw) + 1. / (sw2) + 524. / 729. - 128. * M_PI * M_PI / 243.
2250 - 16. * L / 3. - 128. * L * L / 81.);
2251 CWBsmmArrayNLOewt2[7] = sw2 * ((1. / (sw2)) * (B1t(x, Muw) - C1t(x, Muw)) - 1. / (sw2));
2252
2253 case NLO_QED11:
2257 CWBsmmArrayNLOew[6] = sw2 * (Y0(x) / (sw2) + Wt(x) + 4. / 9. - 4. * 2 * log(Muw / mt) / 9.);
2258 CWBsmmArrayNLOew[7] = sw2 * (-Y0(x) / (sw2));
2259
2260 break;
2261 default:
2262 std::stringstream out;
2263 out << order_qed;
2264 throw std::runtime_error("order_qed" + out.str() + "not implemeted");
2265 }
2266
2267 switch (order_qed) {
2268 case NLO_QED22:
2269 return (CWBsmmArrayNLOewt4[i]);
2270 break;
2271 case NLO_QED21:
2272 return (CWBsmmArrayNLOewt2[i]);
2273 break;
2274 case NLO_QED11:
2275 return (CWBsmmArrayNLOew[i]);
2276 break;
2277 default:
2278 std::stringstream out;
2279 out << order_qed;
2280 throw std::runtime_error("order_qed" + out.str() + "not implemeted");
2281 }
2282}
2283
2284/*******************************************************************************
2285 * Wilson coefficients calculus, MISIAK base for Bd to mu mu decay *
2286 * ****************************************************************************/
2287double StandardModelMatching::setWCBdmm(int i, double x, orders order) {
2288
2289 sw2 = (M_PI * Ale) / (sqrt(2.) * GF * Mw * Mw);
2290
2291 if (sw2b == sw2 && xcacheb == x) {
2292 switch (order) {
2293 case NNLO:
2294 return (CWBdmmArrayNNLOqcd[i]);
2295 break;
2296 case NLO:
2297 return (CWBdmmArrayNLOqcd[i]);
2298 break;
2299 case LO:
2300 return (CWBdmmArrayLOqcd[i]);
2301 break;
2302 default:
2303 std::stringstream out;
2304 out << order;
2305 throw std::runtime_error("order" + out.str() + "not implemeted");
2306 }
2307 }
2308
2309 sw2b = sw2;
2310 xcacheb = x;
2311
2312 switch (order) {
2313 case NNLO:
2314 CWBdmmArrayNNLOqcd[0] = sw2 * (-Tt(x) + 7987. / 72. + 17. * M_PI * M_PI / 3. + 475. * L / 6. + 17. * L * L);
2315 CWBdmmArrayNNLOqcd[1] = sw2 * (127. / 18. + 4. * M_PI * M_PI / 3. + 46. * L / 3. + 4. * L * L);
2316 CWBdmmArrayNNLOqcd[2] = sw2 * (G1t(x, Muw) - 680. / 243. - 20. * M_PI * M_PI / 81. - 68. * L / 81. - 20. * L * L / 27.);
2317 CWBdmmArrayNNLOqcd[3] = sw2 * (E1t(x, Muw) + 950. / 243. + 10. * M_PI * M_PI / 81. + 124. * L / 27. + 10. * L * L / 27.);
2318 CWBdmmArrayNNLOqcd[4] = sw2 * (-G1t(x, Muw) / 10. + 2. * E0t(x) / 15. + 68. / 243. + 2. * M_PI * M_PI / 81. + 14. * L / 81. + 2. * L * L / 27.);
2319 CWBdmmArrayNNLOqcd[5] = sw2 * (-3. * G1t(x, Muw) / 16. + E0t(x) / 4. + 85. / 162. + 5. * M_PI * M_PI / 108. + 35. * L / 108. + 5. * L * L / 36.);
2320
2321 case NLO:
2322 CWBdmmArrayNLOqcd[0] = sw2 * (15. + 6. * L);
2323 CWBdmmArrayNLOqcd[3] = sw2 * (Eet(x) - 2. / 3. + 2. * L / 3.);
2324
2325 case LO:
2326 CWBdmmArrayLOqcd[1] = sw2 * 1.;
2327
2328 break;
2329 default:
2330 std::stringstream out;
2331 out << order;
2332 throw std::runtime_error("order" + out.str() + "not implemeted");
2333 }
2334 switch (order) {
2335 case NNLO:
2336 return (CWBdmmArrayNNLOqcd[i]);
2337
2338 break;
2339 case NLO:
2340 return (CWBdmmArrayNLOqcd[i]);
2341
2342 break;
2343 case LO:
2344 return (CWBdmmArrayLOqcd[i]);
2345
2346 break;
2347 default:
2348 std::stringstream out;
2349 out << order;
2350 throw std::runtime_error("order" + out.str() + "not implemeted");
2351 }
2352}
2353
2354double StandardModelMatching::setWCBdmmEW(int i, double x, orders_qed order_qed) {
2355 sw2 = (M_PI * Ale) / (sqrt(2.) * GF * Mw * Mw);
2356
2357 double mt = SM.Mrun(Muw, SM.getQuarks(QCD::TOP).getMass_scale(),
2358 SM.getQuarks(QCD::TOP).getMass(), QCD::TOP, FULLNNLO);
2359
2360 if (sw2c == sw2 && xcachec == x) {
2361 switch (order_qed) {
2362 case NLO_QED22:
2363 return (CWBdmmArrayNLOewt4[i]);
2364 break;
2365 case NLO_QED21:
2366 return (CWBdmmArrayNLOewt2[i]);
2367 break;
2368 case NLO_QED11:
2369 return (CWBdmmArrayNLOew[i]);
2370 break;
2371 default:
2372 std::stringstream out;
2373 out << order_qed;
2374 throw std::runtime_error("order_qed" + out.str() + "not implemeted");
2375 }
2376
2377 }
2378
2379
2380 sw2c = sw2;
2381 xcachec = x;
2382
2383 switch (order_qed) {
2384 case NLO_QED22:
2385 CWBdmmArrayNLOewt4[7] = sw2 * (1. / (sw2)) * Rest(x, Muw);
2386
2387 case NLO_QED21:
2388 CWBdmmArrayNLOewt2[6] = sw2 * ((1. - 4. * sw2) * C1t(x, Muw) / (sw2) - B1t(x, Muw) / (sw2)
2389 - D1t(x, Muw) + 1. / (sw2) + 524. / 729. - 128. * M_PI * M_PI / 243.
2390 - 16. * L / 3. - 128. * L * L / 81.);
2391 CWBdmmArrayNLOewt2[7] = sw2 * ((1. / (sw2)) * (B1t(x, Muw) - C1t(x, Muw)) - 1. / (sw2));
2392
2393 case NLO_QED11:
2394 CWBdmmArrayNLOew[6] = sw2 * (Y0(x) / (sw2) + Wt(x) + 4. / 9. - 4. * 2 * log(Muw / mt) / 9.);
2395 CWBdmmArrayNLOew[7] = sw2 * (-Y0(x) / (sw2));
2396
2397 break;
2398 default:
2399 std::stringstream out;
2400 out << order_qed;
2401 throw std::runtime_error("order_qed" + out.str() + "not implemeted");
2402 }
2403
2404 switch (order_qed) {
2405 case NLO_QED22:
2406 return (CWBdmmArrayNLOewt4[i]);
2407 break;
2408 case NLO_QED02:
2409 return (CWBdmmArrayNLOewt2[i]);
2410 break;
2411 case NLO_QED11:
2412 return (CWBdmmArrayNLOew[i]);
2413 break;
2414 default:
2415 std::stringstream out;
2416 out << order_qed;
2417 throw std::runtime_error("order_qed" + out.str() + "not implemeted");
2418 }
2419}
2420
2421/*******************************************************************************
2422 * Wilson coefficients calculus, Buras base for nonlep. b decays *
2423 * ****************************************************************************/
2424double StandardModelMatching::setWCbnlep(int i, double x, orders order) {
2425 if (xcacheb == x) {
2426 switch (order) {
2427 case NNLO:
2428 case NLO:
2429 return (CWbnlepArrayNLOqcd[i]);
2430 break;
2431 case LO:
2432 return (CWbnlepArrayLOqcd[i]);
2433 break;
2434 default:
2435 std::stringstream out;
2436 out << order;
2437 throw std::runtime_error("order" + out.str() + "not implemeted");
2438 }
2439 }
2440
2441 xcacheb = x;
2442
2443 switch (order) {
2444 case NNLO:
2445 case NLO:
2446 CWbnlepArrayNLOqcd[0] = 11. / 2.;
2447 CWbnlepArrayNLOqcd[1] = -11. / 6.;
2448 CWbnlepArrayNLOqcd[2] = -1. / 6. * (E0b(x) - 2. / 3.);
2449 CWbnlepArrayNLOqcd[3] = 0.5 * (E0b(x) - 2. / 3.);
2450 CWbnlepArrayNLOqcd[4] = -1. / 6. * (E0b(x) - 2. / 3.);
2451 CWbnlepArrayNLOqcd[5] = 0.5 * (E0b(x) - 2. / 3.);
2452 case LO:
2453 CWbnlepArrayLOqcd[1] = 1.;
2454 break;
2455 default:
2456 std::stringstream out;
2457 out << order;
2458 throw std::runtime_error("order" + out.str() + "not implemeted");
2459 }
2460
2461 switch (order) {
2462 case NNLO:
2463 case NLO:
2464 return (CWbnlepArrayNLOqcd[i]);
2465 break;
2466 case LO:
2467 return (CWbnlepArrayLOqcd[i]);
2468 break;
2469 default:
2470 std::stringstream out;
2471 out << order;
2472 throw std::runtime_error("order" + out.str() + "not implemeted");
2473 }
2474}
2475
2476double StandardModelMatching::setWCbnlepEW(int i, double x) {
2477 sw2 = (M_PI * Ale) / (sqrt(2) * GF * Mw * Mw);
2478
2479 if (sw2b == sw2 && xcacheb == x) {
2480 return (CWbnlepArrayNLOew[i]);
2481 }
2482
2483 sw2b = sw2;
2484 xcacheb = x;
2485
2486 CWbnlepArrayNLOew[1] = -35. / 18.;
2487 CWbnlepArrayNLOew[2] = 2. / (3. * sw2) * (2. * B0b(x) + C0b(x));
2488 CWbnlepArrayNLOew[6] = 2. / 3. * (4. * C0b(x) + D0b(x) - 4. / 9.);
2489 CWbnlepArrayNLOew[8] = 2. / 3. * (4. * C0b(x) + D0b(x) - 4. / 9. + (1. / (sw2)) *
2490 (10. * B0b(x) - 4. * C0b(x)));
2491
2492 return (CWbnlepArrayNLOew[i]);
2493}
2494
2495gslpp::complex StandardModelMatching::S0c() const {
2496 double xc = x_c(SM.getMuc());
2497 gslpp::complex co = GF / 2. / M_PI * Mw * SM.getCKM().computelamc().conjugate();
2498#if SUSYFIT_DEBUG & 2
2499 std::cout << "im lambdac = " << (SM.getCKM().computelamc() * SM.getCKM().computelamc()).imag() << std::endl;
2500#endif
2501 return (co * co * S0(xc, xc));
2502}
2503
2504gslpp::complex StandardModelMatching::S0ct() const {
2505 double xc = SM.Mrun4(SM.getMuc(), SM.getQuarks(QCD::CHARM).getMass_scale(), SM.getQuarks(QCD::CHARM).getMass()) / Mw;
2506 xc *= xc;
2507 double xt = SM.Mrun4(SM.getMuw(), SM.getQuarks(QCD::TOP).getMass_scale(), SM.getQuarks(QCD::TOP).getMass()) / Mw;
2508 xt *= xt;
2509 double co = GF / 2. / M_PI * Mw;
2510#if SUSYFIT_DEBUG & 2
2511 std::cout << "im lamc lamt = " << (SM.getCKM().computelamc() * SM.getCKM().computelamt()).imag() << std::endl;
2512#endif
2513
2514 return ( co * co * 2. * SM.getCKM().computelamc().conjugate() * lam_t.conjugate() * S0(xc, xt));
2515}
2516
2517gslpp::complex StandardModelMatching::S0tt() const {
2518 double xt = x_t(Mut);
2519 gslpp::complex co = GF / 2. / M_PI * Mw * lam_t.conjugate();
2520#if SUSYFIT_DEBUG & 2
2521 double pino = SM.Mrun(Mut, SM.Mp2Mbar(SM.getMtpole(), FULLNNLO),
2522 SM.Mp2Mbar(SM.getMtpole(), FULLNNLO), FULLNNLO);
2523 std::cout << "mt(" << Mut << ")" << pino << std::endl;
2524 double poldo = pino * pino / SM.Mw() / SM.Mw();
2525 std::cout << "S0(" << poldo << ") = " << S0(poldo, poldo) << std::endl;
2526 std::cout << "S0(" << xt << ") = " << S0(xt, xt) << std::endl;
2527 std::cout << "im lamt = " << (SM.getCKM().computelamt() * SM.getCKM().computelamt()).imag() << std::endl;
2528#endif
2529
2530 return ( co * co * S0(xt, xt));
2531}
2532
2533/*******************************************************************************
2534 * Wilson coefficients for Lepton Flavour Violation *
2535 * ****************************************************************************/
2536
2537std::vector<WilsonCoefficient>& StandardModelMatching::CMDLij(int li_lj) {
2538
2539 vmcDLij.clear();
2540
2541 mcDLij.setMu(Muw);
2542
2543 switch (mcDLij.getOrder()) {
2544 case NNLO:
2545 case NLO:
2546 case LO:
2547 mcDLij.setCoeff(0, 0., LO);
2548 mcDLij.setCoeff(1, 0., LO);
2549 break;
2550 default:
2551 std::stringstream out;
2552 out << mcDLij.getOrder();
2553 throw std::runtime_error("StandardModelMatching::CMDLij(): order " + out.str() + " not implemented.\nFor lepton flavour violating observables only Leading Order (LO) necessary.");
2554 }
2555
2556 vmcDLij.push_back(mcDLij);
2557 return (vmcDLij);
2558
2559}
2560
2561std::vector<WilsonCoefficient>& StandardModelMatching::CMDLi3j(int li_lj) {
2562
2563 vmcDLi3j.clear();
2564
2565 mcDLi3j.setMu(Muw);
2566
2567 switch (mcDLi3j.getOrder()) {
2568 case NNLO:
2569 case NLO:
2570 case LO:
2571 mcDLi3j.setCoeff(0, 0., LO);
2572 mcDLi3j.setCoeff(1, 0., LO);
2573 mcDLi3j.setCoeff(2, 0., LO);
2574 mcDLi3j.setCoeff(3, 0., LO);
2575 mcDLi3j.setCoeff(4, 0., LO);
2576 mcDLi3j.setCoeff(5, 0., LO);
2577 mcDLi3j.setCoeff(6, 0., LO);
2578 mcDLi3j.setCoeff(7, 0., LO);
2579 mcDLi3j.setCoeff(8, 0., LO);
2580 mcDLi3j.setCoeff(9, 0., LO);
2581 mcDLi3j.setCoeff(10, 0., LO);
2582 mcDLi3j.setCoeff(11, 0., LO);
2583 mcDLi3j.setCoeff(12, 0., LO);
2584 mcDLi3j.setCoeff(13, 0., LO);
2585 mcDLi3j.setCoeff(14, 0., LO);
2586 mcDLi3j.setCoeff(15, 0., LO);
2587 mcDLi3j.setCoeff(16, 0., LO);
2588 mcDLi3j.setCoeff(17, 0., LO);
2589 mcDLi3j.setCoeff(18, 0., LO);
2590 mcDLi3j.setCoeff(19, 0., LO);
2591 break;
2592 default:
2593 std::stringstream out;
2594 out << mcDLi3j.getOrder();
2595 throw std::runtime_error("StandardModelMatching::CMDLi3j(): order " + out.str() + " not implemented.\nFor lepton flavour violating observables only Leading Order (LO) necessary.");
2596 }
2597
2598 vmcDLi3j.push_back(mcDLi3j);
2599 return (vmcDLi3j);
2600
2601}
2602
2603std::vector<WilsonCoefficient>& StandardModelMatching::CMmueconv() {
2604
2605 vmcmueconv.clear();
2606
2607 mcmueconv.setMu(Muw);
2608
2609 switch (mcmueconv.getOrder()) {
2610 case NNLO:
2611 case NLO:
2612 case LO:
2613 mcmueconv.setCoeff(0, 0., LO);
2614 mcmueconv.setCoeff(1, 0., LO);
2615 mcmueconv.setCoeff(2, 0., LO);
2616 mcmueconv.setCoeff(3, 0., LO);
2617 mcmueconv.setCoeff(4, 0., LO);
2618 mcmueconv.setCoeff(5, 0., LO);
2619 mcmueconv.setCoeff(6, 0., LO);
2620 mcmueconv.setCoeff(7, 0., LO);
2621 break;
2622 default:
2623 std::stringstream out;
2624 out << mcmueconv.getOrder();
2625 throw std::runtime_error("StandardModelMatching::CMmueconv(): order " + out.str() + " not implemented.\nFor lepton flavour violating observables only Leading Order (LO) necessary.");
2626 }
2627
2628 vmcmueconv.push_back(mcmueconv);
2629 return (vmcmueconv);
2630
2631}
2632
2633std::vector<WilsonCoefficient>& StandardModelMatching::CMgminus2mu() {
2634
2635 vmcgminus2mu.clear();
2636
2637 mcgminus2mu.setMu(Muw);
2638
2639 switch (mcgminus2mu.getOrder()) {
2640 case NNLO:
2641 case NLO:
2642 mcgminus2mu.setCoeff(0, 0., NLO);
2643 mcgminus2mu.setCoeff(1, 0., NLO);
2644 break;
2645 case LO:
2646 mcgminus2mu.setCoeff(0, 0., LO);
2647 mcgminus2mu.setCoeff(1, 0., LO);
2648 break;
2649 default:
2650 std::stringstream out;
2651 out << mcgminus2mu.getOrder();
2652 throw std::runtime_error("StandardModelMatching::CMgminus2mu(): order " + out.str() + " not implemented.");
2653 }
2654
2655 vmcgminus2mu.push_back(mcgminus2mu);
2656 return (vmcgminus2mu);
2657
2658}
2659
2660WilsonCoefficientNew& StandardModelMatching::mc_C() {
2661 switch (mcC.getScheme()) {
2662 case NDR:
2663 //case HV:
2664 //case LRI:
2665 break;
2666 default:
2667 std::stringstream out;
2668 out << mcC.getScheme();
2669 throw "StandardModel::mc_C(): scheme " + out.str() + "not implemented";
2670 }
2671
2672 mcC.setMu(Muw); // cleared too
2673
2674 switch (mcC.getOrder_QED()) {
2675 case QED2:
2676 case QED1:
2677 switch (mcC.getOrder_QCD()) {
2678 case QCD2:
2679 case QCD1:
2680 mcC.setCoeff(1, aletilde * (-22. / 9. - 4. / 3. * Lz + 1. / 9.), QCD1, QED1);
2681 case QCD0:
2682 break;
2683 default:
2684 std::stringstream out;
2685 out << mcC.getOrder_QCD();
2686 throw "StandardModelMatching::mc_C(): order " + out.str() + " not implemented";
2687 }
2688 case QED0:
2689 switch (mcC.getOrder_QCD()) {
2690 case QCD2:
2691 mcC.setCoeff(0, alstilde * alstilde * (-Tt(x_t(Muw)) + 7987. / 72. + 17. / 3. * M_PI * M_PI +
2692 L * 475. / 6. + 17. * L * L), QCD2, QED0);
2693 mcC.setCoeff(1, alstilde * alstilde * (127. / 18. + 4. / 3. * M_PI * M_PI + 46. / 3. * L + 4. * L * L),
2694 QCD2);
2695 case QCD1:
2696 mcC.setCoeff(0, alstilde * (15. + 6. * L), QCD1, QED0);
2697 case QCD0:
2698 mcC.setCoeff(1, 1., QCD0, QED0);
2699 break;
2700 default:
2701 std::stringstream out;
2702 out << mcC.getOrder_QCD();
2703 throw "StandardModelMatching::mc_C(): order " + out.str() + " not implemented";
2704 }
2705 break;
2706 default:
2707 std::stringstream out;
2708 out << mcC.getOrder_QED();
2709 throw "StandardModelMatching::mc_C(): order " + out.str() + "not implemented";
2710 }
2711
2712 return (mcC);
2713}
2714
2715double StandardModelMatching::C3funNNLO(double x) {
2716 return (G1t(x, Muw) - 680. / 243. - 20. / 81. * M_PI * M_PI - 68. / 81. * L - 20. / 27 * L * L);
2717}
2718
2719double StandardModelMatching::C4fun(double x, orders ord) {
2720 switch (ord) {
2721 case NNLO:
2722 return (E1t(x, Muw) + 950. / 243. + 10. / 81. * M_PI * M_PI + 124. / 27. * L + 10. / 27. * L * L);
2723 case NLO:
2724 return (E0t(x) - 7. / 9. + 2. / 3. * L);
2725 default:
2726 std::stringstream out;
2727 out << ord;
2728 throw "StandardModelMatching::C4fun(): order " + out.str() + "not implemented";
2729
2730 }
2731}
2732
2733double StandardModelMatching::C5funNNLO(double x) {
2734 return (-0.1 * G1t(x, Muw) + 2. / 15. * E0t(x) + 68. / 243. + 2. / 81. * M_PI * M_PI + 14. / 81. * L + 2. / 27. * L * L);
2735}
2736
2737double StandardModelMatching::C6funNNLO(double x) {
2738 return (-3. / 16. * G1t(x, Muw) + 0.25 * E0t(x) + 85. / 162. + 5. / 108. * M_PI * M_PI + 35. / 108. * L + 5. / 36 * L * L);
2739}
2740
2741WilsonCoefficientNew& StandardModelMatching::mc_P() {
2742 // double sW2 = (M_PI * SM.getAle() ) / ( sqrt(2.) * SM.getGF() * SM.Mw() * SM.Mw() ); // ******* FOR TEST *********
2743 double xt = x_t(Muw);
2744 double xz = SM.getMz() * SM.getMz() / Mw / Mw;
2745
2746 switch (mcP.getScheme()) {
2747 case NDR:
2748 //case HV:
2749 //case LRI:
2750 break;
2751 default:
2752 std::stringstream out;
2753 out << mcP.getScheme();
2754 throw "StandardModelMatching::mc_P(): scheme " + out.str() + "not implemented";
2755 }
2756
2757 mcP.setMu(Muw); // cleared too
2758
2759 switch (mcP.getOrder_QED()) {
2760 case QED2:
2761 case QED1:
2762 switch (mcP.getOrder_QCD()) {
2763 case QCD2:
2764 // mcP.setCoeff(0, aletilde * alstilde * (1. / sW2 * (4. / 9. * B1d(xt, Muw) + 4. / 27. * B1d_tilde(xt, Muw) + 2. / 9. * B1u(xt, Muw) +
2765 // 2. / 27. * B1u_tilde(xt, Muw) - 2. / 9. * C1ew(xt) + 320. / 27. * B0b(xt) +
2766 // 160. / 27. * C0b(xt))), NLO_QED21);
2767 // mcP.setCoeff(1, aletilde * alstilde * (16. / 27. * C0b(xt) + 1. / sW2 * (8. / 9. * B1d_tilde(xt, Muw) + 4. / 9. * B1u_tilde(xt, Muw) -
2768 // 2. / 9. * Gew(xt, xz, Muw) - 88. / 9. * B0b(xt) - 184. / 27. * C0b(xt))), NLO_QED21);
2769 // mcP.setCoeff(2, aletilde * alstilde * (1. / sW2 * (-1. / 9. * B1d(xt, Muw) - 1. / 27. * B1d_tilde(xt, Muw) - 1. / 18. * B1u(xt, Muw) -
2770 // 1. / 54. * B1u_tilde(xt, Muw) + 1. / 18. * C1ew(xt) - 32. / 27. * B0b(xt) -
2771 // 16. / 27. * C0b(xt))), NLO_QED21);
2772 // mcP.setCoeff(3, aletilde * alstilde * (1. / sW2 * (-2. / 9. * B1d_tilde(xt, Muw) - 1. / 9. * B1u_tilde(xt, Muw) + 1. / 18. * Gew(xt, xz, Muw) +
2773 // 4. / 3. * B0b(xt) + 2. / 3. * C0b(xt))), NLO_QED21);
2774
2775 // Expressions obtained from Haisch, Misiak notes (with corrected factors on Buras basis)
2776 mcP.setCoeff(0, aletilde * alstilde * (1. / sW2 * (
2777 4. / 9. * B1d(xt, Muw) + 4. / 27. * B1d_tilde(xt, Muw) +
2778 2. / 9. * B1u(xt, Muw) + 2. / 27. * B1u_tilde(xt, Muw) -
2779 2. / 9. * C1ew(xt) + 320. / 27. * B0b(xt) + 160. / 27. * C0b(xt))), QCD2, QED1);
2780 mcP.setCoeff(1, aletilde * alstilde * (1. / sW2 * (
2781 8. / 9. * B1d_tilde(xt, Muw) + 4. / 9. * B1u_tilde(xt, Muw) -
2782 2. / 9. * Gew(xt, xz, Muw) - 88. / 9. * B0b(xt) - 56. / 27. * C0b(xt))), QCD2, QED1);
2783 mcP.setCoeff(2, aletilde * alstilde * (1. / sW2 * (
2784 -1. / 9. * B1d(xt, Muw) - 1. / 27. * B1d_tilde(xt, Muw) -
2785 1. / 18. * B1u(xt, Muw) - 1. / 54. * B1u_tilde(xt, Muw) +
2786 1. / 18. * C1ew(xt) - 32. / 27. * B0b(xt) - 16. / 27. * C0b(xt))), QCD2, QED1);
2787 mcP.setCoeff(3, aletilde * alstilde * (1. / sW2 * (
2788 -2. / 9. * B1d_tilde(xt, Muw) - 1. / 9. * B1u_tilde(xt, Muw) +
2789 1. / 18. * Gew(xt, xz, Muw) + 4. / 3. * B0b(xt) + 2. / 3. * C0b(xt))), QCD2, QED1);
2790 case QCD1:
2791 mcP.setCoeff(0, aletilde * (-2. / 9. / sW2 * (2. * B0b(xt) + C0b(xt))), QCD1, QED1);
2792 // mcP.setCoeff(0, aletilde*(-2./9./sW2*(2.*Y0(xt) - X0t(xt))), NLO, QED1);
2793 mcP.setCoeff(2, aletilde * (1. / 9. / sW2 * (B0b(xt) + 1. / 2. * C0b(xt))), QCD1, QED1);
2794 // mcP.setCoeff(2, aletilde*(1./9./sW2*(Y0(xt) - 1./2.*X0t(xt))), NLO, QED1);
2795 case QCD0:
2796 break;
2797 default:
2798 std::stringstream out;
2799 out << mcP.getOrder_QCD();
2800 throw "StandardModelMatching::mc_P(): order " + out.str() + "not implemented";
2801 }
2802 case QED0:
2803 switch (mcP.getOrder_QCD()) {
2804 case QCD2:
2805 mcP.setCoeff(0, alstilde * alstilde * C3funNNLO(xt), QCD2, QED0);
2806 mcP.setCoeff(1, alstilde * alstilde * C4fun(xt, NNLO), QCD2, QED0);
2807 mcP.setCoeff(2, alstilde * alstilde * C5funNNLO(xt), QCD2, QED0);
2808 mcP.setCoeff(3, alstilde * alstilde * C6funNNLO(xt), QCD2, QED0);
2809 case QCD1:
2810 mcP.setCoeff(1, alstilde * C4fun(xt, NLO), QCD1, QED0);
2811 case QCD0:
2812 break;
2813 default:
2814 std::stringstream out;
2815 out << mcP.getOrder_QCD();
2816 throw "StandardModelMatching::mc_P(): order " + out.str() + "not implemented";
2817 }
2818 break;
2819 default:
2820 std::stringstream out;
2821 out << mcP.getOrder_QED();
2822 throw "StandardModelMatching::mc_P(): order " + out.str() + "not implemented";
2823 }
2824
2825 return (mcP);
2826}
2827
2828double StandardModelMatching::C7funLO(double x) {
2829 return (-0.5 * A0t(x) - 23. / 36.);
2830}
2831
2832double StandardModelMatching::C8funLO(double x) {
2833 return (-0.5 * F0t(x) - 1. / 3.);
2834}
2835
2836WilsonCoefficientNew& StandardModelMatching::mc_M() {
2837 double xt = x_t(Muw);
2838 double mH = SM.getMHl();
2839
2840 switch (mcM.getScheme()) {
2841 case NDR:
2842 //case HV:
2843 //case LRI:
2844 break;
2845 default:
2846 std::stringstream out;
2847 out << mcM.getScheme();
2848 throw "StandardModelMatching::mc_M(): scheme " + out.str() + "not implemented";
2849 }
2850
2851 mcM.setMu(Muw); // cleared too
2852
2853 switch (mcM.getOrder_QED()) {
2854 case QED2:
2855 case QED1:
2856 switch (mcM.getOrder_QCD()) {
2857 case QCD2:
2858 case QCD1:
2859 mcM.setCoeff(0, aletilde * (1. / sW2 * (1.11 - 1.15 * (1. - Mt_muw * Mt_muw / 170. / 170.) - 0.444 * log(mH / 100.) -
2860 0.21 * log(mH / 100.) * log(mH / 100.) - 0.513 * log(mH / 100.) * log(Mt_muw / 170.)) +
2861 (8. / 9. * C7funLO(xt) - 104. / 243.) * L), QCD1, QED1);
2862 mcM.setCoeff(1, aletilde * (1. / sW2 * (-0.143 + 0.156 * (1. - Mt_muw * Mt_muw / 170. / 170.) - 0.129 * log(mH / 100.) -
2863 0.0244 * log(mH / 100.) * log(mH / 100.) - 0.037 * log(mH / 100.) * log(Mt_muw / 170.)) +
2864 (4. / 9. * C8funLO(xt) - 4. / 3. * C7funLO(xt) - 58. / 81.) * L), QCD1, QED1);
2865 case QCD0:
2866 break;
2867 default:
2868 std::stringstream out;
2869 out << mcM.getOrder_QCD();
2870 throw "StandardModelMatching::mc_M(): order " + out.str() + "not implemented";
2871 }
2872
2873 case QED0:
2874 switch (mcM.getOrder_QCD()) {
2875 case QCD2:
2876 mcM.setCoeff(0, alstilde * alstilde * (C7t_3L_at_mt(xt) + C7t_3L_func(xt, Muw)-(C7c_3L_at_mW(xt) + 13763. / 2187. * L + 814. / 729. * L * L)
2877 - 1. / 3. * C3funNNLO(xt) - 4. / 9. * C4fun(xt, NNLO) - 20. / 3. * C5funNNLO(xt) - 80. / 9. * C6funNNLO(xt)), QCD2, QED0);
2878 mcM.setCoeff(1, alstilde * alstilde * (C8t_3L_at_mt(xt) + C8t_3L_func(xt, Muw)-(C8c_3L_at_mW(xt) + 16607. / 5832. * L + 397. / 486. * L * L)
2879 + C3funNNLO(xt) - 1. / 6. * C4fun(xt, NNLO) - 20. * C5funNNLO(xt) - 10. / 3. * C6funNNLO(xt)), QCD2, QED0);
2880 case QCD1:
2881 mcM.setCoeff(0, alstilde * (-0.5 * A1t(xt, Muw) + 713. / 243. + 4. / 81. * L - 4. / 9. * C4fun(xt, NLO)), QCD1, QED0);
2882 mcM.setCoeff(1, alstilde * (-0.5 * F1t(xt, Muw) + 91. / 324. - 4. / 27. * L - 1. / 6. * C4fun(xt, NLO)), QCD1, QED0);
2883 case QCD0:
2884 mcM.setCoeff(0, C7funLO(xt), QCD0, QED0);
2885 mcM.setCoeff(1, C8funLO(xt), QCD0, QED0);
2886 break;
2887 default:
2888 std::stringstream out;
2889 out << mcM.getOrder_QCD();
2890 throw "StandardModelMatching::mc_M(): order " + out.str() + "not implemented";
2891 }
2892 break;
2893
2894 default:
2895 std::stringstream out;
2896 out << mcM.getOrder_QED();
2897 throw "StandardModelMatching::mc_M(): order " + out.str() + "not implemented";
2898 }
2899
2900 return (mcM);
2901}
2902
2903double StandardModelMatching::fbb(double x) // from hep-ph/9707243
2904{
2905 return (-0.0380386 - 2.24502 * x + 3.8472 * x * x - 7.80586 * x * x * x + 10.0763 * x * x * x * x - 6.9751 * x * x * x * x * x
2906 + 1.95163 * x * x * x * x * x * x - 0.00550335 * log(x)); // fitted between 0 and 1
2907}
2908
2909double StandardModelMatching::gbb(double x) // from hep-ph/9707243
2910{
2911 if (x > 4.)
2912 return (sqrt(x - 4.) * log((1. - sqrt(1. - 4. / x)) / (1. + sqrt(1. - 4. / x))));
2913 else if (x >= 0.)
2914 return (2. * sqrt(4. - x) * acos(sqrt(x / 4.)));
2915 else
2916 throw "StandardModelMatching::gbb(): defined for non-negative argument only.";
2917}
2918
2919double StandardModelMatching::taub2(double x) // from hep-ph/9707243
2920{
2921 double ll = log(x);
2922 return ( 9. - 13. / 4. * x - 2. * x * x - x / 4. * (19. + 6. * x) * ll - x * x / 4. * (7. - 6. * x) * ll * ll - (1. / 4. + 7. / 2. * x * x - 3. * x * x * x) * M_PI * M_PI / 6. +
2923 (x / 2. - 2.) * sqrt(x) * gbb(x)+(x - 1.)*(x - 1.)*(4. * x - 7. / 4.) * gslpp_special_functions::dilog(1. - x)-(x * x * x - 33. / 4. * x * x +
2924 18. * x - 7.) * fbb(x));
2925}
2926
2927double StandardModelMatching::C10_OS1(double x, double mu) {
2928 double Mw_2 = Mw * Mw;
2929 double mt_2 = x * Mw_2;
2930 double mt = sqrt(mt_2);
2931
2932 double a1 = 69354.0995830457;
2933 double b1 = 114.072536156018;
2934 double c1 = -1802.9029763021;
2935 double d1 = -0.6880489462364;
2936 double e1 = 11.7037717083906;
2937 double f1 = -1.422845251416;
2938 double g1 = 0.00856609262141;
2939 double h1 = 0.00005469961928;
2940 double a2 = 4144.6231891191;
2941 double b2 = 0.1706756075896;
2942 double c2 = -104.32434492409;
2943 double d2 = 0.00072836578577;
2944 double e2 = 0.65260451121855;
2945 double f2 = -0.0007261241908;
2946 double a3 = -218.96716792134;
2947 double b3 = -0.039653543101;
2948 double c3 = 5.57241334587797;
2949 double d3 = 2.8289915e-6;
2950 double e3 = -0.0351477930751;
2951 double f3 = 0.00048165797959;
2952
2953 return (a1 + b1 * mt + c1 * Mw + d1 * mt_2 + e1 * Mw_2 + f1 * mt * Mw + g1 * mt_2 * Mw + h1 * mt_2 * mt +
2954 (a2 + b2 * mt + c2 * Mw + d2 * mt_2 + e2 * Mw_2 + f2 * mt * Mw) * log(mu) +
2955 (a3 + b3 * mt + c3 * Mw + d3 * mt_2 + e3 * Mw_2 + f3 * mt * Mw) * log(mu) * log(mu));
2956}
2957
2958double StandardModelMatching::Delta_t(double mu, double x) // from hep-ph/9707243
2959{
2960 return (18. * log(mu / Mt_muw) + 11. - x / 2. + x * (x - 6.) / 2. * log(x)+(x - 4.) / 2. * sqrt(x) * gbb(x));
2961}
2962
2963WilsonCoefficientNew& StandardModelMatching::mc_L() {
2964 // double sW2 = (M_PI * SM.getAle() ) / ( sqrt(2.) * SM.getGF() * SM.Mw() * SM.Mw() ); // ******* FOR TEST *********
2965 double xt = x_t(Muw);
2966 // double xht = SM.getMHl() * SM.getMHl() / Mt_muw / Mt_muw;
2967
2968
2969 switch (mcL.getScheme()) {
2970 case NDR:
2971 //case HV:
2972 //case LRI:
2973 break;
2974 default:
2975 std::stringstream out;
2976 out << mcL.getScheme();
2977 throw "StandardModel::mc_L(): scheme " + out.str() + "not implemented";
2978 }
2979
2980 mcL.setMu(Muw); // cleared too
2981
2982 switch (mcL.getOrder_QED()) {
2983 case QED2:
2984 switch (mcL.getOrder_QCD()) {
2985 case QCD2:
2986 //Eqs. (32-33) of ref. Huber et al.
2987 //Delta_t and tau_b in hep-ph/9707243
2988 // mcL.setCoeff(0, aletilde*aletilde*(-xt*xt/32/sW2/sW2*(4.*sW2-1.)*(3.+taub2(xht)-Delta_t(Muw,xht))), NNLO, QED2);
2989 //mcL.setCoeff(1, aletilde * aletilde * (-xt * xt / 32. / sW2 / sW2 * (3. + taub2(xht) - Delta_t(Muw, xht))), QCD2, QED2);
2990 //mcL.setCoeff(0, (4. * sW2 - 1.) * (mcL.getCoeff(QCD2, QED2)(1)), QCD2, QED2);
2991 // mcL.setCoeff(1, aletilde*aletilde*Rest(xt, Muw)/sW2, NNLO, QED2);
2992 mcL.setCoeff(1, aletilde * aletilde * C10_OS1(xt, Muw), QCD2, QED2);
2993 case QCD1:
2994 case QCD0:
2995 break;
2996 default:
2997 std::stringstream out;
2998 out << mcL.getOrder_QCD();
2999 throw "StandardModelMatching::mc_L(): order " + out.str() + "not implemented";
3000 }
3001
3002 case QED1:
3003 switch (mcL.getOrder_QCD()) {
3004 case QCD2:
3005 mcL.setCoeff(0, aletilde * alstilde * ((1. - 4. * sW2) / sW2 * C1t(xt, Muw) - 1. / sW2 * B1t(xt, Muw) - D1t(xt, Muw) + 1. / sW2 +
3006 524. / 729. - 128. / 243. * M_PI * M_PI - 16. / 3. * L - 128. / 81. * L * L), QCD2, QED1);
3007 mcL.setCoeff(1, aletilde * alstilde * (1. / sW2 * (B1t(xt, Muw) - C1t(xt, Muw)) - 1. / sW2), QCD2, QED1);
3008 case QCD1:
3009 mcL.setCoeff(0, aletilde * (1. / sW2 * Y0(xt) + Wt(xt) + 4. / 9. - 4. / 9. * 2. * log(Muw / Mt_muw)), QCD1, QED1);
3010 mcL.setCoeff(1, aletilde * (-1. / sW2 * Y0(xt)), QCD1, QED1);
3011 case QCD0:
3012 break;
3013 default:
3014 std::stringstream out;
3015 out << mcL.getOrder_QCD();
3016 throw "StandardModelMatching::mc_L(): order " + out.str() + "not implemented";
3017 }
3018
3019 case QED0:
3020 break;
3021
3022 default:
3023 std::stringstream out;
3024 out << mcL.getOrder_QED();
3025 throw "StandardModelMatching::mc_L(): order " + out.str() + "not implemented";
3026
3027 }
3028
3029 return (mcL);
3030}
3031
3032WilsonCoefficientNew& StandardModelMatching::mc_Q() {
3033 double xt = x_t(Muw);
3034 double xz = SM.getMz() * SM.getMz() / Mw / Mw;
3035
3036 switch (mcQ.getScheme()) {
3037 case NDR:
3038 //case HV:
3039 //case LRI:
3040 break;
3041 default:
3042 std::stringstream out;
3043 out << mcQ.getScheme();
3044 throw "StandardModel::mc_Q(): scheme " + out.str() + "not implemented";
3045 }
3046
3047 mcQ.setMu(Muw); // cleared too
3048
3049 switch (mcQ.getOrder_QED()) {
3050 case QED2:
3051 case QED1:
3052 switch (mcQ.getOrder_QCD()) {
3053 case QCD2:
3054 // mcQ.setCoeff(0, aletilde*alstilde*(4.*C1ew(xt) + D1t(xt,Muw) + 320./9.*C0b(xt) +
3055 // 1./sW2*(-2./3.*B1d(xt,Muw) - 2./9.*B1d_tilde(xt,Muw) +
3056 // 2./3.*B1u(xt,Muw) + 2./9.*B1u_tilde(xt,Muw) +
3057 // 4./3.*C1ew(xt) + 800./9.*B0b(xt) - 640./9.*C0b(xt))), NLO_QED21);
3058 // mcQ.setCoeff(1, aletilde*alstilde*(-4./3.*Gew(xt,xz,Muw) - 16./3.*Hew(xt,xz,Muw) - 32.*C0b(xt) +
3059 // 1./sW2*(-4./3.*B1d_tilde(xt,Muw) + 4./3.*B1u_tilde(xt,Muw) +
3060 // 4./3.*Gew(xt,xz,Muw) - 80.*B0b(xt) + 112./3.*C0b(xt))), NLO_QED21);
3061 // mcQ.setCoeff(2, aletilde*alstilde*(-32./9.*C0b(xt) +
3062 // 1./sW2*(1./6.*B1d(xt,Muw) + 1./18.*B1d_tilde(xt,Muw) -
3063 // 1./6.*B1u(xt,Muw) - 1./18.*B1u_tilde(xt,Muw) -
3064 // 1./3.*C1ew(xt) - 80./9.*B0b(xt) + 64./9.*C0b(xt))), NLO_QED21);
3065 // mcQ.setCoeff(3, aletilde*alstilde*(1./3.*Gew(xt,xz,Muw) + 1./3.*Hew(xt,xz,Muw) + 4.*C0b(xt) +
3066 // 1./sW2*(1./3.*B1d_tilde(xt,Muw) - 1./3.*B1u_tilde(xt,Muw) -
3067 // 1./3.*Gew(xt,xz,Muw) +10.*B0b(xt) - 16./3.*C0b(xt))), NLO_QED21);
3068 // Expressions obtained from Haisch, Misiak notes (with corrected factors on Buras basis)
3069 mcQ.setCoeff(0, aletilde * alstilde * (4. * C1ew(xt) + D1t(xt, Muw) + 1. / sW2 * (
3070 -2. / 3. * B1d(xt, Muw) - 2. / 9. * B1d_tilde(xt, Muw) +
3071 2. / 3. * B1u(xt, Muw) + 2. / 9. * B1u_tilde(xt, Muw) +
3072 4. / 3. * C1ew(xt) + 800. / 9. * B0b(xt) - 320. / 9. * C0b(xt))), QCD2, QED1);
3073 mcQ.setCoeff(1, aletilde * alstilde * (-4. / 3. * Gew(xt, xz, Muw) - 16. / 3. * Hew(xt, xz, Muw) -
3074 80. / 3. * C0b(xt) + 1. / sW2 * (-4. / 3. * B1d_tilde(xt, Muw) + 4. / 3. * B1u_tilde(xt, Muw) +
3075 4. / 3. * Gew(xt, xz, Muw) - 80. * B0b(xt) + 32. * C0b(xt))), QCD2, QED1);
3076 mcQ.setCoeff(2, aletilde * alstilde * (1. / sW2 * (1. / 6. * B1d(xt, Muw) + 1. / 18. * B1d_tilde(xt, Muw) -
3077 1. / 6. * B1u(xt, Muw) - 1. / 18. * B1u_tilde(xt, Muw) -
3078 1. / 3. * C1ew(xt) - 80. / 9. * B0b(xt) + 32. / 9. * C0b(xt))), QCD2, QED1);
3079 mcQ.setCoeff(3, aletilde * alstilde * (1. / 3. * Gew(xt, xz, Muw) + 1. / 3. * Hew(xt, xz, Muw) +
3080 8. / 3. * C0b(xt) + 1. / sW2 * (1. / 3. * B1d_tilde(xt, Muw) - 1. / 3. * B1u_tilde(xt, Muw) -
3081 1. / 3. * Gew(xt, xz, Muw) + 10. * B0b(xt) - 4. * C0b(xt))), QCD2, QED1);
3082 case QCD1:
3083 mcQ.setCoeff(0, aletilde * (4. * C0b(xt) + D0b_tilde(xt) + 4. / 9. * L - 1. / sW2 * (10. / 3. * B0b(xt) - 4. / 3. * C0b(xt))), QCD1, QED1); // log from Misiak's notes
3084 mcQ.setCoeff(2, aletilde * (1. / sW2 * (5. / 6. * B0b(xt) - 1. / 3. * C0b(xt))), QCD1, QED1);
3085 case QCD0:
3086 break;
3087 default:
3088 std::stringstream out;
3089 out << mcQ.getOrder_QCD();
3090 throw "StandardModelMatching::mc_Q(): order " + out.str() + "not implemented";
3091 }
3092
3093 case QED0:
3094 break;
3095
3096 default:
3097 std::stringstream out;
3098 out << mcQ.getOrder_QED();
3099 throw "StandardModelMatching::mc_Q(): order " + out.str() + "not implemented";
3100
3101 }
3102
3103 return (mcQ);
3104}
3105
3106WilsonCoefficientNew& StandardModelMatching::mc_B() {
3107 double xt = x_t(Muw);
3108
3109 switch (mcB.getScheme()) {
3110 case NDR:
3111 //case HV:
3112 //case LRI:
3113 break;
3114 default:
3115 std::stringstream out;
3116 out << mcB.getScheme();
3117 throw "StandardModelMatching::mc_B(): scheme " + out.str() + "not implemented";
3118 }
3119
3120 mcB.setMu(Muw); // cleared too
3121
3122 switch (mcB.getOrder_QED()) {
3123 case QED2:
3124 case QED1:
3125 switch (mcB.getOrder_QCD()) {
3126 case QCD2:
3127 case QCD1:
3128 mcB.setCoeff(0, aletilde * (-1. / 2. / sW2 * S0(xt)), QCD1, QED1);
3129 case QCD0:
3130 break;
3131 default:
3132 std::stringstream out;
3133 out << mcB.getOrder_QCD();
3134 throw "StandardModelMatching::mc_B(): order " + out.str() + "not implemented";
3135 }
3136
3137 case QED0:
3138 break;
3139
3140 default:
3141 std::stringstream out;
3142 out << mcB.getOrder_QED();
3143 throw "StandardModelMatching::mc_B(): order " + out.str() + "not implemented";
3144
3145 }
3146
3147 return (mcB);
3148}
3149
3150unsigned int StandardModelMatching::setCMDF1(WilsonCoefficientNew& CMDF1, WilsonCoefficientNew& DF1block,
3151 unsigned int tot, schemes scheme, qcd_orders order_qcd, qed_orders order_qed) {
3152 unsigned int j, nops = DF1block.getSize();
3153
3154 for (qcd_orders ord_qcd = QCD0; ord_qcd <= order_qcd; ord_qcd = (qcd_orders) (ord_qcd + 1))
3155 for (qed_orders ord_qed = QED0; ord_qed <= order_qed; ord_qed = (qed_orders) (ord_qed + 1))
3156 for (j = 0; j < nops; j++)
3157 CMDF1.setCoeff(j + tot, DF1block.getCoeff(ord_qcd, ord_qed)(j), ord_qcd, ord_qed);
3158
3159 return (nops + tot);
3160}
3161
3162std::vector<WilsonCoefficientNew>& StandardModelMatching::CMDF1(std::string blocks, unsigned int nops) {
3163 unsigned int tot = 0;
3164 schemes scheme = mcC.getScheme();
3165 qcd_orders order_qcd = mcC.getOrder_QCD();
3166 qed_orders order_qed = mcC.getOrder_QED();
3167 vmcDF1.clear();
3168 WilsonCoefficientNew mcDF1(nops, scheme, order_qcd, order_qed);
3169
3170 typedef WilsonCoefficientNew & (StandardModelMatching::*BlockM)();
3171 std::vector< std::pair<std::string, BlockM> > Methods = {
3172 std::make_pair("C", &StandardModelMatching::mc_C),
3173 std::make_pair("P", &StandardModelMatching::mc_P),
3174 std::make_pair("M", &StandardModelMatching::mc_M),
3175 std::make_pair("L", &StandardModelMatching::mc_L),
3176 std::make_pair("Q", &StandardModelMatching::mc_Q),
3177 std::make_pair("B", &StandardModelMatching::mc_B)
3178 };
3179
3180 vmcDF1.clear();
3181 mcDF1.setMu(Muw);
3182 for (std::vector< std::pair<std::string, BlockM> >::iterator it = Methods.begin(); it != Methods.end(); it++)
3183 if (blocks.find(it->first) != std::string::npos)
3184 tot = setCMDF1(mcDF1, (this->*(it->second))(), tot, scheme, order_qcd, order_qed);
3185
3186 vmcDF1.push_back(mcDF1);
3187 return (vmcDF1);
3188}
FULLNNNLO
@ FULLNNNLO
Definition: OrderScheme.h:40
NNLO
@ NNLO
Definition: OrderScheme.h:36
LO
@ LO
Definition: OrderScheme.h:34
NLO
@ NLO
Definition: OrderScheme.h:35
FULLNNLO
@ FULLNNLO
Definition: OrderScheme.h:39
FULLNLO
@ FULLNLO
Definition: OrderScheme.h:38
HV
@ HV
Definition: OrderScheme.h:22
LRI
@ LRI
Definition: OrderScheme.h:23
NDR
@ NDR
Definition: OrderScheme.h:21
QED1
@ QED1
Definition: OrderScheme.h:92
QED0
@ QED0
Definition: OrderScheme.h:91
QED2
@ QED2
Definition: OrderScheme.h:93
NLO_QED11
@ NLO_QED11
Definition: OrderScheme.h:59
NLO_QED02
@ NLO_QED02
Definition: OrderScheme.h:61
NLO_QED12
@ NLO_QED12
Definition: OrderScheme.h:62
LO_QED
@ LO_QED
Definition: OrderScheme.h:58
NLO_QED22
@ NLO_QED22
Definition: OrderScheme.h:63
NLO_QED21
@ NLO_QED21
Definition: OrderScheme.h:60
QCD1
@ QCD1
Definition: OrderScheme.h:76
QCD0
@ QCD0
Definition: OrderScheme.h:75
QCD2
@ QCD2
Definition: OrderScheme.h:77
ModelMatching::CMBMll
virtual std::vector< WilsonCoefficient > & CMBMll(QCD::lepton lepton)=0
ModelMatching::CMBXsnn
virtual std::vector< WilsonCoefficient > & CMBXsnn(QCD::lepton lepton)=0
ModelMatching::CMDF1
virtual std::vector< WilsonCoefficientNew > & CMDF1(std::string blocks, unsigned int nops)=0
ModelMatching::CMd1
virtual std::vector< WilsonCoefficient > & CMd1()=0
ModelMatching::CMd1Buras
virtual std::vector< WilsonCoefficient > & CMd1Buras()=0
ModelMatching::CMprimebsg
virtual std::vector< WilsonCoefficient > & CMprimebsg()=0
ModelMatching::CMbnlep
virtual std::vector< WilsonCoefficient > & CMbnlep(int a)=0
ModelMatching::CMbnlepCC
virtual std::vector< WilsonCoefficient > & CMbnlepCC(const int a)=0
ModelMatching::CMprimeBMll
virtual std::vector< WilsonCoefficient > & CMprimeBMll(QCD::lepton lepton)=0
ModelMatching::CMbsg
virtual std::vector< WilsonCoefficient > & CMbsg()=0
Mw
An observable class for the -boson mass.
Definition: Mw.h:22
QCD::TOP
@ TOP
Definition: QCD.h:328
QCD::CHARM
@ CHARM
Definition: QCD.h:326
QCD::lepton
lepton
An enum type for leptons.
Definition: QCD.h:310
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:509
StandardModelMatching::updateSMParameters
void updateSMParameters()
Updates to new Standard Model parameter sets.
Definition: StandardModelMatching.cpp:110
StandardModelMatching::CMdbs2
virtual std::vector< WilsonCoefficient > & CMdbs2()
,
Definition: StandardModelMatching.cpp:1016
StandardModelMatching::CMdbd2
virtual std::vector< WilsonCoefficient > & CMdbd2()
,
Definition: StandardModelMatching.cpp:966
StandardModelMatching::StandardModelMatching
StandardModelMatching(const StandardModel &SM_i)
Definition: StandardModelMatching.cpp:16
StandardModelMatching::CMdd2
virtual std::vector< WilsonCoefficient > & CMdd2()
,
Definition: StandardModelMatching.cpp:1169
WilsonCoefficientNew
A class for the Wilson coefficients.
Definition: WilsonCoefficientNew.h:22
WilsonCoefficientNew::getCoeff
gslpp::vector< gslpp::complex > getCoeff(qcd_orders order_qcd_i, qed_orders order_qed_i=QED0) const
Definition: WilsonCoefficientNew.cpp:58
WilsonTemplateNew::getSize
unsigned int getSize() const
Definition: WilsonTemplateNew.h:86
lambda1
parameter of the Higgs potential
Definition: THDMquantities.h:382
orders
orders
An enum type for orders in QCD.
Definition: OrderScheme.h:33
schemes
schemes
An enum type for regularization schemes.
Definition: OrderScheme.h:20
qed_orders
qed_orders
An enum type for qed_orders in electroweak.
Definition: OrderScheme.h:90
orders_qed
orders_qed
An enum type for orders in electroweak.
Definition: OrderScheme.h:56
qcd_orders
qcd_orders
An enum type for qcd_orders in QCD.
Definition: OrderScheme.h:74
Li2
cd Li2(cd x)
Definition: hpl.h:1011