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