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