a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
BXqll.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2016 HEPfit Collaboration
3 *
4 *
5 * For the licensing terms see doc/COPYING.
6 */
7
8#include <gsl/gsl_sf.h>
9#include <gslpp.h>
10#include <boost/bind/bind.hpp>
11#include "std_make_vector.h"
12//#include <limits>
13#include "BXqll.h"
14#include "StandardModel.h"
15#include "gslpp_function_adapter.h"
16
17using namespace boost::placeholders;
18
19#define MPI2 (M_PI * M_PI)
20
21BXqll::BXqll(const StandardModel& SM_i, QCD::quark quark_i, QCD::lepton lep_i)
22: mySM(SM_i), myF_1(), myF_2(), myHeff("CPMLQB", SM_i, QCD2, QED2), WC(15, 9, 0.)
23{
24 lep = lep_i;
25 quark = quark_i;
26 w_H = gsl_integration_cquad_workspace_alloc(100);
27 CF = mySM.getCF();
28 phi1 = 50./3. - 8. * MPI2 / 3.; // functions for Phi_u at nh = 2 and nl = 3
29 phi2 = 2. * (- 2048. / 9. * gslpp_special_functions::zeta(3) + 16987. / 54. -340. / 81. * MPI2)
30 + 3. * (256. / 9. * gslpp_special_functions::zeta(3) -1009. / 27. + 308. / 81. * MPI2)
31 - 41848. / 81. * gslpp_special_functions::zeta(3) + 578. / 81. * MPI2 * MPI2
32 - 104480. / 729. * MPI2 + 1571095. / 1458. - 848. / 27. * MPI2 * log(2.);
33 QCD_max = 3;
34 QED_max = 3;
35}
36
37BXqll::~BXqll()
38{
39}
40
41std::vector<std::string> BXqll::initializeBXqllParameters()
42{
43 BXqllParameters = make_vector<std::string>() << "BR_BXcenu" << "C_ratio" << "BI_lambda1" << "BI_lambda2";
44
45 return (BXqllParameters);
46}
47
48void BXqll::updateParameters()
49{
50 GF = mySM.getGF();
51 Mlep = mySM.getLeptons(lep).getMass();
52 mu_b = mySM.getMub();
53 mu_c = mySM.getMuc();
54 Mb = mySM.getQuarks(QCD::BOTTOM).getMass(); // add the PS b mass
55 Mc = mySM.getQuarks(QCD::CHARM).getMass();
56 Ms = mySM.getQuarks(QCD::STRANGE).getMass();
57 MW = mySM.Mw();
58 abslambdat_over_Vcb = mySM.getCKM().computelamt_s().abs()/mySM.getCKM().getV_cb().abs();
59 alsmu = mySM.Als(mu_b, FULLNNNLO, true, true);
60 alsmuc = mySM.Als(mu_c, FULLNNNLO, true, true);
61 ale = mySM.Ale(mu_b, FULLNLO);
62 alstilde = alsmu / 4. / M_PI;
63 aletilde = ale / 4. / M_PI;
64 kappa = ale / alsmu;
65 Mtau = mySM.getLeptons(QCD::TAU).getMass(); // pole mass?
66 Mb_pole = mySM.Mbar2Mp(Mb, QCD::BOTTOM, FULLNNLO);
67 //Mc_pole = mySM.Mbar2Mp(Mc, FULLNNLO); //*** Mbar2Mp does not receive Mc ***/
68 Mc_pole = Mc*(1.+4.*alsmuc/3./M_PI+alsmuc*alsmuc/M_PI/M_PI*(-1.0414*(1.-4.*Ms/3.*Mc)+13.4434));
69 muh = mu_b/Mb_pole; // log(muh) uses the pole mass as stated in hep-ph/9910220
70 z = Mc_pole*Mc_pole/Mb_pole/Mb_pole; //****** Must be pole masses ****/
71 Lbl = 2. * log(Mb_pole / Mlep); // Mb,pole ?
72
73 BR_BXcenu = mySM.getOptionalParameter("BR_BXcenu"); // Branching ratio of B -> Xc e nu
74 C_ratio = mySM.getOptionalParameter("C_ratio"); // Ratio of branching ratios as defined by Gambino, Misiak, arXiv:hep-ph/0104034
75
76 pre = BR_BXcenu * abslambdat_over_Vcb * abslambdat_over_Vcb * 4. / C_ratio;
77
78 lambda_1 = mySM.getOptionalParameter("BI_lambda1");
79 lambda_2 = mySM.getOptionalParameter("BI_lambda2");
80
81 // Auxiliary variables for Phi_u_inv
82 phi00 = 1. + (lambda_1 - 9. * lambda_2) / 2. / Mb_pole / Mb_pole;
83 phi00_2 = phi00 * phi00;
84 phi20 = phi2 + 2. * mySM.Beta0(5) * phi1 * log(muh);
85 phi01 = 12. / 23. * (1. - alsmu / mySM.Als(mySM.getMuw(), FULLNNNLO, true, true));
86
87 phinv00 = 1. / phi00;
88 phinv10 = alstilde * (- phi1 / phi00_2);
89 phinv20 = alstilde * alstilde * (phi1 * phi1 / phi00_2 / phi00 - phi20 / phi00_2);
90 phinv01 = kappa * (- phi01 / phi00_2);
91 phinv11 = alstilde * kappa * 2. * phi1 * phi01 / phi00_2 / phi00;
92 phinv21 = alstilde * alstilde * kappa * (2. * phi20 * phi01 / phi00_2 / phi00 -
93 3. * phi1 * phi1 * phi01 / phi00_2 / phi00_2);
94
95 phinv_00_01 = phinv00 + phinv01;
96
97 //ISIDORI VALUES
98// z = 0.29*0.29;
99// mu_b = 5.0;
100// Mb = 4.9;
101// Mtau = 1.77;
102// muh = mu_b/Mb;
103// ale = 0.0078125;
104// abslambdat_over_Vcb = 0.97;
105// Vts_over_Vcb = 0.97;
106// alsmu = 0.215;
107
108 //RYOUTARO'S VALUES
109// mu_b = 5.0;
110// alstilde = 0.0170686027;
111// ale = 0.00757373;
112// aletilde = ale / 4. / M_PI;
113// kappa = aletilde / alstilde;
114// Mb_pole = 4.8;
115// Mc_pole = 2.04545;
116
117// allcoeff_smm = mySM.getFlavour().ComputeCoeffsmumu(mu_b, NDR);
118 allcoeff = myHeff.ComputeCoeff(mu_b);
119
120 double alpha_kappa;
121
122 for(unsigned int qed_ord = QED0; qed_ord <= QED2; qed_ord++)
123 for(unsigned int qcd_ord = QCD0; qcd_ord <= QCD2; qcd_ord++)
124 {
125 alpha_kappa = pow(alstilde, qcd_ord) * pow(kappa, qed_ord);
126 for (unsigned int i = 0; i < 15; i++)
127 if (i != 8 && i != 9 && qed_ord < 2 && qcd_ord * qed_ord < 2)
128 WC.assign(i, qcd_ord + 3*qed_ord, alpha_kappa * (myHeff.LowScaleCoeff((qcd_orders) qcd_ord, (qed_orders) qed_ord))(i));
129 WC.assign(8, qcd_ord + 3*qed_ord, alpha_kappa * (myHeff.LowScaleCoeff((qcd_orders) qcd_ord, (qed_orders) qed_ord))(8));
130 WC.assign(9, qcd_ord + 3*qed_ord, alpha_kappa * (myHeff.LowScaleCoeff((qcd_orders) qcd_ord, (qed_orders) qed_ord))(9));
131 }
132}
133
134double BXqll::getR_LOWQ2(double sh)
135{
136 updateParameters();
137
138 //To test HeffDF1 Wilson coefficients and Expanded multiplications
139 Test_WC_DF1();
140 return 0.;
141
142// computeMi(sh);
143// return H_A(sh);
144}
145
146gslpp::complex BXqll::F19(double sh)
147{
148 gslpp::complex i = gslpp::complex::i();
149
150 return (F_19re(muh, z, sh) + i*F_19im(muh, z, sh));
151}
152
153gslpp::complex BXqll::F29(double sh)
154{
155 gslpp::complex i = gslpp::complex::i();
156
157 return (F_29re(muh, z, sh) + i*F_29im(muh, z, sh));
158}
159
160gslpp::complex BXqll::F17(double sh)
161{
162 gslpp::complex i = gslpp::complex::i();
163
164 return (F_17re(muh, z, sh) + i*F_17im(muh, z, sh));
165}
166
167gslpp::complex BXqll::F27(double sh)
168{
169 gslpp::complex i = gslpp::complex::i();
170
171 return (F_27re(muh, z, sh) + i*F_27im(muh, z, sh));
172}
173
174gslpp::complex BXqll::F87(double sh)
175{
176 gslpp::complex i = gslpp::complex::i();
177 double ash = asin(sqrt(sh)/2.);
178 double umsh = 1.-sh;
179
180 return (4.*M_PI*M_PI/27.*(2.+sh)/umsh/umsh/umsh/umsh-4./9.*(11.-16.*sh+8.*sh*sh)/umsh/umsh-
181 8./9.*sqrt(sh*(4.-sh))/umsh/umsh/umsh*(9.-5.*sh+2.*sh*sh)*ash-16./3.*(2+sh)/
182 umsh/umsh/umsh/umsh*ash*ash-8./9.*sh/umsh*log(sh)-32./9.*log(muh)-i*8./9.*M_PI);
183}
184
185double BXqll::F89(double sh)
186{
187 double ash = asin(sqrt(sh)/2.);
188 double umsh = 1.-sh;
189
190 return (-8.*M_PI*M_PI/27.*(4.-sh)/umsh/umsh/umsh/umsh+8./9.*(5.-2.*sh)/umsh/umsh+
191 16./9.*sqrt((4.-sh)/sh)/umsh/umsh/umsh*(4.+3.*sh-sh*sh)*ash+32./3.*(4.-sh)/
192 umsh/umsh/umsh/umsh*ash*ash+16./9./umsh*log(sh));
193}
194
195double BXqll::F_17re(double muh, double z, double sh, int maxpow)
196{
197 return myF_1.F_17re(muh, z, sh, maxpow);
198};
199
200double BXqll::F_17im(double muh, double z, double sh, int maxpow)
201{
202 return myF_1.F_17im(muh, z, sh, maxpow);
203};
204
205double BXqll::F_19re(double muh, double z, double sh, int maxpow)
206{
207 return myF_1.F_19re(muh, z, sh, maxpow);
208};
209
210double BXqll::F_19im(double muh, double z, double sh, int maxpow)
211{
212 return myF_1.F_19im(muh, z, sh, maxpow);
213};
214
215double BXqll::F_27re(double muh, double z, double sh, int maxpow)
216{
217 return myF_2.F_27re(muh, z, sh, maxpow);
218};
219
220double BXqll::F_27im(double muh, double z, double sh, int maxpow)
221{
222 return myF_2.F_27im(muh, z, sh, maxpow);
223};
224
225double BXqll::F_29re(double muh, double z, double sh, int maxpow)
226{
227 return myF_2.F_29re(muh, z, sh, maxpow);
228};
229
230double BXqll::F_29im(double muh, double z, double sh, int maxpow)
231{
232 return myF_2.F_29im(muh, z, sh, maxpow);
233};
234
235double BXqll::DeltaF_19re(double muh, double z, double sh, int maxpow)
236{
237 return myF_1.DeltaF_19re(muh, z, sh, maxpow);
238};
239
240double BXqll::DeltaF_19im(double muh, double z, double sh, int maxpow)
241{
242 return myF_1.DeltaF_19im(muh, z, sh, maxpow);
243};
244
245double BXqll::DeltaF_29re(double muh, double z, double sh, int maxpow)
246{
247 return myF_2.DeltaF_29re(muh, z, sh, maxpow);
248};
249
250double BXqll::DeltaF_29im(double muh, double z, double sh, int maxpow)
251{
252 return myF_2.DeltaF_29im(muh, z, sh, maxpow);
253};
254
255
256/*********************************************************
257 * Implementation of the notation of @cite Huber:2015sra *
258 *********************************************************/
259
260double BXqll::integrateH(std::string obs, double q_min, double q_max)
261{
262 updateParameters();
263
264 old_handler = gsl_set_error_handler_off();
265
266 double sh_min = q_min/Mb_pole/Mb_pole, sh_max = q_max/Mb_pole/Mb_pole; // pole mass as explicitly stated in hep-ph/051206
267
268 FH = convertToGslFunction(bind(&BXqll::getH, &(*this), obs, _1));
269
270 if (gsl_integration_cquad(&FH, sh_min, sh_max, 1.e-5, 1.e-4, w_H, &aveH, &errH, NULL) != 0)
271 return std::numeric_limits<double>::quiet_NaN();
272 return aveH;
273
274 gsl_set_error_handler(old_handler);
275}
276
277double BXqll::getH(std::string obs, double sh)
278{
279 updateParameters();
280 computeMi(sh);
281
282 if (obs == "T")
283 return H_T(sh);
284 else if (obs == "L")
285 return H_L(sh);
286 else if (obs == "A")
287// return H_A(sh);
288 return (H_A(sh) + pre * PhiA_brems(sh) * phinv_00_01);
289 else if (obs == "TL")
290// return (H_T(sh) + H_L(sh));
291 return (H_T(sh) + H_L(sh) + pre * PhiTL_brems(sh) * phinv_00_01);
292 else
293 throw std::runtime_error("BXqll::getH: Angular observable not implemented");
294}
295
296double BXqll::H_T(double sh)
297{
298 // Clears the vector for a new value of sh
299 Hij_T.clear();
300
301 gslpp::matrix<gslpp::complex> Hij(15, 15, 0.);
302
303 // LO
304 Hij.assign(8, 8, M_9[QCD0 + 3*QED0](8).abs2() * S99_T(sh, LO) );
305 Hij.assign(9, 9, M_10[QCD0 + 3*QED0](9).abs2() * S1010_T(sh, LO) );
306
307 Hij_T.push_back(Hij);
308
309
310 // NLO
311 Hij.reset(); // To clear Hij
312
313 Hij.assign(8, 8, M_9[QCD0 + 3*QED0](8).abs2() * S99_T(sh, NLO) );
314 Hij.assign(9, 9, M_10[QCD0 + 3*QED0](9).abs2() * S1010_T(sh, NLO) );
315
316 Hij_T.push_back(Hij);
317
318
319 // NNLO
320 Hij.reset(); // To clear Hij
321
322 Hij.assign(8, 8, M_9[QCD0 + 3*QED0](8).abs2() * S99_T(sh, NNLO) );
323 Hij.assign(9, 9, M_10[QCD0 + 3*QED0](9).abs2() * S1010_T(sh, NNLO) );
324
325 Hij_T.push_back(Hij);
326
327
328 // LO_QED
329 Hij.reset(); // To clear Hij
330 Hij_T.push_back(Hij);
331
332
333 // NLO_QED11
334 for (unsigned int j = 0; j < 15; j++)
335 {
336 for (unsigned int i = 0; i <= j; i++)
337 {
338 if (i == j)
339 Hij.assign(i, j,
340 M_9[QCD0 + 3*QED0](i) * (M_9[QCD1 + 3*QED1](i).conjugate()).real() * 2. * S99_T(sh, LO) );
341
342 else
343 Hij.assign(i, j,
344 2. * (M_9[QCD0 + 3*QED0](i) * M_9[QCD1 + 3*QED1](j).conjugate() +
345 M_9[QCD1 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate()) * S99_T(sh, LO) +
346 (M_7[QCD1 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate() +
347 M_9[QCD0 + 3*QED0](i) * M_7[QCD1 + 3*QED1](j).conjugate()) * S79_T(sh, LO) );
348 }
349 }
350
351 Hij_T.push_back(Hij);
352
353
354 // NLO_QED21
355 Hij.reset(); // To clear Hij
356
357 for (unsigned int j = 0; j < 15; j++)
358 {
359 for (unsigned int i = 0; i <= j; i++)
360 {
361 if (i == j)
362 Hij.assign(i, j,
363 (M_9[QCD0 + 3*QED0](i) * M_9[QCD1 + 3*QED1](i).conjugate()).real() * 2. * S99_T(sh, NLO) );
364
365 else
366 Hij.assign(i, j,
367 2. * (M_9[QCD0 + 3*QED0](i) * M_9[QCD2 + 3*QED1](j).conjugate() +
368 M_9[QCD2 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate()) * S99_T(sh, LO) +
369 2. * (M_9[QCD0 + 3*QED0](i) * M_9[QCD1 + 3*QED1](j).conjugate() +
370 M_9[QCD1 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate()) * S99_T(sh, NLO) +
371 (M_7[QCD2 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate() +
372 M_9[QCD0 + 3*QED0](i) * M_7[QCD2 + 3*QED1](j).conjugate()) * S79_T(sh, LO) +
373 (M_7[QCD1 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate() +
374 M_9[QCD0 + 3*QED0](i) * M_7[QCD1 + 3*QED1](j).conjugate()) * S79_T(sh, NLO) );
375 }
376 }
377
378 Hij_T.push_back(Hij);
379
380
381 // NLO_QED02 -> NLO_QED12
382 Hij.reset(); // To clear Hij
383 Hij_T.push_back(Hij);
384 Hij_T.push_back(Hij);
385
386
387 // NLO_QED22
388 for (unsigned int j = 0; j < 15; j++)
389 {
390 for (unsigned int i = 0; i <= j; i++)
391 {
392 if (i == j)
393 Hij.assign(i, j,
394 M_7[QCD1 + 3*QED1](i).abs2() * S77_T(sh, LO) +
395 M_9[QCD1 + 3*QED1](i).abs2() * S99_T(sh, LO) );
396
397 else
398 Hij.assign(i, j,
399 2. * M_7[QCD1 + 3*QED1](i) * M_7[QCD1 + 3*QED1](j).conjugate() * S77_T(sh, LO) +
400 2. * M_9[QCD1 + 3*QED1](i) * M_9[QCD1 + 3*QED1](j).conjugate() * S99_T(sh, LO) +
401 (M_7[QCD1 + 3*QED1](i) * M_9[QCD1 + 3*QED1](j).conjugate() +
402 M_9[QCD1 + 3*QED1](i) * M_7[QCD1 + 3*QED1](j).conjugate()) * S79_T(sh, LO) );
403 }
404 }
405
406 Hij_T.push_back(Hij);
407
408
409 // Additional entry for 31
410 Hij.reset(); // To clear Hij
411
412 for (unsigned int j = 0; j < 15; j++)
413 {
414 for (unsigned int i = 0; i <= j; i++)
415 {
416 if (i == j)
417 Hij.assign(i, j,
418 (M_9[QCD0 + 3*QED0](i) * M_9[QCD1 + 3*QED1](i).conjugate()).real() * 2. * S99_T(sh, NNLO) );
419
420 else
421 Hij.assign(i, j,
422 2. * (M_9[QCD0 + 3*QED0](i) * M_9[QCD2 + 3*QED1](j).conjugate() +
423 M_9[QCD2 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate()) * S99_T(sh, NLO) +
424 2. * (M_9[QCD0 + 3*QED0](i) * M_9[QCD1 + 3*QED1](j).conjugate() +
425 M_9[QCD1 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate()) * S99_T(sh, NNLO) +
426 (M_7[QCD2 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate() +
427 M_9[QCD0 + 3*QED0](i) * M_7[QCD2 + 3*QED1](j).conjugate()) * S79_T(sh, NLO) );
428 }
429 }
430
431 Hij_T.push_back(Hij);
432
433
434 // Additional entry for 32
435 Hij.reset(); // To clear Hij
436
437 for (unsigned int j = 0; j < 15; j++)
438 {
439 for (unsigned int i = 0; i <= j; i++)
440 {
441 if (i == j)
442 Hij.assign(i, j,
443 M_7[QCD1 + 3*QED1](i).abs2() * S77_T(sh, NLO) +
444 M_9[QCD1 + 3*QED1](i).abs2() * S99_T(sh, NLO) +
445 (M_9[QCD1 + 3*QED1](i) * M_9[QCD2 + 3*QED1](i).conjugate()).real() * 2. * S99_T(sh, LO) +
446 (M_7[QCD2 + 3*QED1](i) * M_9[QCD1 + 3*QED1](i)).real() * S79_T(sh, LO));
447
448 else
449 Hij.assign(i, j,
450 2. * M_7[QCD1 + 3*QED1](i) * M_7[QCD1 + 3*QED1](j).conjugate() * S77_T(sh, NLO) +
451 2. * (M_7[QCD1 + 3*QED1](i) * M_7[QCD2 + 3*QED1](j).conjugate() +
452 M_7[QCD2 + 3*QED1](i) * M_7[QCD1 + 3*QED1](j).conjugate()) * S77_T(sh, LO) +
453 2. * M_9[QCD1 + 3*QED1](i) * M_9[QCD1 + 3*QED1](j).conjugate() * S99_T(sh, NLO) +
454 2. * (M_9[QCD1 + 3*QED1](i) * M_9[QCD2 + 3*QED1](j).conjugate() +
455 M_9[QCD2 + 3*QED1](i) * M_9[QCD1 + 3*QED1](j).conjugate()) * S99_T(sh, LO) +
456 (M_7[QCD2 + 3*QED1](i) * M_9[QCD1 + 3*QED1](j).conjugate() +
457 M_7[QCD1 + 3*QED1](i) * M_9[QCD2 + 3*QED1](j).conjugate() +
458 M_9[QCD2 + 3*QED1](i) * M_7[QCD1 + 3*QED1](j).conjugate() +
459 M_9[QCD1 + 3*QED1](i) * M_7[QCD2 + 3*QED1](j).conjugate()) * S79_T(sh, LO) +
460 (M_7[QCD1 + 3*QED1](i) * M_9[QCD1 + 3*QED1](j).conjugate() +
461 M_9[QCD1 + 3*QED1](i) * M_7[QCD1 + 3*QED1](j).conjugate()) * S79_T(sh, NLO) );
462 }
463 }
464
465 Hij_T.push_back(Hij);
466
467 // Additional entry for 33
468 Hij.reset(); // To clear Hij
469 Hij_T.push_back(Hij);
470
471 // 1/mc^2 corrections
472 for(unsigned int i = 0; i < 2; i++)
473 for(unsigned int j = 0; j < 15; j++)
474 {
475 Hij_T[QCD1 + 3*QED1].assign(i, j, Hij_T[QCD1 + 3*QED1](i, j) + cij_T(i, j, sh, 11));
476 Hij_T[QCD2 + 3*QED2].assign(i, j, Hij_T[QCD2 + 3*QED2](i, j) + cij_T(i, j, sh, 22));
477 Hij_T[QCD2 + 3*QED2 + 2].assign(i, j, Hij_T[QCD2 + 3*QED2 + 2](i, j) + cij_T(i, j, sh, 32));
478 }
479
480 // log-enhanced electromagnetic corrections
481 for(unsigned int j = 0; j < 7; j++)
482 for(unsigned int i = 0; i <= j; i++)
483 Hij_T[QCD2 + 3*QED2 + 3].assign(i, j, Hij_T[QCD2 + 3*QED2 + 3](i, j) + eij_T(i, j, sh));
484
485 for(unsigned int i = 0; i < 7; i++)
486 Hij_T[QCD2 + 3*QED2].assign(i, 8, Hij_T[QCD2 + 3*QED2](i, 8) + eij_T(i, 8, sh));
487
488 Hij_T[QCD1 + 3*QED1].assign(8, 8, Hij_T[QCD1 + 3*QED1](8, 8) + eij_T(8, 8, sh));
489 Hij_T[QCD1 + 3*QED1].assign(9, 9, Hij_T[QCD1 + 3*QED1](9, 9) + eij_T(9, 9, sh));
490
491
492 return pre * CCH_multiplication(Hij_T);
493}
494
495double BXqll::H_L(double sh)
496{
497 // Clears the vector for a new value of sh
498 Hij_L.clear();
499
500 gslpp::matrix<gslpp::complex> Hij(15, 15, 0.);
501
502 // LO
503 Hij.assign(8, 8, M_9[QCD0 + 3*QED0](8).abs2() * S99_L(sh, LO) );
504 Hij.assign(9, 9, M_10[QCD0 + 3*QED0](9).abs2() * S1010_L(sh, LO) );
505
506 Hij_L.push_back(Hij);
507
508
509 // NLO
510 Hij.reset(); // To clear Hij
511
512 Hij.assign(8, 8, M_9[QCD0 + 3*QED0](8).abs2() * S99_L(sh, NLO) );
513 Hij.assign(9, 9, M_10[QCD0 + 3*QED0](9).abs2() * S1010_L(sh, NLO) );
514
515 Hij_L.push_back(Hij);
516
517
518 // NNLO
519 Hij.reset(); // To clear Hij
520
521 Hij.assign(8, 8, M_9[QCD0 + 3*QED0](8).abs2() * S99_L(sh, NNLO) );
522 Hij.assign(9, 9, M_10[QCD0 + 3*QED0](9).abs2() * S1010_L(sh, NNLO) );
523
524 Hij_L.push_back(Hij);
525
526
527 // LO_QED
528 Hij.reset(); // To clear Hij
529 Hij_L.push_back(Hij);
530
531
532 // NLO_QED11
533 for (unsigned int j = 0; j < 15; j++)
534 {
535 for (unsigned int i = 0; i <= j; i++)
536 {
537 if (i == j)
538 Hij.assign(i, j,
539 M_9[QCD0 + 3*QED0](i) * (M_9[QCD1 + 3*QED1](i).conjugate()).real() * 2. * S99_L(sh, LO) );
540
541 else
542 Hij.assign(i, j,
543 2. * (M_9[QCD0 + 3*QED0](i) * M_9[QCD1 + 3*QED1](j).conjugate() +
544 M_9[QCD1 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate()) * S99_L(sh, LO) +
545 (M_7[QCD1 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate() +
546 M_9[QCD0 + 3*QED0](i) * M_7[QCD1 + 3*QED1](j).conjugate()) * S79_L(sh, LO) );
547 }
548 }
549
550 Hij_L.push_back(Hij);
551
552
553 // NLO_QED21
554 Hij.reset(); // To clear Hij
555
556 for (unsigned int j = 0; j < 15; j++)
557 {
558 for (unsigned int i = 0; i <= j; i++)
559 {
560 if (i == j)
561 Hij.assign(i, j,
562 (M_9[QCD0 + 3*QED0](i) * M_9[QCD1 + 3*QED1](i).conjugate()).real() * 2. * S99_L(sh, NLO) );
563
564 else
565 Hij.assign(i, j,
566 2. * (M_9[QCD0 + 3*QED0](i) * M_9[QCD2 + 3*QED1](j).conjugate() +
567 M_9[QCD2 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate()) * S99_L(sh, LO) +
568 2. * (M_9[QCD0 + 3*QED0](i) * M_9[QCD1 + 3*QED1](j).conjugate() +
569 M_9[QCD1 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate()) * S99_L(sh, NLO) +
570 (M_7[QCD2 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate() +
571 M_9[QCD0 + 3*QED0](i) * M_7[QCD2 + 3*QED1](j).conjugate()) * S79_L(sh, LO) +
572 (M_7[QCD1 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate() +
573 M_9[QCD0 + 3*QED0](i) * M_7[QCD1 + 3*QED1](j).conjugate()) * S79_L(sh, NLO) );
574 }
575 }
576
577 Hij_L.push_back(Hij);
578
579
580 // NLO_QED02 -> NLO_QED12
581 Hij.reset(); // To clear Hij
582 Hij_L.push_back(Hij);
583 Hij_L.push_back(Hij);
584
585
586 // NLO_QED22
587 for (unsigned int j = 0; j < 15; j++)
588 {
589 for (unsigned int i = 0; i <= j; i++)
590 {
591 if (i == j)
592 Hij.assign(i, j,
593 M_7[QCD1 + 3*QED1](i).abs2() * S77_L(sh, LO) +
594 M_9[QCD1 + 3*QED1](i).abs2() * S99_L(sh, LO) );
595
596 else
597 Hij.assign(i, j,
598 2. * M_7[QCD1 + 3*QED1](i) * M_7[QCD1 + 3*QED1](j).conjugate() * S77_L(sh, LO) +
599 2. * M_9[QCD1 + 3*QED1](i) * M_9[QCD1 + 3*QED1](j).conjugate() * S99_L(sh, LO) +
600 (M_7[QCD1 + 3*QED1](i) * M_9[QCD1 + 3*QED1](j).conjugate() +
601 M_9[QCD1 + 3*QED1](i) * M_7[QCD1 + 3*QED1](j).conjugate()) * S79_L(sh, LO) );
602 }
603 }
604
605 Hij_L.push_back(Hij);
606
607
608 // Additional entry for 31
609 Hij.reset(); // To clear Hij
610
611 for (unsigned int j = 0; j < 15; j++)
612 {
613 for (unsigned int i = 0; i <= j; i++)
614 {
615 if (i == j)
616 Hij.assign(i, j,
617 (M_9[QCD0 + 3*QED0](i) * M_9[QCD1 + 3*QED1](i).conjugate()).real() * 2. * S99_L(sh, NNLO) );
618
619 else
620 Hij.assign(i, j,
621 2. * (M_9[QCD0 + 3*QED0](i) * M_9[QCD2 + 3*QED1](j).conjugate() +
622 M_9[QCD2 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate()) * S99_L(sh, NLO) +
623 2. * (M_9[QCD0 + 3*QED0](i) * M_9[QCD1 + 3*QED1](j).conjugate() +
624 M_9[QCD1 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate()) * S99_L(sh, NNLO) +
625 (M_7[QCD2 + 3*QED1](i) * M_9[QCD0 + 3*QED0](j).conjugate() +
626 M_9[QCD0 + 3*QED0](i) * M_7[QCD2 + 3*QED1](j).conjugate()) * S79_L(sh, NLO) );
627 }
628 }
629
630 Hij_L.push_back(Hij);
631
632
633 // Additional entry for 32
634 Hij.reset(); // To clear Hij
635
636 for (unsigned int j = 0; j < 15; j++)
637 {
638 for (unsigned int i = 0; i <= j; i++)
639 {
640 if (i == j)
641 Hij.assign(i, j,
642 M_7[QCD1 + 3*QED1](i).abs2() * S77_L(sh, NLO) +
643 M_9[QCD1 + 3*QED1](i).abs2() * S99_L(sh, NLO) +
644 (M_9[QCD1 + 3*QED1](i) * M_9[QCD2 + 3*QED1](i).conjugate()).real() * 2. * S99_L(sh, LO) +
645 (M_7[QCD2 + 3*QED1](i) * M_9[QCD1 + 3*QED1](i)).real() * S79_L(sh, LO));
646
647 else
648 Hij.assign(i, j,
649 2. * M_7[QCD1 + 3*QED1](i) * M_7[QCD1 + 3*QED1](j).conjugate() * S77_L(sh, NLO) +
650 2. * (M_7[QCD1 + 3*QED1](i) * M_7[QCD2 + 3*QED1](j).conjugate() +
651 M_7[QCD2 + 3*QED1](i) * M_7[QCD1 + 3*QED1](j).conjugate()) * S77_L(sh, LO) +
652 2. * M_9[QCD1 + 3*QED1](i) * M_9[QCD1 + 3*QED1](j).conjugate() * S99_L(sh, NLO) +
653 2. * (M_9[QCD1 + 3*QED1](i) * M_9[QCD2 + 3*QED1](j).conjugate() +
654 M_9[QCD2 + 3*QED1](i) * M_9[QCD1 + 3*QED1](j).conjugate()) * S99_L(sh, LO) +
655 (M_7[QCD1 + 3*QED1](i) * M_9[QCD1 + 3*QED1](j).conjugate() +
656 M_9[QCD1 + 3*QED1](i) * M_7[QCD1 + 3*QED1](j).conjugate()) * S79_L(sh, NLO) +
657 (M_7[QCD2 + 3*QED1](i) * M_9[QCD1 + 3*QED1](j).conjugate() +
658 M_7[QCD1 + 3*QED1](i) * M_9[QCD2 + 3*QED1](j).conjugate() +
659 M_9[QCD2 + 3*QED1](i) * M_7[QCD1 + 3*QED1](j).conjugate() +
660 M_9[QCD1 + 3*QED1](i) * M_7[QCD2 + 3*QED1](j).conjugate()) * S79_L(sh, LO) );
661 }
662 }
663
664 Hij_L.push_back(Hij);
665
666 // Additional entry for 33
667 Hij.reset(); // To clear Hij
668 Hij_L.push_back(Hij);
669
670 // 1/mc^2 corrections
671 for(unsigned int i = 0; i < 2; i++)
672 for(unsigned int j = 0; j < 15; j++)
673 {
674 Hij_L[QCD1 + 3*QED1].assign(i, j, Hij_L[QCD1 + 3*QED1](i, j) + cij_L(i, j, sh, 11));
675 Hij_L[QCD2 + 3*QED2].assign(i, j, Hij_L[QCD2 + 3*QED2](i, j) + cij_L(i, j, sh, 22));
676 Hij_L[QCD2 + 3*QED2 + 2].assign(i, j, Hij_L[QCD2 + 3*QED2 + 2](i, j) + cij_L(i, j, sh, 32));
677 }
678
679 // log-enhanced electromagnetic corrections
680 for(unsigned int j = 0; j < 7; j++)
681 for(unsigned int i = 0; i <= j; i++)
682 Hij_L[QCD2 + 3*QED2 + 3].assign(i, j, Hij_L[QCD2 + 3*QED2 + 3](i, j) + eij_L(i, j, sh));
683
684 for(unsigned int i = 0; i < 7; i++)
685 Hij_L[QCD2 + 3*QED2].assign(i, 8, Hij_L[QCD2 + 3*QED2](i, 8) + eij_L(i, 8, sh));
686
687 Hij_L[QCD1 + 3*QED1].assign(8, 8, Hij_L[QCD1 + 3*QED1](8, 8) + eij_L(8, 8, sh));
688 Hij_L[QCD1 + 3*QED1].assign(9, 9, Hij_L[QCD1 + 3*QED1](9, 9) + eij_L(9, 9, sh));
689
690
691 return pre * CCH_multiplication(Hij_L);
692}
693
694double BXqll::H_A(double sh)
695{
696 // Clears the vector for a new value of sh
697 Hij_A.clear();
698
699 gslpp::matrix<gslpp::complex> Hij(15, 15, 0.);
700
701 //REMINDER: M_10(i) = 1, for i = 9, zero otherwise
702 // LO
703 Hij.assign(8, 9, M_9[QCD0 + 3*QED0](8) * S910_A(sh, LO) );
704
705 Hij_A.push_back(Hij);
706
707
708 // NLO
709 Hij.reset(); // To clear Hij
710
711 Hij.assign(8, 9, M_9[QCD0 + 3*QED0](8) * S910_A(sh, NLO) );
712
713 Hij_A.push_back(Hij);
714
715
716 // NNLO
717 Hij.reset(); // To clear Hij
718
719 Hij.assign(8, 9, M_9[QCD0 + 3*QED0](8) * S910_A(sh, NNLO) );
720
721 Hij_A.push_back(Hij);
722
723
724 // LO_QED
725 Hij.reset(); // To clear Hij
726 Hij_A.push_back(Hij);
727
728
729 // NLO_QED11
730 Hij.assign(6, 9, M_7[QCD1 + 3*QED1](6) * S710_A(sh, LO) );
731
732 for (unsigned int i = 0; i < 9; i++)
733 if (i != 6)
734 Hij.assign(i, 9, M_9[QCD1 + 3*QED1](i) * S910_A(sh, LO) );
735
736 for (unsigned int j = 10; j < 15; j++)
737 Hij.assign(9, j, M_9[QCD1 + 3*QED1](j).conjugate() * S910_A(sh, LO) );
738
739 Hij_A.push_back(Hij);
740
741
742 // NLO_QED21
743 Hij.reset(); // To clear Hij
744
745 Hij.assign(0, 9, M_7[QCD2 + 3*QED1](0) * S710_A(sh, LO) +
746 M_9[QCD2 + 3*QED1](0) * S910_A(sh, LO) + M_9[QCD1 + 3*QED1](0) * S910_A(sh, NLO) );
747 Hij.assign(1, 9, M_7[QCD2 + 3*QED1](1) * S710_A(sh, LO) +
748 M_9[QCD2 + 3*QED1](1) * S910_A(sh, LO) + M_9[QCD1 + 3*QED1](1) * S910_A(sh, NLO) );
749 Hij.assign(6, 9, M_7[QCD1 + 3*QED1](6) * S710_A(sh, NLO) );
750 Hij.assign(7, 9, M_7[QCD2 + 3*QED1](7) * S710_A(sh, LO) + M_9[QCD2 + 3*QED1](7) * S910_A(sh, LO) );
751
752 for (unsigned int i = 2; i < 9; i++)
753 if (i != 6 && i != 7)
754 Hij.assign(i, 9, M_9[QCD1 + 3*QED1](i) * S910_A(sh, NLO) );
755
756 for (unsigned int j = 10; j < 15; j++)
757 Hij.assign(9, j, M_9[QCD1 + 3*QED1](j).conjugate() * S910_A(sh, NLO) );
758
759 Hij_A.push_back(Hij);
760
761
762 // NLO_QED02 -> NLO_QED22
763 Hij.reset(); // To clear Hij
764 Hij_A.push_back(Hij);
765 Hij_A.push_back(Hij);
766 Hij_A.push_back(Hij);
767
768 // Additional entry for 31
769 Hij.assign(0, 9, M_7[QCD2 + 3*QED1](0) * S710_A(sh, NLO) +
770 M_9[QCD2 + 3*QED1](0) * S910_A(sh, NLO) + M_9[QCD1 + 3*QED1](0) * S910_A(sh, NNLO) );
771 Hij.assign(1, 9, M_7[QCD2 + 3*QED1](1) * S710_A(sh, NLO) +
772 M_9[QCD2 + 3*QED1](1) * S910_A(sh, NLO) + M_9[QCD1 + 3*QED1](1) * S910_A(sh, NNLO) );
773 Hij.assign(7, 9, M_7[QCD2 + 3*QED1](7) * S710_A(sh, NLO) + M_9[QCD2 + 3*QED1](7) * S910_A(sh, NLO) );
774
775 for (unsigned int i = 2; i < 9; i++)
776 if (i != 7)
777 Hij.assign(i, 9, M_9[QCD1 + 3*QED1](i) * S910_A(sh, NNLO) );
778
779 for (unsigned int j = 10; j < 15; j++)
780 Hij.assign(9, j, M_9[QCD1 + 3*QED1](j).conjugate() * S910_A(sh, NNLO) );
781
782 Hij_A.push_back(Hij);
783
784 // Additional entries for 32, 33
785 Hij.reset(); // To clear Hij
786 Hij_A.push_back(Hij);
787 Hij_A.push_back(Hij);
788
789 // 1/mc^2 corrections associated with C10
790 Hij_A[QCD1 + 3*QED1].assign(0, 9, Hij_A[QCD1 + 3*QED1](0, 9) + cij_A(0, 9, sh));
791 Hij_A[QCD1 + 3*QED1].assign(1, 9, Hij_A[QCD1 + 3*QED1](1, 9) + cij_A(1, 9, sh));
792
793 // log-enhanced electromagnetic corrections
794 Hij_A[QCD2 + 3*QED2].assign(0, 9, Hij_A[QCD2 + 3*QED2](0, 9) + eij_A(0, 9, sh));
795 Hij_A[QCD2 + 3*QED2].assign(1, 9, Hij_A[QCD2 + 3*QED2](1, 9) + eij_A(1, 9, sh));
796 Hij_A[QCD2 + 3*QED2].assign(6, 9, Hij_A[QCD2 + 3*QED2](6, 9) + eij_A(6, 9, sh));
797 Hij_A[QCD1 + 3*QED1].assign(8, 9, Hij_A[QCD1 + 3*QED1](8, 9) + eij_A(8, 9, sh));
798
799
800 return pre * CCH_multiplication(Hij_A);
801}
802
803void BXqll::computeMi(double sh)
804{
805 // Clears the vectors for a new value of sh
806 M_7.clear();
807 M_9.clear();
808 M_10.clear();
809
810 gslpp::vector<gslpp::complex> M7i(15, 0.);
811 gslpp::vector<gslpp::complex> M9i(15, 0.);
812 gslpp::vector<gslpp::complex> M10i(15, 0.);
813
814 // M_7: LO -> LO_QED
815 for (unsigned int i = LO; i <= int_qed(LO_QED); i++)
816 M_7.push_back(M7i);
817
818 // M_7: NLO_QED11
819 M7i.assign(6, aletilde);
820 M_7.push_back(M7i);
821
822 // M_7: NLO_QED21
823 M7i.reset(); // To clear M7i
824 M7i.assign(0, -alstilde * aletilde * F17(sh));
825 M7i.assign(1, -alstilde * aletilde * F27(sh));
826 M7i.assign(7, -alstilde * aletilde * F87(sh));
827 M_7.push_back(M7i);
828
829 // M_9: LO
830 M9i.assign(8, 1.);
831 M_9.push_back(M9i);
832
833 // M_9: NLO -> LO_QED
834 M9i.reset(); // To clear M9i
835 for (unsigned int i = NLO; i <= int_qed(LO_QED); i++)
836 M_9.push_back(M9i);
837
838 // M_9: NLO_QED11
839 M9i.assign(0, aletilde * f_Huber(sh, -32./27., 4./3., 0., 0., -16./27.));
840 M9i.assign(1, aletilde * f_Huber(sh, -8./9., 1., 0., 0., -4./9.));
841 M9i.assign(2, aletilde * f_Huber(sh, -16./9., 6., -7./2., 2./9., 2./27.));
842 M9i.assign(3, aletilde * f_Huber(sh, 32./27., 0., -2./3., 8./27., 8./81.));
843 M9i.assign(4, aletilde * f_Huber(sh, -112./9., 60., -38., 32./9., -136./27.));
844 M9i.assign(5, aletilde * f_Huber(sh, 512./27., 0., -32./3., 128./27., 320./81.));
845 M9i.assign(10, aletilde * f_Huber(sh, -272./27., 4., 7./6., -74./27., 358./81.));
846 M9i.assign(11, aletilde * f_Huber(sh, -32./81., 0., 2./9., -8./81., -8./243.));
847 M9i.assign(12, aletilde * f_Huber(sh, -2768./27., 40., 38./3., -752./27., 1144./81.));
848 M9i.assign(13, aletilde * f_Huber(sh, -512./81., 0., 32./9., -128./81., -320./243.));
849 M9i.assign(14, aletilde * f_Huber(sh, 16./9., 0., -2., 0., 26./27.));
850 M9i.assign(8, aletilde * f9pen_Huber(sh));
851 M_9.push_back(M9i);
852
853 // M_9: NLO_QED21
854 M9i.reset(); // To clear M9i
855 M9i.assign(0, -alstilde * aletilde * F19(sh));
856 M9i.assign(1, -alstilde * aletilde * F29(sh));
857 M9i.assign(7, -alstilde * aletilde * F89(sh));
858 M_9.push_back(M9i);
859
860 // M_10: LO
861 M10i.assign(9, 1.);
862 M_10.push_back(M10i);
863
864 // M_10: NLO -> NLO_QED21
865 M10i.reset(); // To clear M7i
866 for (unsigned int i = NLO; i <= int_qed(NLO_QED21); i++)
867 M_10.push_back(M10i);
868}
869
870double BXqll::S77_T(double sh, orders order)
871{
872 double umsh = 1. - sh;
873 double sigma = 8.*umsh*umsh/sh;
874 double chi_1 = 4.*umsh*(5.*sh + 3.)/3./sh;
875 double chi_2 = 4.*(3.*sh*sh + 2.*sh - 9.)/sh;
876 double deltaMb2 = (lambda_1*chi_1 + lambda_2*chi_2)/Mb_pole/Mb_pole; // pole mass as stated in hep-ph/9801456
877
878 switch(order)
879 {
880 case LO:
881 return sigma + deltaMb2;
882 case NLO:
883 return sigma*8.*alstilde*omega77_T(sh);
884 default:
885 throw std::runtime_error("BXqll::S77_T: order not implemented");
886 }
887}
888
889double BXqll::S79_T(double sh, orders order)
890{
891 double umsh = 1. - sh;
892 double sigma = 8.*umsh*umsh;
893 double chi_1 = 4.*umsh*umsh;
894 double chi_2 = 4.*(9.*sh*sh - 6.*sh - 7.);
895 double deltaMb2 = (lambda_1*chi_1 + lambda_2*chi_2)/Mb_pole/Mb_pole;
896
897 switch(order)
898 {
899 case LO:
900 return sigma + deltaMb2;
901 case NLO:
902 return sigma*8.*alstilde*omega79_T(sh);
903 default:
904 throw std::runtime_error("BXqll::S79_T: order not implemented");
905 }
906}
907
908double BXqll::S99_T(double sh, orders order)
909{
910 double umsh = 1. - sh;
911 double sigma = 2.*sh*umsh*umsh;
912 double chi_1 = -sh*umsh*(3.*sh + 5.)/3.;
913 double chi_2 = sh*(15.*sh*sh - 14.*sh - 5.);
914 double deltaMb2 = (lambda_1*chi_1 + lambda_2*chi_2)/Mb_pole/Mb_pole;
915 double omega_2 = mySM.Beta0(5.)*log(muh)*omega99_T(sh) + 54.919*umsh*umsh*umsh*umsh -
916 136.374*umsh*umsh*umsh + 119.344*umsh*umsh - 15.6175*umsh - 31.1706;
917
918 switch(order)
919 {
920 case LO:
921 return sigma + deltaMb2;
922 case NLO:
923 return sigma*8.*alstilde*omega99_T(sh);
924 case NNLO:
925 return sigma*16.*alstilde*alstilde*omega_2;
926 default:
927 throw std::runtime_error("BXqll::S99_T: order not implemented");
928 }
929}
930
931double BXqll::S1010_T(double sh, orders order)
932{
933 return S99_T(sh,order);
934}
935
936double BXqll::S77_L(double sh, orders order)
937{
938 double umsh = 1. - sh;
939 double sigma = 4.*umsh*umsh;
940 double chi_1 = -2.*umsh*(3.*sh + 13.)/3.;
941 double chi_2 = 2.*(15.*sh*sh - 6.*sh - 13.);
942 double deltaMb2 = (lambda_1*chi_1 + lambda_2*chi_2)/Mb_pole/Mb_pole;
943
944 switch(order)
945 {
946 case LO:
947 return sigma + deltaMb2;
948 case NLO:
949 return sigma*8.*alstilde*omega77_L(sh);
950 default:
951 throw std::runtime_error("BXqll::S77_L: order not implemented");
952 }
953}
954
955double BXqll::S79_L(double sh, orders order)
956{
957 double umsh = 1. - sh;
958 double sigma = 4.*umsh*umsh;
959 double chi_1 = 2.*umsh*umsh;
960 double chi_2 = 2.*(3.*sh*sh - 6.*sh - 1.);
961 double deltaMb2 = (lambda_1*chi_1 + lambda_2*chi_2)/Mb_pole/Mb_pole;
962
963 switch(order)
964 {
965 case LO:
966 return sigma + deltaMb2;
967 case NLO:
968 return sigma*8.*alstilde*omega79_L(sh);
969 default:
970 throw std::runtime_error("BXqll::S79_L: order not implemented");
971 }
972}
973
974double BXqll::S99_L(double sh, orders order)
975{
976 double umsh = 1. - sh;
977 double sigma = umsh*umsh;
978 double chi_1 = umsh*(13.*sh + 3.)/6.;
979 double chi_2 = (-17.*sh*sh + 10.*sh + 3.)/2.;
980 double deltaMb2 = (lambda_1*chi_1 + lambda_2*chi_2)/Mb_pole/Mb_pole;
981 double omega_2 = mySM.Beta0(5.)*log(muh)*omega99_L(sh) - 5.95974*umsh*umsh*umsh +
982 11.7493*umsh*umsh + 12.2293*umsh - 38.6457;
983
984 switch(order)
985 {
986 case LO:
987 return sigma + deltaMb2;
988 case NLO:
989 return sigma*8.*alstilde*omega99_L(sh);
990 case NNLO:
991 return sigma*16.*alstilde*alstilde*omega_2;
992 default:
993 throw std::runtime_error("BXqll::S99_L: order not implemented");
994 }
995}
996
997double BXqll::S1010_L(double sh, orders order)
998{
999 return S99_L(sh,order);
1000}
1001
1002double BXqll::S710_A(double sh, orders order)
1003{
1004 double umsh = 1. - sh;
1005 double sigma = -8.*umsh*umsh;
1006 double chi_1 = -4.*(3.*sh*sh + 2.*sh + 3.)/3.;
1007 double chi_2 = -4.*(9.*sh*sh - 10.*sh - 7.);
1008 double deltaMb2 = (lambda_1*chi_1 + lambda_2*chi_2)/Mb_pole/Mb_pole;
1009
1010 switch(order)
1011 {
1012 case LO:
1013 return sigma + deltaMb2;
1014 case NLO:
1015 return sigma*8.*alstilde*omega710_A(sh);
1016 default:
1017 throw std::runtime_error("BXqll::S710_A: order not implemented");
1018 }
1019}
1020
1021double BXqll::S910_A(double sh, orders order)
1022{
1023 double umsh = 1. - sh;
1024 double sigma = -4.*sh*umsh*umsh;
1025 double chi_1 = -2.*sh*(3.*sh*sh + 2.*sh + 3.)/3.;
1026 double chi_2 = -2.*sh*(15.*sh*sh - 14.*sh - 9.);
1027 double deltaMb2 = (lambda_1*chi_1 + lambda_2*chi_2)/Mb_pole/Mb_pole;
1028 double omega_2 = mySM.Beta0(5.)*log(muh)*omega910_A(sh) + 74.3717*umsh*umsh*umsh*umsh -
1029 183.885*umsh*umsh*umsh + 158.739*umsh*umsh - 29.0124*umsh - 30.8056;
1030
1031 switch(order)
1032 {
1033 case LO:
1034 return sigma + deltaMb2;
1035 case NLO:
1036 return sigma*8.*alstilde*omega910_A(sh);
1037 case NNLO:
1038 return sigma*16.*alstilde*alstilde*omega_2;
1039 default:
1040 throw std::runtime_error("BXqll::S910_A: order not implemented");
1041 }
1042}
1043
1044gslpp::complex BXqll::cij_T(unsigned int i, unsigned int j, double sh, unsigned int ord)
1045{
1046 double umsh = 1. - sh, uptsh = 1. + 3.*sh;
1047 double r = sh * Mb_pole * Mb_pole / 4. / Mc_pole / Mc_pole;
1048 double pre = aletilde * 8. * lambda_2 / 9. / Mc_pole / Mc_pole;
1049 gslpp::complex F_M7_M9;
1050
1051 switch(ord)
1052 {
1053 case 11:
1054 F_M7_M9 = F_BIR(r) * (M_7[LO](j) / sh + M_9[LO](j) / 2.).conjugate();
1055 break;
1056 case 22:
1057 F_M7_M9 = F_BIR(r) * (M_7[int_qed(NLO_QED11)](j) / sh + M_9[int_qed(NLO_QED11)](j) / 2.).conjugate();
1058 if (i == 0 && j == 1) {
1059 F_M7_M9 = F_BIR(r).conjugate() * (M_7[int_qed(NLO_QED11)](i) / sh +
1060 M_9[int_qed(NLO_QED11)](i) / 2.) - F_M7_M9 / 6.;
1061 }
1062 break;
1063 case 32:
1064 F_M7_M9 = F_BIR(r) * (M_7[int_qed(NLO_QED21)](j) / sh + M_9[int_qed(NLO_QED21)](j) / 2.).conjugate();
1065 if (i == 0 && j == 1) {
1066 F_M7_M9 = F_BIR(r).conjugate() * (M_7[int_qed(NLO_QED21)](i) / sh +
1067 M_9[int_qed(NLO_QED21)](i) / 2.) - F_M7_M9 / 6.;
1068 }
1069 break;
1070 default:
1071 throw std::runtime_error("BXqll::cij_T: order not implemented");
1072 }
1073
1074 if ((i == 1 && j >= i) || 10*(i+1) + j+1 == 12)
1075 return (- pre * umsh * umsh * uptsh * F_M7_M9);
1076 else if (i == 0)
1077 return (pre / 6. * umsh * umsh * uptsh * F_M7_M9);
1078 else
1079 return (0.);
1080}
1081
1082gslpp::complex BXqll::cij_L(unsigned int i, unsigned int j, double sh, unsigned int ord)
1083{
1084 double umsh = 1. - sh, tmsh = 3. - sh;
1085 double r = sh * Mb_pole * Mb_pole / 4. / Mc_pole / Mc_pole;
1086 double pre = aletilde * 8. * lambda_2 / 9. / Mc_pole / Mc_pole;
1087 gslpp::complex F_M7_M9;
1088
1089 switch(ord)
1090 {
1091 case 11:
1092 F_M7_M9 = F_BIR(r) * (M_7[LO](j) + M_9[LO](j) / 2.).conjugate();
1093 break;
1094 case 22:
1095 F_M7_M9 = F_BIR(r) * (M_7[int_qed(NLO_QED11)](j) + M_9[int_qed(NLO_QED11)](j) / 2.).conjugate();
1096 if (i == 0 && j == 1) {
1097 F_M7_M9 = F_BIR(r).conjugate() * (M_7[int_qed(NLO_QED11)](i) +
1098 M_9[int_qed(NLO_QED11)](i) / 2.) - F_M7_M9 / 6.;
1099 }
1100 break;
1101 case 32:
1102 F_M7_M9 = F_BIR(r) * (M_7[int_qed(NLO_QED21)](j) + M_9[int_qed(NLO_QED21)](j) / 2.).conjugate();
1103 if (i == 0 && j == 1) {
1104 F_M7_M9 = F_BIR(r).conjugate() * (M_7[int_qed(NLO_QED21)](i) +
1105 M_9[int_qed(NLO_QED21)](i) / 2.) - F_M7_M9 / 6.;
1106 }
1107 break;
1108 default:
1109 throw std::runtime_error("BXqll::cij_L: order not implemented");
1110 }
1111
1112 if ((i == 1 && j >= i) || 10*(i+1) + j+1 == 12)
1113 return (- pre * umsh * umsh * tmsh * F_M7_M9);
1114 else if (i == 0)
1115 return (pre / 6. * umsh * umsh * tmsh * F_M7_M9);
1116 else
1117 return (0.);
1118}
1119
1120gslpp::complex BXqll::cij_A(unsigned int i, unsigned int j, double sh)
1121{
1122 unsigned int ij = 100*(i + 1) + (j + 1);
1123 double umsh = 1. - sh, uptsh = 1. + 3.*sh;
1124 double r = sh * Mb_pole * Mb_pole / 4. / Mc_pole / Mc_pole;
1125
1126 switch(ij)
1127 {
1128 case 110:
1129 return (-aletilde*4.*lambda_2/54./Mc_pole/Mc_pole*umsh*umsh*uptsh * F_BIR(r));
1130 case 210:
1131 return (aletilde*4.*lambda_2/9./Mc_pole/Mc_pole*umsh*umsh*uptsh * F_BIR(r));
1132 default:
1133 return (0.);
1134 }
1135}
1136
1137gslpp::complex BXqll::eij_T(unsigned int i, unsigned int j, double sh)
1138{
1139 unsigned int ij;
1140 double umsh = (1. - sh);
1141 double sigma77_T = 8.*umsh*umsh/sh, sigma79_T = 8.*umsh*umsh, sigma99_T = 2.*sh*umsh*umsh;
1142
1143 if (j == 9)
1144 ij = 100*(i + 1) + (j + 1);
1145 else
1146 ij = 10*(i + 1) + (j + 1);
1147
1148 switch(ij)
1149 {
1150 case 11:
1151 return (128./9. * aletilde*aletilde*aletilde * sigma99_T*omega22em_T(sh));
1152 case 12:
1153 return (64./3. * aletilde*aletilde*aletilde * sigma99_T*omega22em_T(sh));
1154 case 17:
1155 return (32./3. * aletilde*aletilde*aletilde * sigma79_T*omega27em_T(sh));
1156 case 19:
1157 return (32./3. * aletilde*aletilde * sigma99_T*omega29em_T(sh));
1158 case 22:
1159 return (8. * aletilde*aletilde*aletilde * sigma99_T*omega22em_T(sh));
1160 case 27:
1161 return (8. * aletilde*aletilde*aletilde * sigma79_T*omega27em_T(sh));
1162 case 29:
1163 return (8. * aletilde*aletilde * sigma99_T*omega29em_T(sh));
1164 case 77:
1165 return (8. * aletilde*aletilde*aletilde * sigma77_T*omega77em_T(sh));
1166 case 79:
1167 return (8. * aletilde*aletilde * sigma79_T*omega79em_T(sh));
1168 case 99:
1169 return (8. * aletilde * sigma99_T*omega99em_T(sh));
1170 case 1010:
1171 return (8. * aletilde * sigma99_T*omega99em_T(sh));
1172 default:
1173 return (0.);
1174 }
1175}
1176
1177gslpp::complex BXqll::eij_L(unsigned int i, unsigned int j, double sh)
1178{
1179 unsigned int ij;
1180 double umsh = (1. - sh);
1181 double sigma77_L = 4.*umsh*umsh, sigma79_L = 4.*umsh*umsh, sigma99_L = umsh*umsh;
1182
1183 if (j == 9)
1184 ij = 100* (i + 1) + (j + 1);
1185 else
1186 ij = 10*(i + 1) + (j + 1);
1187
1188 switch(ij)
1189 {
1190 case 11:
1191 return (128./9. * aletilde*aletilde*aletilde * sigma99_L*omega22em_L(sh));
1192 case 12:
1193 return (64./3. * aletilde*aletilde*aletilde * sigma99_L*omega22em_L(sh));
1194 case 17:
1195 return (32./3. * aletilde*aletilde*aletilde * sigma79_L*omega27em_L(sh));
1196 case 19:
1197 return (32./3. * aletilde*aletilde * sigma99_L*omega29em_L(sh));
1198 case 22:
1199 return (8. * aletilde*aletilde*aletilde * sigma99_L*omega22em_L(sh));
1200 case 27:
1201 return (8. * aletilde*aletilde*aletilde * sigma79_L*omega27em_L(sh));
1202 case 29:
1203 return (8. * aletilde*aletilde * sigma99_L*omega29em_L(sh));
1204 case 77:
1205 return (8. * aletilde*aletilde*aletilde * sigma77_L*omega77em_L(sh));
1206 case 79:
1207 return (8. * aletilde*aletilde * sigma79_L*omega79em_L(sh));
1208 case 99:
1209 return (8. * aletilde * sigma99_L*omega99em_L(sh));
1210 case 1010:
1211 return (8. * aletilde * sigma99_L*omega99em_L(sh));
1212 default:
1213 return (0.);
1214 }
1215}
1216
1217gslpp::complex BXqll::eij_A(unsigned int i, unsigned int j, double sh)
1218{
1219 unsigned int ij;
1220 double umsh = (1. - sh);
1221 double sigma710_A = -8.*umsh*umsh, sigma910_A = -4.*sh*umsh*umsh;
1222
1223 if (j == 9)
1224 ij = 100*(i + 1) + (j + 1);
1225 else
1226 ij = 10*(i + 1) + (j + 1);
1227
1228 switch(ij)
1229 {
1230 case 110:
1231 return (32./3. * aletilde*aletilde * sigma910_A*omega210em_A(sh));
1232 case 210:
1233 return (8. * aletilde*aletilde * sigma910_A*omega210em_A(sh));
1234 case 710:
1235 return (8. * aletilde*aletilde * sigma710_A*omega710em_A(sh));
1236 case 910:
1237 return (8. * aletilde * sigma910_A*omega910em_A(sh));
1238 default:
1239 return (0.);
1240 }
1241}
1242
1243double BXqll::omega77_T(double sh)
1244{
1245 double umsh = 1.-sh;
1246 double umsqrt = 1.-sqrt(sh);
1247 double dilog_umsh = gslpp_special_functions::dilog(umsh);
1248 double dilog_umsqrt = gslpp_special_functions::dilog(umsqrt);
1249
1250 return (-8./3.*log(muh) - (sqrt(sh)+1.)*(sqrt(sh)+1.)*(pow(sh,1.5)-10.*sh+13.*sqrt(sh)-8.)*
1251 dilog_umsh/6./umsh/umsh + 2.*sqrt(sh)*(sh*sh-6.*sh-3.)*dilog_umsqrt/3./umsh/umsh -
1252 M_PI*M_PI*(3.*pow(sh,1.5)+22.*sh+23.*sqrt(sh)+16.)*umsqrt*umsqrt/36./umsh/umsh +
1253 (5.*sh*sh*sh-54.*sh*sh+57.*sh-8.)/18./umsh/umsh - log(umsh) +
1254 sh*(5.*sh+1.)*log(sh)/3./umsh/umsh + 2./3.*log(umsh)*log(sh));
1255}
1256
1257double BXqll::omega79_T(double sh)
1258{
1259 double umsh = 1.-sh;
1260 double umsqrt = 1.-sqrt(sh);
1261 double dilog_umsh = gslpp_special_functions::dilog(umsh);
1262 double dilog_umsqrt = gslpp_special_functions::dilog(umsqrt);
1263
1264 return (-4./3.*log(muh) - 2.*sqrt(sh)*(sh+3.)*dilog_umsqrt/3./umsh/umsh - M_PI*M_PI*
1265 (16.*sh+29.*sqrt(sh)+19.)*umsqrt*umsqrt/36./umsh/umsh + (sh*sh-6.*sh+5.)/6./umsh/umsh +
1266 (sqrt(sh)+1.)*(sqrt(sh)+1.)*(8.*sh-15.*sqrt(sh)+9.)*dilog_umsh/6./umsh/umsh - (5.*sh+1.)*
1267 log(umsh)/6./sh + sh*(3.*sh+1.)*log(sh)/6./umsh/umsh + 2./3.*log(umsh)*log(sh));
1268}
1269
1270double BXqll::omega99_T(double sh)
1271{
1272 double umsh = 1.-sh;
1273 double umsqrt = 1.-sqrt(sh);
1274 double dilog_umsh = gslpp_special_functions::dilog(umsh);
1275 double dilog_umsqrt = gslpp_special_functions::dilog(umsqrt);
1276
1277 return ((sqrt(sh)+1.)*(sqrt(sh)+1.)*(8.*pow(sh,1.5)-15.*sh+4.*sqrt(sh)-5.)*dilog_umsh/
1278 6./umsh/umsh/sqrt(sh) - 2.*(sh*sh-12.*sh-5.)*dilog_umsqrt/3./umsh/umsh/sqrt(sh) -
1279 M_PI*M_PI*(16.*pow(sh,1.5)+29.*sh+4.*sqrt(sh)+15.)*umsqrt*umsqrt/36./umsh/umsh/sqrt(sh) +
1280 (2.*sh*sh-7.*sh-5.)*log(sh)/3./umsh/umsh + (sh*sh+18.*sh-19.)/6./umsh/umsh -
1281 (2.*sh+1)*log(umsh)/3./sh + 2./3.*log(umsh)*log(sh));
1282}
1283
1284double BXqll::omega77_L(double sh)
1285{
1286 double umsh = 1.-sh;
1287 double umsqrt = 1.-sqrt(sh);
1288 double dilog_umsh = gslpp_special_functions::dilog(umsh);
1289 double dilog_umsqrt = gslpp_special_functions::dilog(umsqrt);
1290
1291 return (-8./3.*log(muh) + (sqrt(sh)+1.)*(sqrt(sh)+1.)*(4.*pow(sh,1.5)-7.*sh+2.*sqrt(sh)-3.)*
1292 dilog_umsh/3./umsh/umsh/sqrt(sh) - (9.*sh*sh-38.*sh+29.)/6./umsh/umsh -
1293 4.*(sh*sh-6.*sh-3.)*dilog_umsqrt/3./umsh/umsh/sqrt(sh) - M_PI*M_PI*
1294 (8.*pow(sh,1.5)+13.*sh+2.*sqrt(sh)+9.)*umsqrt*umsqrt/18./umsh/umsh/sqrt(sh) - (sh*sh*sh-3.*sh+2.)*
1295 log(umsh)/3./umsh/umsh/sh + 2.*(sh*sh-3.*sh-3.)*log(sh)/3./umsh/umsh + 2./3.*log(umsh)*log(sh));
1296}
1297
1298double BXqll::omega79_L(double sh)
1299{
1300 double umsh = 1.-sh;
1301 double umsqrt = 1.-sqrt(sh);
1302 double dilog_umsh = gslpp_special_functions::dilog(umsh);
1303 double dilog_umsqrt = gslpp_special_functions::dilog(umsqrt);
1304
1305 return (-4./3.*log(muh) + 4.*sqrt(sh)*(sh+3.)*dilog_umsqrt/3./umsh/umsh + (sqrt(sh)+1.)*
1306 (sqrt(sh)+1.)*(4.*sh-9.*sqrt(sh)+3.)*dilog_umsh/3./umsh/umsh + (7.*sh*sh-2.*sh-5.)/
1307 6./umsh/umsh - M_PI*M_PI*(8.*sh+19.*sqrt(sh)+5.)*umsqrt*umsqrt/18./umsh/umsh -
1308 (2.*sh+1.)*log(umsh)/3./sh + (sh-7.)*sh*log(sh)/3./umsh/umsh + 2./3.*log(umsh)*log(sh));
1309}
1310
1311double BXqll::omega99_L(double sh)
1312{
1313 double umsh = 1.-sh;
1314 double umsqrt = 1.-sqrt(sh);
1315 double dilog_umsh = gslpp_special_functions::dilog(umsh);
1316 double dilog_umsqrt = gslpp_special_functions::dilog(umsqrt);
1317
1318 return (-(sqrt(sh)+1.)*(sqrt(sh)+1.)*(pow(sh,1.5)-8.*sh+3.*sqrt(sh)-4.)*dilog_umsh/
1319 3./umsh/umsh + 4.*sqrt(sh)*(sh*sh-12.*sh-5.)*dilog_umsqrt/3./umsh/umsh -
1320 M_PI*M_PI*(3.*pow(sh,1.5)+20.*sh+sqrt(sh)+8.)*umsqrt*umsqrt/18./umsh/umsh +
1321 (4.*sh*sh*sh-51.*sh*sh+42.*sh+5.)/6./umsh/umsh - log(umsh) +
1322 8.*sh*(2.*sh+1.)*log(sh)/3./umsh/umsh + + 2./3.*log(umsh)*log(sh));
1323}
1324
1325double BXqll::omega710_A(double sh)
1326{
1327 double umsh = 1.-sh;
1328 double umsqrt = 1.-sqrt(sh);
1329
1330 // Expression from 1503.04849, which appears to be wrong
1331// double num = 3.*umsh*umsh;
1332// double dilog_umsh = gslpp_special_functions::dilog(umsh);
1333// double dilog_umsqrt = gslpp_special_functions::dilog(umsqrt);
1334//
1335// return (-4./3.*log(muh) + 2.*(4.*sh*sh-13.*sh-1.)*dilog_umsqrt/num - (2.*sh*sh-9.*sh-3.)*dilog_umsh/num -
1336// (3.*sh*sh-16.*sh+13.)*log(umsqrt)/num + (4.*sh*sh-13.*sh-1.)*log(umsqrt)*log(sh)/num -
1337// (2.*sh*sh-9.*sh-3.)*log(umsh)*log(sh)/num + (sh*sh*sh-23.*sh*sh+23.*sh-1.)*log(umsh)/2./num/sh +
1338// (sh-20.*sqrt(sh)+5.)*umsqrt*umsqrt/2./num - M_PI*M_PI/3.);
1339
1340 // Expression from hep-ph/0612156
1341 double upsqrt = 1. + sqrt(sh);
1342 double tdsqrt = 3. - 2. * sqrt(sh);
1343 double dilog_sh = gslpp_special_functions::dilog(sh);
1344 double dilog_sqrt = gslpp_special_functions::dilog(sqrt(sh));
1345
1346 return ((3. + sh*(9. - 2.*sh)) * dilog_sh / umsh / umsh + (1. - sh*(22. - sh)) * log(umsh) / 2. / sh / umsh -
1347 2.*(1. + sh*(13. - 4.*sh)) * dilog_sqrt / umsh / umsh + (13. - 3.*sh) * log(umsqrt) / umsh +
1348 5.*(1. + sh)*M_PI*M_PI / 6. / umsh / umsh - (sh + tdsqrt*tdsqrt) / 2. / upsqrt / upsqrt)*(-1./3.);
1349}
1350
1351double BXqll::omega910_A(double sh)
1352{
1353 double umsh = 1.-sh;
1354 double num = 3.*umsh*umsh;
1355 double umsqrt = 1.-sqrt(sh);
1356 double dilog_umsh = gslpp_special_functions::dilog(umsh);
1357 double dilog_umsqrt = gslpp_special_functions::dilog(umsqrt);
1358
1359 return (-2.*(sh*sh-3.*sh-1.)*dilog_umsh/num - 4.*(5.-2.*sh)*sh*dilog_umsqrt/num -
1360 (4.*sqrt(sh)-3.)*umsqrt*umsqrt/num - 2.*(2.*sh*sh-7.*sh+5.)*log(umsqrt)/num -
1361 2.*(sh*sh-3.*sh-1.)*log(umsh)*log(sh)/num + (2.*sh*sh*sh-11.*sh*sh+10.*sh-1.)*
1362 log(umsh)/num/sh + 2.*sh*(2.*sh-5.)*log(umsqrt)*log(sh)/num - M_PI*M_PI/3.);
1363}
1364
1365double BXqll::omega77em_T(double sh)
1366{
1367 double umsh = 1. - sh;
1368
1369 return (Lbl*(1.54986 - 1703.72*sh*sh*sh*sh*sh + 1653.38*sh*sh*sh*sh - 683.608*sh*sh*sh +
1370 179.279*sh*sh - 35.5047*sh)/8./umsh/umsh);
1371}
1372
1373double BXqll::omega79em_T(double sh)
1374{
1375 double umsh = 1. - sh;
1376
1377 return (Lbl*(19.063 + 2158.03*sh*sh*sh*sh - 2062.92*sh*sh*sh + 830.53*sh*sh -
1378 186.12*sh + 0.324236/sh)/8./umsh/umsh);
1379}
1380
1381double BXqll::omega99em_T(double sh)
1382{
1383 double umsh = 1. - sh;
1384
1385 return (Lbl*(2.2596 + 157.984*sh*sh*sh*sh - 141.281*sh*sh*sh + 52.8914*sh*sh -
1386 13.5377*sh + 0.0284049/sh)/2./sh/umsh/umsh);
1387}
1388
1389double BXqll::omega22em_T(double sh)
1390{
1391 double umsh = 1. - sh;
1392 double Lmub = log(mu_b/5.);
1393
1394 return (Lbl*((2.84257 + 269.974*sh*sh*sh*sh - 194.443*sh*sh*sh + 48.4535*sh*sh -
1395 8.24929*sh + 0.0111118/sh)/2./sh/umsh/umsh +
1396 Lmub*4.*(4.54727 + 330.182*sh*sh*sh*sh - 258.194*sh*sh*sh + 79.8713*sh*sh -
1397 19.6855*sh + 0.0371348/sh)/9./sh/umsh/umsh) +
1398 64./81.*omega99em_T(sh)*Lmub*Lmub);
1399}
1400
1401gslpp::complex BXqll::omega27em_T(double sh)
1402{
1403 double umsh = 1. - sh;
1404 double Lmub = log(mu_b/5.);
1405 gslpp::complex i = gslpp::complex::i();
1406
1407 return (Lbl*((21.5291 + 3044.94*sh*sh*sh*sh - 2563.05*sh*sh*sh + 874.074*sh*sh -
1408 175.874*sh + 0.121398/sh)/8./umsh/umsh +
1409 i*(2.49475 + 598.376*sh*sh*sh*sh - 456.831*sh*sh*sh + 117.683*sh*sh -
1410 9.90525*sh - 0.0116501/sh)/8./umsh/umsh) +
1411 8./9.*omega79em_T(sh)*Lmub);
1412}
1413
1414gslpp::complex BXqll::omega29em_T(double sh)
1415{
1416 double umsh = 1. - sh;
1417 double Lmub = log(mu_b/5.);
1418 gslpp::complex i = gslpp::complex::i();
1419
1420 return (Lbl*((4.54727 + 330.182*sh*sh*sh*sh - 258.194*sh*sh*sh + 79.8713*sh*sh -
1421 19.6855*sh + 0.0371348/sh)/2./sh/umsh/umsh +
1422 i*(73.9149*sh*sh*sh*sh - 61.1338*sh*sh*sh + 14.6517*sh*sh - 0.102331*sh +
1423 0.710037)/2./sh/umsh/umsh) +
1424 16./9.*omega99em_T(sh)*Lmub);
1425}
1426
1427double BXqll::omega77em_L(double sh)
1428{
1429 double umsh = 1. - sh;
1430
1431 return (Lbl*(9.73761 + 647.747*sh*sh*sh*sh - 642.637*sh*sh*sh + 276.839*sh*sh -
1432 68.3562*sh - 1.6755/sh)/4./umsh/umsh);
1433}
1434
1435double BXqll::omega79em_L(double sh)
1436{
1437 double umsh = 1. - sh;
1438
1439 return (Lbl*(-6.03641 - 896.643*sh*sh*sh*sh + 807.349*sh*sh*sh - 278.559*sh*sh +
1440 47.6636*sh - 0.190701/sh)/4./umsh/umsh);
1441}
1442
1443double BXqll::omega99em_L(double sh)
1444{
1445 double umsh = 1. - sh;
1446
1447 return (Lbl*(-0.768521 - 80.8068*sh*sh*sh*sh + 70.0821*sh*sh*sh - 21.2787*sh*sh +
1448 2.9335*sh - 0.0180809/sh)/umsh/umsh);
1449}
1450
1451double BXqll::omega22em_L(double sh)
1452{
1453 double umsh = 1. - sh;
1454 double Lmub = log(mu_b/5.);
1455
1456 return (Lbl*((-1.71832 - 234.11*sh*sh*sh*sh + 162.126*sh*sh*sh - 37.2361*sh*sh +
1457 6.29949*sh - 0.00810233/sh)/umsh/umsh +
1458 Lmub*8.*(224.662*sh*sh*sh - 2.27221 - 298.369*sh*sh*sh*sh - 65.1375*sh*sh +
1459 11.5686*sh - 0.0233098/sh)/9./umsh/umsh) +
1460 64./81.*omega99em_L(sh)*Lmub*Lmub);
1461}
1462
1463gslpp::complex BXqll::omega27em_L(double sh)
1464{
1465 double umsh = 1. - sh;
1466 double Lmub = log(mu_b/5.);
1467 gslpp::complex i = gslpp::complex::i();
1468
1469 return (Lbl*((-8.01684 - 1121.13*sh*sh*sh*sh + 882.711*sh*sh*sh - 280.866*sh*sh +
1470 54.1943*sh - 0.128988/sh)/4./umsh/umsh +
1471 i*(-2.14058 - 588.771*sh*sh*sh*sh + 483.997*sh*sh*sh - 124.579*sh*sh +
1472 12.3282*sh + 0.0145059/sh)/4./umsh/umsh) +
1473 8./9.*omega79em_L(sh)*Lmub);
1474}
1475
1476gslpp::complex BXqll::omega29em_L(double sh)
1477{
1478 double umsh = 1. - sh;
1479 double Lmub = log(mu_b/5.);
1480 gslpp::complex i = gslpp::complex::i();
1481
1482 return (Lbl*((-2.27221 - 298.369*sh*sh*sh*sh + 224.662*sh*sh*sh - 65.1375*sh*sh +
1483 11.5686*sh - 0.0233098/sh)/umsh/umsh +
1484 i*(-0.666157 - 120.303*sh*sh*sh*sh + 109.315*sh*sh*sh - 28.2734*sh*sh +
1485 2.44527*sh + 0.00279781/sh)/umsh/umsh) +
1486 16./9.*omega99em_L(sh)*Lmub);
1487}
1488
1489double BXqll::omega710em_A(double sh)
1490{
1491 double umsh = 1. - sh;
1492
1493 return (Lbl*((7. - 16.*sqrt(sh) + 9.*sh)/4./umsh + log(1. - sqrt(sh)) +
1494 (1. + 3.*sh)/umsh*log((1. + sqrt(sh))/2.) - sh*log(sh)/umsh));
1495}
1496
1497double BXqll::omega910em_A(double sh)
1498{
1499 double umsh = 1. - sh;
1500
1501 return (Lbl*(log(1. - sqrt(sh)) - (5. - 16.*sqrt(sh) + 11.*sh)/4./umsh +
1502 (1. - 5.*sh)/umsh*log((1. + sqrt(sh))/2.) - (1. - 3.*sh)*log(sh)/umsh));
1503}
1504
1505gslpp::complex BXqll::omega210em_A(double sh)
1506{
1507 double umsh = 1. - sh;
1508 double a = 16.*Mc*Mc*Mc*Mc/Mb/Mb/Mb/Mb;
1509 double shma = sh - a;
1510 double Lmub = log(mu_b/5.);
1511 gslpp::complex omega, i = gslpp::complex::i();
1512
1513 omega = Lbl*((-351.322*sh*sh*sh*sh + 378.173*sh*sh*sh - 160.158*sh*sh + 24.2096*sh +
1514 0.305176)/24./sh/umsh/umsh) +
1515 8./9.*omega910em_A(sh)*Lmub;
1516 if(sh > a)
1517 omega += Lbl*i*(7.98625 + 238.507*shma - 766.869*shma*shma)*shma*shma/24./sh/umsh/umsh;
1518
1519 return omega;
1520}
1521
1522gslpp::complex BXqll::f_Huber(double sh, double gamma_9, double rho_c, double rho_b, double rho_0, double rho_num)
1523{
1524 gslpp::complex i = gslpp::complex::i();
1525
1526// return (gamma_9 * log(Mb/mu_b) + rho_c * (g_Huber(4.*Mc_pole*Mc_pole/Mb_pole/Mb_pole/sh) +
1527// 8./9. * log(Mb/Mc)) + rho_b * g_Huber(4./sh) +
1528// rho_0 * (log(sh) - i*M_PI) + rho_num);
1529
1530 return (-gamma_9 * log(muh) + rho_c * (g_Huber(4.*z/sh) - 4./9. * log(z) + KS_cc(sh)) +
1531 rho_b * g_Huber(4./sh) + rho_0 * (log(sh) - i*M_PI) + rho_num);
1532}
1533
1534gslpp::complex BXqll::f9pen_Huber(double sh)
1535{
1536 gslpp::complex i = gslpp::complex::i();
1537
1538// return (8. * log(Mb/mu_b) - 3. * g_Huber(4.*Mtau*Mtau/Mb_pole/Mb_pole/sh) - 8./3. * log(Mb/Mtau) +
1539// 8./3. * (log(sh) - i*M_PI) - 40./9.);
1540
1541 return (-8. * log(muh) - 3. * g_Huber(4.*Mtau*Mtau/Mb_pole/Mb_pole/sh) - 8./3. * log(Mb_pole/Mtau) +
1542 8./3. * (log(sh) - i*M_PI) - 40./9.);
1543}
1544
1545gslpp::complex BXqll::g_Huber(double y)
1546{
1547 gslpp::complex i = gslpp::complex::i();
1548 gslpp::complex g_y;
1549
1550 g_y = -2./9.*(2.+y)*sqrt(fabs(1.-y));
1551 if(y < 1.)
1552 g_y *= log(fabs((1.+sqrt(1.-y))/(1.-sqrt(1.-y))))-i*M_PI;
1553 else
1554 g_y *= 2.*atan(1./sqrt(y-1.));
1555
1556 g_y += 20./27. + 4./9.*y;
1557
1558 return (g_y);
1559}
1560
1561gslpp::complex BXqll::KS_cc(double sh)
1562{
1563 double pre = 3. * sh / ale / ale;
1564 gslpp::complex i = gslpp::complex::i();
1565 gslpp::complex kscc;
1566
1567 // Contributions of J/Psi, Psi(2S)
1568 // Masses, total decay widths, and branching fractions taken from Kruger:1996cv
1569 kscc = KS_aux(sh, 3.097, 0.088e-3, 6.0e-2, 0.877);
1570 kscc += KS_aux(sh, 3.686, 0.277e-3, 8.3e-3, 0.9785);
1571
1572 kscc = pre * kscc;
1573
1574 // Contributions from the continuum for the large-q^2 region
1575 // 0.6 < sh < 0.69 region divergent, not implemented
1576 if (sh > 0.69)
1577 kscc += (i * M_PI * 1.02 - log(1. - sh / 4. / z)) / 3.;
1578
1579 return (kscc);
1580}
1581
1582gslpp::complex BXqll::KS_aux(double sh, double m, double Gamma, double Br_ll, double Br_had)
1583{
1584 gslpp::complex i = gslpp::complex::i();
1585 double Gamma_had = Br_had * Gamma / Mb_pole;
1586 double a = m * m / Mb_pole / Mb_pole;
1587 double b = a * Gamma * Gamma / Mb_pole / Mb_pole;
1588 double sqrtb = sqrt(b);
1589 double amsh = a - sh;
1590 double am4z = a - 4. * z;
1591
1592 return (Br_ll * Gamma / Mb_pole * Gamma_had * ((M_PI * amsh + 2. * amsh * atan(am4z / sqrtb) +
1593 sqrtb * (log(b + am4z * am4z) - 2. * log(4. * z - sh))) /
1594 (2. * sqrtb * (b + amsh * amsh)) + i * M_PI / (amsh * amsh + b)));
1595}
1596
1597gslpp::complex BXqll::F_BIR(double r)
1598{
1599 gslpp::complex i = gslpp::complex::i();
1600
1601 if(r > 0. && r < 1.)
1602 return (1.5 / r * (atan(sqrt(r / (1. - r))) / sqrt(r - r*r) - 1.));
1603 else if (r > 1.)
1604 return (1.5 / r * ((log((1. - sqrt(1. - 1./r)) / (1. + sqrt(1. - 1./r))) + i * M_PI) /
1605 2. / sqrt(r*r - r) - 1.));
1606 else
1607 throw std::runtime_error("BXqll::F_BIR(): 1/mc^2 corrections diverge at q^2 = 4*mc^2");
1608}
1609
1610double BXqll::Phi_u(orders ord)
1611{
1612 switch(ord)
1613 {
1614 case LO:
1615 return(1.);
1616 case NLO:
1617 return(alstilde * phi1);
1618 case NNLO:
1619 return(alstilde * alstilde * (phi2 + 2. * mySM.Beta0(5) * phi1 * log(muh))
1620 + (lambda_1 / 2. - 9. / 2. * lambda_2) / Mb_pole / Mb_pole );
1621 case FULLNNLO:
1622 return (1. + alstilde * phi1 +
1623 alstilde * alstilde * (phi2 + 2. * mySM.Beta0(5) * phi1 * log(muh)) +
1624 (lambda_1 / 2. - 9. / 2. * lambda_2) / Mb_pole / Mb_pole);
1625 default:
1626 throw std::runtime_error("BXqll::Phi_u(): order not implemented.");
1627 }
1628}
1629
1630double BXqll::Phi_u(orders_qed ord_qed)
1631{
1632 switch(ord_qed)
1633 {
1634 case LO_QED:
1635 return (kappa * (12. / 23. * (1. - alsmu / mySM.Als(mySM.getMuw(), FULLNNNLO, true, true))));
1636 case FULLNLO_QED:
1637 return (Phi_u(FULLNNLO) +
1638 kappa * (12. / 23. * (1. - alsmu / mySM.Als(mySM.getMuw(), FULLNNNLO, true, true))));
1639 default:
1640 throw std::runtime_error("BXqll::Phi_u(): order not implemented.");
1641 }
1642}
1643
1644double BXqll::Phi_u_inv(unsigned int ord_qcd, unsigned int ord_qed)
1645{
1646 if (ord_qcd == QCD0 && ord_qed == QED0)
1647 return (phinv00);
1648 else if (ord_qcd == QCD1 && ord_qed == QED0)
1649 return (phinv10);
1650 else if (ord_qcd == QCD2 && ord_qed == QED0)
1651 return (phinv20);
1652 else if (ord_qcd == QCD0 && ord_qed == QED1)
1653 return (phinv01);
1654 else if (ord_qcd == QCD1 && ord_qed == QED1)
1655 return (phinv11);
1656 else if (ord_qcd == QCD2 && ord_qed == QED1)
1657 return (phinv21);
1658 else
1659 return (0.);
1660}
1661
1662unsigned int BXqll::int_qed(orders_qed order_qed)
1663{
1664 // For LO_QED to come right after NNLO
1665 return (order_qed - NO_QED + NNLO);
1666}
1667
1668double BXqll::CCH_multiplication(std::vector< gslpp::matrix<gslpp::complex> >& Hij)
1669{
1670 double Phi_ll_u = 0.;
1671
1672 // C_i (qcd_a, qed_a) * C_j (qcd_b, qed_b) * H_ij (qcd_c, qed_c) / Phi_u (qcd_u, qed_u)
1673
1674 for(unsigned int j = 0; j < 15; j++)
1675 for(unsigned int i = 0; i <= j; i++)
1676 for(unsigned int qcd_u = QCD0; qcd_u <= QCD2; qcd_u++)
1677 for(unsigned int qed_u = QED0; qed_u <= QED1; qed_u++)
1678 for(unsigned int qcd_a = QCD0; qcd_a <= QCD2; qcd_a++)
1679 for(unsigned int qed_a = QED0; qed_a <= QED2; qed_a++)
1680 for(unsigned int qcd_b = QCD0; qcd_b <= QCD2; qcd_b++)
1681 for(unsigned int qed_b = QED0; qed_b <= QED2; qed_b++)
1682 {
1683 for(unsigned int qcd_c = QCD0; qcd_c <= QCD2; qcd_c++)
1684 for(unsigned int qed_c = QED0; qed_c <= QED2; qed_c++)
1685 {
1686 if (qcd_a + qcd_b + qcd_c + qcd_u <= QCD_max && qed_a + qed_b + qed_c + qed_u <= QED_max)
1687 Phi_ll_u += (WC(i, qcd_a + 3*qed_a).conjugate() * WC(j, qcd_b + 3*qed_b) *
1688 Hij[qcd_c + 3*qed_c](i,j) ).real() * Phi_u_inv(qcd_u, qed_u);
1689 }
1690
1691 if (qcd_a + qcd_b + 3 + qcd_u <= QCD_max && qed_a + qed_b + 1 + qed_u <= QED_max)
1692 Phi_ll_u += (WC(i, qcd_a + 3*qed_a).conjugate() * WC(j, qcd_b + 3*qed_b) *
1693 Hij[QCD2 + 3*QED2 + 1](i,j) ).real() * Phi_u_inv(qcd_u, qed_u);
1694
1695 if (qcd_a + qcd_b + 3 + qcd_u <= QCD_max && qed_a + qed_b + 2 + qed_u <= QED_max)
1696 Phi_ll_u += (WC(i, qcd_a + 3*qed_a).conjugate() * WC(j, qcd_b + 3*qed_b) *
1697 Hij[QCD2 + 3*QED2 + 2](i,j) ).real() * Phi_u_inv(qcd_u, qed_u);
1698
1699 if (qcd_a + qcd_b + 3 + qcd_u <= QCD_max && qed_a + qed_b + 3 + qed_u <= QED_max)
1700 Phi_ll_u += (WC(i, qcd_a + 3*qed_a).conjugate() * WC(j, qcd_b + 3*qed_b) *
1701 Hij[QCD2 + 3*QED2 + 3](i,j) ).real() * Phi_u_inv(qcd_u, qed_u);
1702 }
1703
1704 return Phi_ll_u;
1705}
1706
1707
1708/***********************************
1709 * Temporary method used for tests *
1710 * TO BE REMOVED *
1711 ***********************************/
1712void BXqll::Test_WC_DF1()
1713{
1714 double muw = mySM.getMuw();
1715// double alstilde_2 = alstilde * alstilde;
1716// double kappa_2 = kappa * kappa;
1717// double GF_2 = GF * GF;
1718// double MW_2 = MW * MW;
1720// double MZ = mySM.getMz();
1721// double sw2 = 1. - (MW_2 / MZ / MZ);
1722// gslpp::complex Vtb = mySM.getCKM().getV_tb();
1723// gslpp::complex Vts = mySM.getCKM().getV_ts();
1724// gslpp::complex lambda_t = Vtb.conjugate() * Vts;
1725// double xt = mySM.getMatching().x_t(muw);
1726// double alemuw = mySM.Ale(muw, FULLNLO);
1727// double aletildemuw = alemuw / 4. / M_PI;
1728//
1729// std::cout << "MW = " << MW << std::endl;
1730// std::cout << "1/ale(muw) = " << 1./alemuw << std::endl;
1731// std::cout << "1/ale(mub) = " << 1./ale << std::endl;
1732// std::cout << "xt = " << xt << std::endl;
1733//
1734// gslpp::complex c7_df1 = 0.;
1735// gslpp::complex c9_df1 = 0.;
1736// gslpp::complex c10_df1 = 0.;
1737//
1738// for (int i = 0; i < 5; i++)
1739// {
1740// c7_df1 += WC(6, i);
1741// c9_df1 += WC(8, i);
1742// c10_df1 += WC(9, i);
1743// }
1744// for (int i = 5; i < 9; i++)
1745// {
1746// c9_df1 += WC(8, i);
1747// c10_df1 += WC(9, i);
1748// }
1749//
1750// std::cout << std::endl;
1751//
1752// std::cout << "00: " << WC(6,0) << std::endl;
1753// std::cout << "10: " << WC(6,1) / alstilde << std::endl;
1754// std::cout << "20: " << WC(6,2) / alstilde / alstilde << std::endl;
1755// std::cout << "01: " << WC(6,3) / kappa << std::endl;
1756// std::cout << "11: " << WC(6,4) / alstilde / kappa << std::endl;
1757// std::cout << "C7: " << c7_df1 << std::endl << std::endl;
1758//
1759// std::cout << "01: " << WC(8,3) / kappa << std::endl;
1760// std::cout << "11: " << WC(8,4) / alstilde / kappa << std::endl;
1761// std::cout << "21: " << WC(8,5) / alstilde_2 / kappa << std::endl;
1762// std::cout << "02: " << WC(8,6) / kappa_2 << std::endl;
1763// std::cout << "12: " << WC(8,7) / alstilde / kappa_2 << std::endl;
1764// std::cout << "22: " << WC(8,8) / alstilde_2 / kappa_2 << std::endl;
1765// std::cout << "C9: " << c9_df1 / aletilde << std::endl << std::endl;
1766//
1767// std::cout << "11: " << WC(9,4) / alstilde / kappa << std::endl;
1768// std::cout << "21: " << WC(9,5) / alstilde_2 / kappa << std::endl;
1769// std::cout << "02: " << WC(9,6) / kappa_2 << std::endl;
1770// std::cout << "12: " << WC(9,7) / alstilde / kappa_2 << std::endl;
1771// std::cout << "22: " << WC(9,8) / alstilde_2 / kappa_2 << std::endl;
1772// std::cout << "C10: " << c10_df1 / aletilde << std::endl;
1773//
1774// std::cout << "C9 (27): " << (c9_df1 - WC(8, int_qed(NLO_QED22)) + alstilde_2*kappa_2*27.32) / aletilde << std::endl;
1775//
1776//
1777// std::cout << "C10(-36): " << (c10_df1 - WC(9, int_qed(NLO_QED22)) +
1778// alstilde_2 * kappa_2 * (-36.09) ) / aletilde << std::endl;
1779//
1780// gslpp::complex c10_stu = ( (*(allcoeff_smm[NLO_QED11]))(7) +
1781// alstilde * (*(allcoeff_smm[NLO_QED21]))(7) +
1782// kappa / alstilde * (*(allcoeff_smm[NLO_QED02]))(7) +
1783// kappa * (*(allcoeff_smm[NLO_QED12]))(7) +
1784// alstilde * kappa * (*(allcoeff_smm[NLO_QED22]))(7) ) / lambda_t;
1785//
1786// std::cout << "C10_stu: " << c10_stu / sw2 << std::endl;
1787// std::cout << "Should be zero: " << (*(allcoeff_smm[LO]))(7) + (*(allcoeff_smm[NLO]))(7) +
1788// (*(allcoeff_smm[NNLO]))(7) + (*(allcoeff_smm[LO_QED]))(7) << std::endl;
1789//
1790// std::cout << std::endl;
1791// std::cout << "4 GF / sqrt(2) * C10 = " << 4. * GF * c10_df1 / sqrt(2.) << std::endl;
1792// std::cout << "GF^2 MW^2 / pi^2 * C10 = " << GF_2 * MW_2 * c10_stu / M_PI / M_PI << std::endl;
1793//
1794// std::cout << std::endl;
1795// std::cout << "C10_OS1: " << mySM.getMatching().C10_OS1(mt_w * mt_w / MW_2, muw) << std::endl;
1796// std::cout << "Rest: " << mySM.getMatching().Rest(mt_w * mt_w / MW_2, muw) << std::endl;
1797//
1798// const std::vector<WilsonCoefficientNew>& match_df1 = mySM.getMatching().CMDF1("CPMLQB",15);
1799// const std::vector<WilsonCoefficient>& match_stu = mySM.getMatching().CMbsmm();
1800// gslpp::complex c10_11_df1 = match_df1[0].getCoeff(QCD1, QED1)(9);
1801// gslpp::complex c10_22_df1 = match_df1[0].getCoeff(QCD2, QED2)(9);
1802// gslpp::complex c10_11_stu = (*(match_stu[0].getCoeff(orders_qed(NLO_QED11))))(7) / lambda_t;
1803// gslpp::complex c10_22_stu = (*(match_stu[0].getCoeff(orders_qed(NLO_QED22))))(7) / lambda_t;
1804//
1805// std::cout << std::endl;
1806// std::cout << "c10_11_df1: " << c10_11_df1/aletildemuw << std::endl;
1807// std::cout << "c10_11_stu: " << c10_11_stu << std::endl;
1808// std::cout << "c10_11_stu/sW2: " << c10_11_stu/sw2 << std::endl;
1809// std::cout << "c10_22_df1: " << c10_22_df1/aletildemuw/aletildemuw << std::endl;
1810// std::cout << "c10_22_stu: " << c10_22_stu << std::endl;
1811// std::cout << "c10_22_stu/sW2: " << c10_22_stu/sw2 << std::endl;
1812// std::cout << "GF matching: " << (4. * GF / sqrt(2.)) * (c10_11_df1 + c10_22_df1) << std::endl;
1813// std::cout << "GF^2 matching: " << (GF_2 * MW_2 / M_PI / M_PI) * (c10_11_stu +
1814// (mySM.Ale(muw,FULLNLO) / 4. / M_PI) * c10_22_stu) << std::endl;
1815
1816 //ATTEMPT AT A WORKING EXPANDED * EXPANDED
1817// gslpp::vector<double> vec2(2,0.);
1818// std::vector<gslpp::vector<double> > vtmp;
1819// std::vector<std::vector<gslpp::vector<double> > > vtmp12;
1820// std::vector<std::vector<gslpp::vector<double> > > vtmp910;
1821//
1822// for (int j = 0; j <= QCD2; j++)
1823// vtmp.push_back(vec2);
1824// for (int i = 0; i <= QED1; i++)
1825// vtmp12.push_back(vtmp);
1826// for (int i = 0; i <= QED2; i++)
1827// vtmp910.push_back(vtmp);
1828//
1829// Expanded<gslpp::vector<double> > wilson12(vtmp12);
1830// Expanded<gslpp::vector<double> > wilson910(vtmp910);
1831//
1832// wilson12.setVectorElement(QCD0, QED0, 0, 1.);
1833// wilson12.setVectorElement(QCD1, QED0, 0, 2.);
1834// wilson12.setVectorElement(QCD2, QED0, 0, 3.);
1835// wilson12.setVectorElement(QCD0, QED1, 0, 4.);
1836// wilson12.setVectorElement(QCD0, QED0, 1, 1.);
1837// wilson12.setVectorElement(QCD1, QED0, 1, 2.);
1838// wilson12.setVectorElement(QCD2, QED0, 1, 3.);
1839// wilson12.setVectorElement(QCD0, QED1, 1, 4.);
1840//
1841// wilson910.setVectorElement(QCD1, QED1, 0, 5.);
1842// wilson910.setVectorElement(QCD2, QED2, 0, 6.);
1843// wilson910.setVectorElement(QCD1, QED1, 1, 5.);
1844// wilson910.setVectorElement(QCD2, QED2, 1, 6.);
1845//
1846//
1847// std::cout << "00 " << wilson12.getOrd(QCD0,QED0) << std::endl;
1848// std::cout << "10 " << wilson12.getOrd(QCD1,QED0) << std::endl;
1849// std::cout << "20 " << wilson12.getOrd(QCD2,QED0) << std::endl;
1850// std::cout << "01 " << wilson12.getOrd(QCD0,QED1) << std::endl;
1851//
1852// std::cout << "11 " << wilson910.getOrd(QCD1,QED1) << std::endl;
1853// std::cout << "22 " << wilson910.getOrd(QCD2,QED2) << std::endl;
1854//
1855// std::cout << std::endl;
1856// std::cout << wilson12*wilson12 << std::endl;
1857// std::cout << std::endl;
1858// std::cout << wilson12*wilson910 << std::endl;
1859// std::cout << std::endl;
1860// std::cout << wilson910*wilson910 << std::endl;
1861
1862 //PRINT OF AUXILIARY FUNCTIONS
1863 std::cout << "mu_b = " << mu_b << std::endl;
1864 std::cout << "mu_c = " << mu_c << std::endl;
1865 std::cout << "Mb = " << Mb << std::endl;
1866 std::cout << "Mc = " << Mc << std::endl;
1867 std::cout << "Mb_pole = " << Mb_pole << std::endl;
1868 std::cout << "Mc_pole = " << Mc_pole << std::endl;
1869 std::cout << "Ms = " << Ms << std::endl;
1870 std::cout << "m_e = " << Mlep << std::endl;
1871 std::cout << "m_tau = " << Mtau << std::endl;
1872 std::cout << "alstilde = " << alstilde << std::endl;
1873 std::cout << "aletilde = " << aletilde << std::endl;
1874 std::cout << "1/ale_MZ = " << 1./mySM.alphaMz() << std::endl;
1875 std::cout << "lt/Vcb^2 = " << abslambdat_over_Vcb * abslambdat_over_Vcb << std::endl;
1876 std::cout << "BXcenu = " << BR_BXcenu << std::endl;
1877 std::cout << "C = " << C_ratio << std::endl;
1878 std::cout << "pre = " << pre << std::endl;
1879 std::cout << "lambda_1 = " << lambda_1 << std::endl;
1880 std::cout << "lambda_2 = " << lambda_2 << std::endl;
1881 std::cout << "Phi_u = " << Phi_u(FULLNLO_QED) << std::endl;
1882 std::cout << "alsMuw = " << mySM.Als(muw, FULLNNNLO, true, true) << std::endl;
1883 std::cout << "eta = " << mySM.Als(muw, FULLNNNLO, true, true) / alsmu << std::endl;
1884
1885// std::cout << std::endl;
1886// std::cout << "H_T(0.05) = " << getH("T", 0.05) << std::endl;
1887// std::cout << "H_T(0.15) = " << getH("T", 0.15) << std::endl;
1888// std::cout << "H_T(0.25) = " << getH("T", 0.25) << std::endl;
1889// std::cout << "H_L(0.05) = " << getH("L", 0.05) << std::endl;
1890// std::cout << "H_L(0.15) = " << getH("L", 0.15) << std::endl;
1891// std::cout << "H_L(0.25) = " << getH("L", 0.25) << std::endl;
1892// std::cout << "H_A(0.05) = " << getH("A", 0.05) << std::endl;
1893// std::cout << "H_A(0.15) = " << getH("A", 0.15) << std::endl;
1894// std::cout << "H_A(0.25) = " << getH("A", 0.25) << std::endl;
1895//
1896// double s = 0.15;
1897// computeMi(s);
1898//
1899// std::cout << std::endl;
1900// std::cout << "Hij_T = " << 1.0e7 * (Hij_T[LO] + Hij_T[NLO] + Hij_T[NNLO] + Hij_T[int_qed(LO_QED)] +
1901// Hij_T[int_qed(NLO_QED11)] + Hij_T[int_qed(NLO_QED21)] +
1902// Hij_T[int_qed(NLO_QED02)] + Hij_T[int_qed(NLO_QED12)] +
1903// Hij_T[int_qed(NLO_QED22)]) << std::endl;
1904// std::cout << std::endl;
1905// std::cout << "Hij_L = " << 1.0e7 * (Hij_L[LO] + Hij_L[NLO] + Hij_L[NNLO] + Hij_L[int_qed(LO_QED)] +
1906// Hij_L[int_qed(NLO_QED11)] + Hij_L[int_qed(NLO_QED21)] +
1907// Hij_L[int_qed(NLO_QED02)] + Hij_L[int_qed(NLO_QED12)] +
1908// Hij_L[int_qed(NLO_QED22)]) << std::endl;
1909//
1910// std::cout << std::endl;
1911// std::cout << "Hij_A = " << Hij_A[QCD0 + 3*QED0] + Hij_A[QCD1 + 3*QED0] + Hij_A[QCD2 + 3*QED0] +
1912// Hij_A[QCD1 + 3*QED1] + Hij_A[QCD2 + 3*QED1] + Hij_A[QCD2 + 3*QED2] +
1913// Hij_A[QCD2 + 3*QED2 + 1] << std::endl;
1914//
1915// std::cout << std::endl;
1916// std::cout << "Hij^T_00 = " << Hij_T[QCD0 + 3*QED0] << std::endl << std::endl;
1917// std::cout << "Hij^T_10 = " << Hij_T[QCD1 + 3*QED0] << std::endl << std::endl;
1918// std::cout << "Hij^T_20 = " << Hij_T[QCD2 + 3*QED0] << std::endl << std::endl;
1919// std::cout << "Hij^T_01 = " << Hij_T[QCD0 + 3*QED1] << std::endl << std::endl;
1920// std::cout << "Hij^T_11 = " << Hij_T[QCD1 + 3*QED1] << std::endl << std::endl;
1921// std::cout << "Hij^T_21 = " << Hij_T[QCD2 + 3*QED1] << std::endl << std::endl;
1922// std::cout << "Hij^T_02 = " << Hij_T[QCD0 + 3*QED2] << std::endl << std::endl;
1923// std::cout << "Hij^T_12 = " << Hij_T[QCD1 + 3*QED2] << std::endl << std::endl;
1924// std::cout << "Hij^T_22 = " << Hij_T[QCD2 + 3*QED2] << std::endl << std::endl;
1925// std::cout << "Hij^T_31 = " << Hij_T[QCD2 + 3*QED2 + 1] << std::endl << std::endl;
1926// std::cout << "Hij^T_32 = " << Hij_T[QCD2 + 3*QED2 + 2] << std::endl << std::endl;
1927// std::cout << "Hij^T_33 = " << Hij_T[QCD2 + 3*QED2 + 3] << std::endl << std::endl;
1928//
1929// std::cout << std::endl;
1930// std::cout << "Hij^L_00 = " << Hij_L[QCD0 + 3*QED0] << std::endl << std::endl;
1931// std::cout << "Hij^L_10 = " << Hij_L[QCD1 + 3*QED0] << std::endl << std::endl;
1932// std::cout << "Hij^L_20 = " << Hij_L[QCD2 + 3*QED0] << std::endl << std::endl;
1933// std::cout << "Hij^L_01 = " << Hij_L[QCD0 + 3*QED1] << std::endl << std::endl;
1934// std::cout << "Hij^L_11 = " << Hij_L[QCD1 + 3*QED1] << std::endl << std::endl;
1935// std::cout << "Hij^L_21 = " << Hij_L[QCD2 + 3*QED1] << std::endl << std::endl;
1936// std::cout << "Hij^L_02 = " << Hij_L[QCD0 + 3*QED2] << std::endl << std::endl;
1937// std::cout << "Hij^L_12 = " << Hij_L[QCD1 + 3*QED2] << std::endl << std::endl;
1938// std::cout << "Hij^L_22 = " << Hij_L[QCD2 + 3*QED2] << std::endl << std::endl;
1939// std::cout << "Hij^L_31 = " << Hij_L[QCD2 + 3*QED2 + 1] << std::endl << std::endl;
1940// std::cout << "Hij^L_32 = " << Hij_L[QCD2 + 3*QED2 + 2] << std::endl << std::endl;
1941// std::cout << "Hij^L_33 = " << Hij_L[QCD2 + 3*QED2 + 3] << std::endl << std::endl;
1942//
1943// std::cout << std::endl;
1944// std::cout << "Hij^A_00 = " << Hij_A[QCD0 + 3*QED0] << std::endl << std::endl;
1945// std::cout << "Hij^A_10 = " << Hij_A[QCD1 + 3*QED0] << std::endl << std::endl;
1946// std::cout << "Hij^A_20 = " << Hij_A[QCD2 + 3*QED0] << std::endl << std::endl;
1947// std::cout << "Hij^A_01 = " << Hij_A[QCD0 + 3*QED1] << std::endl << std::endl;
1948// std::cout << "Hij^A_11 = " << Hij_A[QCD1 + 3*QED1] << std::endl << std::endl;
1949// std::cout << "Hij^A_21 = " << Hij_A[QCD2 + 3*QED1] << std::endl << std::endl;
1950// std::cout << "Hij^A_02 = " << Hij_A[QCD0 + 3*QED2] << std::endl << std::endl;
1951// std::cout << "Hij^A_12 = " << Hij_A[QCD1 + 3*QED2] << std::endl << std::endl;
1952// std::cout << "Hij^A_22 = " << Hij_A[QCD2 + 3*QED2] << std::endl << std::endl;
1953// std::cout << "Hij^A_31 = " << Hij_A[QCD2 + 3*QED2 + 1] << std::endl;
1954// std::cout << "Hij^A_32 = " << Hij_A[QCD2 + 3*QED2 + 2] << std::endl;
1955// std::cout << "Hij^A_33 = " << Hij_A[QCD2 + 3*QED2 + 3] << std::endl;
1956//
1957// std::cout << std::endl;
1958// std::cout << "H11^T = (";
1959// std::cout << 1.0e7 * Hij_T[QCD0 + 3*QED0](0,0) << ", " << 1.0e7 * Hij_T[QCD1 + 3*QED0](0,0) << ", ";
1960// std::cout << 1.0e7 * Hij_T[QCD2 + 3*QED0](0,0) << ", " << 1.0e7 * Hij_T[QCD0 + 3*QED1](0,0) << ", ";
1961// std::cout << 1.0e7 * Hij_T[QCD1 + 3*QED1](0,0) << ", " << 1.0e7 * Hij_T[QCD2 + 3*QED1](0,0) << ", ";
1962// std::cout << 1.0e7 * Hij_T[QCD0 + 3*QED2](0,0) << ", " << 1.0e7 * Hij_T[QCD1 + 3*QED2](0,0) << ", ";
1963// std::cout << 1.0e7 * Hij_T[QCD2 + 3*QED2](0,0) << ")" << std::endl;
1964//
1965// std::cout << std::endl;
1966// std::cout << "H12^T = (";
1967// std::cout << 1.0e7 * Hij_T[QCD0 + 3*QED0](0,1) << ", " << 1.0e7 * Hij_T[QCD1 + 3*QED0](0,1) << ", ";
1968// std::cout << 1.0e7 * Hij_T[QCD2 + 3*QED0](0,1) << ", " << 1.0e7 * Hij_T[QCD0 + 3*QED1](0,1) << ", ";
1969// std::cout << 1.0e7 * Hij_T[QCD1 + 3*QED1](0,1) << ", " << 1.0e7 * Hij_T[QCD2 + 3*QED1](0,1) << ", ";
1970// std::cout << 1.0e7 * Hij_T[QCD0 + 3*QED2](0,1) << ", " << 1.0e7 * Hij_T[QCD1 + 3*QED2](0,1) << ", ";
1971// std::cout << 1.0e7 * Hij_T[QCD2 + 3*QED2](0,1) << ")" << std::endl;
1972//
1973// std::cout << std::endl;
1974// std::cout << "H22^T = (";
1975// std::cout << 1.0e7 * Hij_T[QCD0 + 3*QED0](1,1) << ", " << 1.0e7 * Hij_T[QCD1 + 3*QED0](1,1) << ", ";
1976// std::cout << 1.0e7 * Hij_T[QCD2 + 3*QED0](1,1) << ", " << 1.0e7 * Hij_T[QCD0 + 3*QED1](1,1) << ", ";
1977// std::cout << 1.0e7 * Hij_T[QCD1 + 3*QED1](1,1) << ", " << 1.0e7 * Hij_T[QCD2 + 3*QED1](1,1) << ", ";
1978// std::cout << 1.0e7 * Hij_T[QCD0 + 3*QED2](1,1) << ", " << 1.0e7 * Hij_T[QCD1 + 3*QED2](1,1) << ", ";
1979// std::cout << 1.0e7 * Hij_T[QCD2 + 3*QED2](1,1) << ")" << std::endl;
1980//
1981// std::cout << std::endl;
1982// std::cout << "H99^T = (";
1983// std::cout << 1.0e7 * Hij_T[QCD0 + 3*QED0](8,8) << ", " << 1.0e7 * Hij_T[QCD1 + 3*QED0](8,8) << ", ";
1984// std::cout << 1.0e7 * Hij_T[QCD2 + 3*QED0](8,8) << ", " << 1.0e7 * Hij_T[QCD0 + 3*QED1](8,8) << ", ";
1985// std::cout << 1.0e7 * Hij_T[QCD1 + 3*QED1](8,8) << ", " << 1.0e7 * Hij_T[QCD2 + 3*QED1](8,8) << ", ";
1986// std::cout << 1.0e7 * Hij_T[QCD0 + 3*QED2](8,8) << ", " << 1.0e7 * Hij_T[QCD1 + 3*QED2](8,8) << ", ";
1987// std::cout << 1.0e7 * Hij_T[QCD2 + 3*QED2](8,8) << ")" << std::endl;
1988//
1989// std::cout << std::endl;
1990// std::cout << "H1010^T = (";
1991// std::cout << 1.0e7 * Hij_T[QCD0 + 3*QED0](9,9) << ", " << 1.0e7 * Hij_T[QCD1 + 3*QED0](9,9) << ", ";
1992// std::cout << 1.0e7 * Hij_T[QCD2 + 3*QED0](9,9) << ", " << 1.0e7 * Hij_T[QCD0 + 3*QED1](9,9) << ", ";
1993// std::cout << 1.0e7 * Hij_T[QCD1 + 3*QED1](9,9) << ", " << 1.0e7 * Hij_T[QCD2 + 3*QED1](9,9) << ", ";
1994// std::cout << 1.0e7 * Hij_T[QCD0 + 3*QED2](9,9) << ", " << 1.0e7 * Hij_T[QCD1 + 3*QED2](9,9) << ", ";
1995// std::cout << 1.0e7 * Hij_T[QCD2 + 3*QED2](9,9) << ")" << std::endl;
1996//
1997//
1998// std::cout << std::endl;
1999// std::cout << "H11^L = (";
2000// std::cout << 1.0e7 * Hij_L[QCD0 + 3*QED0](0,0) << ", " << 1.0e7 * Hij_L[QCD1 + 3*QED0](0,0) << ", ";
2001// std::cout << 1.0e7 * Hij_L[QCD2 + 3*QED0](0,0) << ", " << 1.0e7 * Hij_L[QCD0 + 3*QED1](0,0) << ", ";
2002// std::cout << 1.0e7 * Hij_L[QCD1 + 3*QED1](0,0) << ", " << 1.0e7 * Hij_L[QCD2 + 3*QED1](0,0) << ", ";
2003// std::cout << 1.0e7 * Hij_L[QCD0 + 3*QED2](0,0) << ", " << 1.0e7 * Hij_L[QCD1 + 3*QED2](0,0) << ", ";
2004// std::cout << 1.0e7 * Hij_L[QCD2 + 3*QED2](0,0) << ")" << std::endl;
2005//
2006// std::cout << std::endl;
2007// std::cout << "H12^L = (";
2008// std::cout << 1.0e7 * Hij_L[QCD0 + 3*QED0](0,1) << ", " << 1.0e7 * Hij_L[QCD1 + 3*QED0](0,1) << ", ";
2009// std::cout << 1.0e7 * Hij_L[QCD2 + 3*QED0](0,1) << ", " << 1.0e7 * Hij_L[QCD0 + 3*QED1](0,1) << ", ";
2010// std::cout << 1.0e7 * Hij_L[QCD1 + 3*QED1](0,1) << ", " << 1.0e7 * Hij_L[QCD2 + 3*QED1](0,1) << ", ";
2011// std::cout << 1.0e7 * Hij_L[QCD0 + 3*QED2](0,1) << ", " << 1.0e7 * Hij_L[QCD1 + 3*QED2](0,1) << ", ";
2012// std::cout << 1.0e7 * Hij_L[QCD2 + 3*QED2](0,1) << ")" << std::endl;
2013//
2014// std::cout << std::endl;
2015// std::cout << "H22^L = (";
2016// std::cout << 1.0e7 * Hij_L[QCD0 + 3*QED0](1,1) << ", " << 1.0e7 * Hij_L[QCD1 + 3*QED0](1,1) << ", ";
2017// std::cout << 1.0e7 * Hij_L[QCD2 + 3*QED0](1,1) << ", " << 1.0e7 * Hij_L[QCD0 + 3*QED1](1,1) << ", ";
2018// std::cout << 1.0e7 * Hij_L[QCD1 + 3*QED1](1,1) << ", " << 1.0e7 * Hij_L[QCD2 + 3*QED1](1,1) << ", ";
2019// std::cout << 1.0e7 * Hij_L[QCD0 + 3*QED2](1,1) << ", " << 1.0e7 * Hij_L[QCD1 + 3*QED2](1,1) << ", ";
2020// std::cout << 1.0e7 * Hij_L[QCD2 + 3*QED2](1,1) << ")" << std::endl;
2021//
2022// std::cout << std::endl;
2023// std::cout << "H99^L = (";
2024// std::cout << 1.0e7 * Hij_L[QCD0 + 3*QED0](8,8) << ", " << 1.0e7 * Hij_L[QCD1 + 3*QED0](8,8) << ", ";
2025// std::cout << 1.0e7 * Hij_L[QCD2 + 3*QED0](8,8) << ", " << 1.0e7 * Hij_L[QCD0 + 3*QED1](8,8) << ", ";
2026// std::cout << 1.0e7 * Hij_L[QCD1 + 3*QED1](8,8) << ", " << 1.0e7 * Hij_L[QCD2 + 3*QED1](8,8) << ", ";
2027// std::cout << 1.0e7 * Hij_L[QCD0 + 3*QED2](8,8) << ", " << 1.0e7 * Hij_L[QCD1 + 3*QED2](8,8) << ", ";
2028// std::cout << 1.0e7 * Hij_L[QCD2 + 3*QED2](8,8) << ")" << std::endl;
2029//
2030// std::cout << std::endl;
2031// std::cout << "H1010^L = (";
2032// std::cout << 1.0e7 * Hij_L[QCD0 + 3*QED0](9,9) << ", " << 1.0e7 * Hij_L[QCD1 + 3*QED0](9,9) << ", ";
2033// std::cout << 1.0e7 * Hij_L[QCD2 + 3*QED0](9,9) << ", " << 1.0e7 * Hij_L[QCD0 + 3*QED1](9,9) << ", ";
2034// std::cout << 1.0e7 * Hij_L[QCD1 + 3*QED1](9,9) << ", " << 1.0e7 * Hij_L[QCD2 + 3*QED1](9,9) << ", ";
2035// std::cout << 1.0e7 * Hij_L[QCD0 + 3*QED2](9,9) << ", " << 1.0e7 * Hij_L[QCD1 + 3*QED2](9,9) << ", ";
2036// std::cout << 1.0e7 * Hij_L[QCD2 + 3*QED2](9,9) << ")" << std::endl;
2037//
2038// std::cout << std::endl;
2039// std::cout << "H110^A = (";
2040// std::cout << 1.0e5 * Hij_A[QCD0 + 3*QED0](0,9) << ", " << 1.0e5 * Hij_A[QCD1 + 3*QED0](0,9) << ", ";
2041// std::cout << 1.0e5 * Hij_A[QCD2 + 3*QED0](0,9) << ", " << 1.0e5 * Hij_A[QCD1 + 3*QED1](0,9) << ", ";
2042// std::cout << 1.0e5 * Hij_A[QCD2 + 3*QED1](0,9) << ", " << 1.0e5 * Hij_A[QCD2 + 3*QED2 + 1](0,9) << ")";
2043// std::cout << std::endl;
2044//
2045// std::cout << std::endl;
2046// std::cout << "H310^A = (";
2047// std::cout << 1.0e5 * Hij_A[QCD0 + 3*QED0](2,9) << ", " << 1.0e5 * Hij_A[QCD1 + 3*QED0](2,9) << ", ";
2048// std::cout << 1.0e5 * Hij_A[QCD2 + 3*QED0](2,9) << ", " << 1.0e5 * Hij_A[QCD1 + 3*QED1](2,9) << ", ";
2049// std::cout << 1.0e5 * Hij_A[QCD2 + 3*QED1](2,9) << ", " << 1.0e5 * Hij_A[QCD2 + 3*QED2 + 1](2,9) << ")";
2050// std::cout << std::endl;
2051//
2052// std::cout << std::endl;
2053// std::cout << "KS_cc = " << KS_cc(s) << std::endl;
2054//
2055// std::cout << std::endl;
2056// for (unsigned int i = 0; i < 10; i++)
2057// {
2058// std::cout << "M_" << i+1 << "^7 = " << M_7[LO](i) + M_7[NLO](i) + M_7[NNLO](i) + M_7[int_qed(LO_QED)](i) +
2059// M_7[int_qed(NLO_QED11)](i) + M_7[int_qed(NLO_QED21)](i) << std::endl;
2060// std::cout << "M_" << i+1 << "^9 = " << M_9[LO](i) + M_9[NLO](i) + M_9[NNLO](i) + M_9[int_qed(LO_QED)](i) +
2061// M_9[int_qed(NLO_QED11)](i) + M_9[int_qed(NLO_QED21)](i) << std::endl;
2062// }
2063//
2064// for (unsigned int i = 10; i < 14; i++)
2065// {
2066// std::cout << "M_" << i-7 << "Q^7 = " << M_7[LO](i) + M_7[NLO](i) + M_7[NNLO](i) + M_7[int_qed(LO_QED)](i) +
2067// M_7[int_qed(NLO_QED11)](i) + M_7[int_qed(NLO_QED21)](i) << std::endl;
2068// std::cout << "M_" << i-7 << "Q^9 = " << M_9[LO](i) + M_9[NLO](i) + M_9[NNLO](i) + M_9[int_qed(LO_QED)](i) +
2069// M_9[int_qed(NLO_QED11)](i) + M_9[int_qed(NLO_QED21)](i) << std::endl;
2070// }
2071//
2072// std::cout << "M_b^7 = " << M_7[LO](14) + M_7[NLO](14) + M_7[NNLO](14) + M_7[int_qed(LO_QED)](14) +
2073// M_7[int_qed(NLO_QED11)](14) + M_7[int_qed(NLO_QED21)](14) << std::endl;
2074// std::cout << "M_b^9 = " << M_9[LO](14) + M_9[NLO](14) + M_9[NNLO](14) + M_9[int_qed(LO_QED)](14) +
2075// M_9[int_qed(NLO_QED11)](14) + M_9[int_qed(NLO_QED21)](14) << std::endl;
2076//
2077// std::cout << std::endl;
2078// std::cout << "S_77^T = " << S77_T(s, LO) + S77_T(s, NLO) << std::endl;
2079// std::cout << "S_99^T = " << S99_T(s, LO) + S99_T(s, NLO) + S99_T(s, NNLO) << std::endl;
2080// std::cout << "S_79^T = " << S79_T(s, LO) + S79_T(s, NLO) << std::endl;
2081// std::cout << "S_1010^T = " << S1010_T(s, LO) + S1010_T(s, NLO) + S1010_T(s, NNLO) << std::endl;
2082//
2083// std::cout << std::endl;
2084// std::cout << "S_77^L = " << S77_L(s, LO) + S77_L(s, NLO) << std::endl;
2085// std::cout << "S_99^L = " << S99_L(s, LO) + S99_L(s, NLO) + S99_L(s, NNLO) << std::endl;
2086// std::cout << "S_79^L = " << S79_L(s, LO) + S79_L(s, NLO) << std::endl;
2087// std::cout << "S_1010^L = " << S1010_L(s, LO) + S1010_L(s, NLO) + S1010_L(s, NNLO) << std::endl;
2088//
2089// std::cout << std::endl;
2090// std::cout << "S_710^A = " << S710_A(s, LO) + S710_A(s, NLO) << std::endl;
2091// std::cout << "S_910^A = " << S910_A(s, LO) + S910_A(s, NLO) + S910_A(s, NNLO) << std::endl;
2092//
2093// std::cout << std::endl;
2094// std::cout << "f_1 = " << f_Huber(s, -32./27., 4./3., 0., 0., -16./27.) << std::endl;
2095// std::cout << "f_2 = " << f_Huber(s, -8./9., 1., 0., 0., -4./9.) << std::endl;
2096// std::cout << "f_3 = " << f_Huber(s, -16./9., 6., -7./2., 2./9., 2./27.) << std::endl;
2097// std::cout << "f_4 = " << f_Huber(s, 32./27., 0., -2./3., 8./27., 8./81.) << std::endl;
2098// std::cout << "f_5 = " << f_Huber(s, -112./9., 60., -38., 32./9., -136./27.) << std::endl;
2099// std::cout << "f_6 = " << f_Huber(s, 512./27., 0., -32./3., 128./27., 320./81.) << std::endl;
2100// std::cout << "f_9 = " << f9pen_Huber(s) << std::endl;
2101// std::cout << "f_3Q = " << f_Huber(s, -272./27., 4., 7./6., -74./27., 358./81.) << std::endl;
2102// std::cout << "f_4Q = " << f_Huber(s, -32./81., 0., 2./9., -8./81., -8./243.) << std::endl;
2103// std::cout << "f_5Q = " << f_Huber(s, -2768./27., 40., 38./3., -752./27., 1144./81.) << std::endl;
2104// std::cout << "f_6Q = " << f_Huber(s, -512./81., 0., 32./9., -128./81., -320./243.) << std::endl;
2105// std::cout << "f_b = " << f_Huber(s, 16./9., 0., -2., 0., 26./27.) << std::endl;
2106//
2107// std::cout << std::endl;
2108// std::cout << "F_1^7 = " << F17(s) << std::endl;
2109// std::cout << "F_2^7 = " << F27(s) << std::endl;
2110// std::cout << "F_8^7 = " << F87(s) << std::endl;
2111// std::cout << "F_1^9 = " << F19(s) << std::endl;
2112// std::cout << "F_2^9 = " << F29(s) << std::endl;
2113// std::cout << "F_8^9 = " << F89(s) << std::endl;
2114//
2115// double umsh = 1. - s;
2116// std::cout << std::endl;
2117// std::cout << "sigma_77^T = " << 8.*umsh*umsh/s << std::endl;
2118// std::cout << "sigma_99^T = " << 2.*s*umsh*umsh << std::endl;
2119// std::cout << "sigma_79^T = " << 8.*umsh*umsh << std::endl;
2120// std::cout << "sigma_77^L = " << 4.*umsh*umsh << std::endl;
2121// std::cout << "sigma_99^L = " << umsh*umsh << std::endl;
2122// std::cout << "sigma_79^L = " << 4.*umsh*umsh << std::endl;
2123//
2124// std::cout << std::endl;
2125// std::cout << "w1_77^T = " << omega77_T(s) << std::endl;
2126// std::cout << "w1_99^T = " << omega99_T(s) << std::endl;
2127// std::cout << "w1_79^T = " << omega79_T(s) << std::endl;
2128// std::cout << "w1_77^L = " << omega77_L(s) << std::endl;
2129// std::cout << "w1_99^L = " << omega99_L(s) << std::endl;
2130// std::cout << "w1_79^L = " << omega79_L(s) << std::endl;
2131//
2132// std::cout << std::endl;
2133// std::cout << "chi1_77^T = " << 4.*umsh*(5.*s + 3.)/3./s << std::endl;
2134// std::cout << "chi1_99^T = " << -s*umsh*(3.*s + 5.)/3. << std::endl;
2135// std::cout << "chi1_79^T = " << 4.*umsh*umsh << std::endl;
2136// std::cout << "chi1_77^L = " << -2.*umsh*(3.*s + 13.)/3. << std::endl;
2137// std::cout << "chi1_99^L = " << umsh*(13.*s + 3.)/6. << std::endl;
2138// std::cout << "chi1_79^L = " << 2.*umsh*umsh << std::endl;
2139//
2140// std::cout << std::endl;
2141// std::cout << "chi2_77^T = " << 4.*(3.*s*s + 2.*s - 9.)/s << std::endl;
2142// std::cout << "chi2_99^T = " << s*(15.*s*s - 14.*s - 5.) << std::endl;
2143// std::cout << "chi2_79^T = " << 4.*(9.*s*s - 6.*s - 7.) << std::endl;
2144// std::cout << "chi2_77^L = " << 2.*(15.*s*s - 6.*s - 13.) << std::endl;
2145// std::cout << "chi2_99^L = " << (-17.*s*s + 10.*s + 3.)/2. << std::endl;
2146// std::cout << "chi2_79^L = " << 2.*(3.*s*s - 6.*s - 1.) << std::endl;
2147//
2148// std::cout << std::endl;
2149// for (unsigned int i = 0; i < 2; i++)
2150// {
2151// for (unsigned int j = 0; j < 15; j++)
2152// std::cout << "c^T_" << i+1 << "," << j+1 << " = " << cij_T(i, j, s, 11) + cij_T(i, j, s, 22) + cij_T(i, j, s, 32) << std::endl;
2153// std::cout << std::endl;
2154// }
2155//
2156// for (unsigned int i = 0; i < 2; i++)
2157// {
2158// for (unsigned int j = 0; j < 15; j++)
2159// std::cout << "c^L_" << i+1 << "," << j+1 << " = " << cij_L(i, j, s, 11) + cij_L(i, j, s, 22) + cij_L(i, j, s, 32) << std::endl;
2160// std::cout << std::endl;
2161// }
2162//
2163// std::cout << std::endl;
2164// std::cout << "c^A_1,10 = " << cij_A(0, 9, s) << std::endl;
2165// std::cout << "c^A_2,10 = " << cij_A(1, 9, s) << std::endl;
2166//
2167// std::cout << std::endl;
2168// std::cout << "e^T_1,1 = " << eij_T(0, 0, s) << std::endl;
2169// std::cout << "e^T_1,2 = " << eij_T(0, 1, s) << std::endl;
2170// std::cout << "e^T_1,7 = " << eij_T(0, 6, s) << std::endl;
2171// std::cout << "e^T_1,9 = " << eij_T(0, 8, s) << std::endl;
2172// std::cout << "e^T_2,2 = " << eij_T(1, 1, s) << std::endl;
2173// std::cout << "e^T_2,7 = " << eij_T(1, 6, s) << std::endl;
2174// std::cout << "e^T_2,9 = " << eij_T(1, 8, s) << std::endl;
2175// std::cout << "e^T_7,7 = " << eij_T(6, 6, s) << std::endl;
2176// std::cout << "e^T_7,9 = " << eij_T(6, 8, s) << std::endl;
2177// std::cout << "e^T_9,9 = " << eij_T(8, 8, s) << std::endl;
2178// std::cout << "e^T_10,10 = " << eij_T(9, 9, s) << std::endl;
2179//
2180// std::cout << std::endl;
2181// std::cout << "e^L_1,1 = " << eij_L(0, 0, s) << std::endl;
2182// std::cout << "e^L_1,2 = " << eij_L(0, 1, s) << std::endl;
2183// std::cout << "e^L_1,7 = " << eij_L(0, 6, s) << std::endl;
2184// std::cout << "e^L_1,9 = " << eij_L(0, 8, s) << std::endl;
2185// std::cout << "e^L_2,2 = " << eij_L(1, 1, s) << std::endl;
2186// std::cout << "e^L_2,7 = " << eij_L(1, 6, s) << std::endl;
2187// std::cout << "e^L_2,9 = " << eij_L(1, 8, s) << std::endl;
2188// std::cout << "e^L_7,7 = " << eij_L(6, 6, s) << std::endl;
2189// std::cout << "e^L_7,9 = " << eij_L(6, 8, s) << std::endl;
2190// std::cout << "e^L_9,9 = " << eij_L(8, 8, s) << std::endl;
2191// std::cout << "e^L_10,10 = " << eij_L(9, 9, s) << std::endl;
2192//
2193// std::cout << std::endl;
2194// std::cout << "e^A_1,10 = " << eij_A(0, 9, s) << std::endl;
2195// std::cout << "e^A_2,10 = " << eij_A(1, 9, s) << std::endl;
2196// std::cout << "e^A_7,10 = " << eij_A(6, 9, s) << std::endl;
2197// std::cout << "e^A_9,10 = " << eij_A(8, 9, s) << std::endl;
2198//
2199// std::cout << std::endl;
2200// for (unsigned int i = 0; i < 15; i++)
2201// {
2202// std::cout << "C_" << i+1 << " = (";
2203// for(unsigned int qed_ord = QED0; qed_ord <= QED2; qed_ord++)
2204// for(unsigned int qcd_ord = QCD0; qcd_ord <= QCD2; qcd_ord++)
2205// std::cout << WC(i, qcd_ord + 3*qed_ord) / pow(alstilde, qcd_ord) / pow(kappa, qed_ord) << ", ";
2206// std::cout << ")" << std::endl;
2207// }
2208//
2209// std::cout << std::endl;
2210// for (unsigned int i = 0; i < 15; i++)
2211// {
2212// std::cout << "C_" << i+1 << " = " << WC(i, QCD0 + 3*QED0) + WC(i, QCD1 + 3*QED0) + WC(i, QCD2 + 3*QED0) +
2213// WC(i, QCD0 + 3*QED1) + WC(i, QCD1 + 3*QED1) + WC(i, QCD2 + 3*QED1) +
2214// WC(i, QCD0 + 3*QED2) + WC(i, QCD1 + 3*QED2) + WC(i, QCD2 + 3*QED2) << std::endl;
2215// }
2216//
2217// double akmu;
2218// std::cout << std::endl;
2219// for (unsigned int qedi = QED0; qedi <= QED2; qedi++)
2220// for (unsigned int qcdi = QCD0; qcdi <= QCD2; qcdi++)
2221// {
2222// akmu = pow(alstilde, qcdi) * pow(kappa, qedi);
2223// std::cout << "WC(" << qcdi << "," << qedi << ") = (";
2224// std::cout << WC(0,qcdi+3*qedi)/akmu << ", " << WC(1,qcdi+3*qedi)/akmu << ", " << WC(2,qcdi+3*qedi)/akmu << ", ";
2225// std::cout << WC(3,qcdi+3*qedi)/akmu << ", " << WC(4,qcdi+3*qedi)/akmu << ", " << WC(5,qcdi+3*qedi)/akmu << ", ";
2226// std::cout << WC(6,qcdi+3*qedi)/akmu << ", " << WC(7,qcdi+3*qedi)/akmu << ", " << WC(8,qcdi+3*qedi)/akmu << ", ";
2227// std::cout << WC(9,qcdi+3*qedi)/akmu << ", " << WC(10,qcdi+3*qedi)/akmu << ", " << WC(11,qcdi+3*qedi)/akmu << ", ";
2228// std::cout << WC(12,qcdi+3*qedi)/akmu << ", " << WC(13,qcdi+3*qedi)/akmu << ", " << WC(14,qcdi+3*qedi)/akmu << ")";
2229// std::cout << std::endl;
2230// }
2231//
2232// double mymt = mySM.Mrun(muw, mySM.getQuarks(QCD::TOP).getMass_scale(),
2233// mySM.getQuarks(QCD::TOP).getMass(), FULLNNLO);
2234// std::cout << std::endl;
2235// std::cout << "Mw = " << MW << std::endl;
2236// std::cout << "Muw = " << muw << std::endl;
2237// std::cout << "xt = " << xt << std::endl;
2238// std::cout << "mt = " << mymt << " = " << sqrt(xt * MW * MW) << std::endl;
2239// std::cout << "sW2 = " << mySM.sW2() << std::endl;
2240//
2241// double alsMuw = mySM.Als(muw, FULLNNNLO, true) / 4. / M_PI;
2242// double aleMuw = mySM.Ale(muw, FULLNLO) / 4. / M_PI;
2243// std::cout << std::endl;
2244// std::cout << "alstilde(Muw) = " << alsMuw << std::endl;
2245// std::cout << "aletilde(Muw) = " << aleMuw << std::endl;
2246//
2247// std::cout << std::endl;
2248// std::cout << "C7t_3L_at_mt(xt) = " << mySM.getMatching().C7t_3L_at_mt(xt) << std::endl;
2249// std::cout << "C7t_3L_func(xt,Muw) = " << mySM.getMatching().C7t_3L_func(xt, muw) << std::endl;
2250// std::cout << "C7c_3L_at_mW(xt) = " << mySM.getMatching().C7c_3L_at_mW(xt) << std::endl;
2251// std::cout << "C8t_3L_at_mt(xt) = " << mySM.getMatching().C8t_3L_at_mt(xt) << std::endl;
2252// std::cout << "C8t_3L_func(xt,Muw) = " << mySM.getMatching().C8t_3L_func(xt, muw) << std::endl;
2253// std::cout << "C8c_3L_at_mW(xt) = " << mySM.getMatching().C8c_3L_at_mW(xt) << std::endl;
2254//
2255// std::vector<WilsonCoefficientNew>& mymc = mySM.getMatching().CMDF1("CPMLQB", 15);
2256// std::cout << std::endl;
2257// for (unsigned int qedi = QED0; qedi <= QED2; qedi++)
2258// for (unsigned int qcdi = QCD0; qcdi <= QCD2; qcdi++)
2259// {
2260// akmu = pow(alsMuw, qcdi) * pow(aleMuw / alsMuw, qedi);
2261// std::cout << "mc(" << qcdi << "," << qedi << ") = ";
2262// std::cout << mymc[0].getCoeff((qcd_orders) qcdi, (qed_orders) qedi) / akmu << std::endl;
2263// }
2264//
2265// std::cout << myHeff.getEvol().DF1Evol(mu_b, muw, NDR).getOrd(QCD0,QED1) * alsMuw / aleMuw << std::endl;
2266//
2267// std::cout << std::endl;
2268// std::cout << "PhiTL_brems(0.05) = " << PhiTL_brems(0.05) << std::endl;
2269// std::cout << "PhiTL_brems(0.15) = " << PhiTL_brems(0.15) << std::endl;
2270// std::cout << "PhiTL_brems(0.25) = " << PhiTL_brems(0.25) << std::endl;
2271// std::cout << "PhiA_brems(0.05) = " << PhiA_brems(0.05) << std::endl;
2272// std::cout << "PhiA_brems(0.15) = " << PhiA_brems(0.15) << std::endl;
2273// std::cout << "PhiA_brems(0.25) = " << PhiA_brems(0.25) << std::endl;
2274}
2275
2276
2277/********************************************************************************************
2278 * Finite bremsstrahlung corrections in the notation of Asatryan:2002iy and Asatrian:2003yk *
2279 ********************************************************************************************/
2280
2281double BXqll::PhiTL_brems(double sh)
2282{
2283 gslpp::complex c78, c89, c17, c27, c18, c28, c19, c29;
2284 double c88, c11, c12, c22;
2285 double Brems_a;
2286 double Brems_b;
2287 double ctau1 = (3. * 3. - 1.) / 8. / 3. / 3. / 3.;
2288 double ctau2 = - (3. * 3. - 1.) / 4. / 3. / 3.;
2289
2290 c78 = CF * WC(6, LO) * WC(7, LO).conjugate();
2291 c89 = CF * WC(7, LO) * C9mod(sh).conjugate();
2292 c88 = CF * WC(7,LO).abs2();
2293
2294 c11 = ctau1 * WC(0, LO).abs2();
2295 c12 = ctau2 * 2. * (WC(0, LO) * WC(1, LO).conjugate()).real();
2296 c22 = CF * WC(1, LO).abs2();
2297
2298 c17 = ctau2 * WC(0, LO) * WC(6, LO).conjugate();
2299 c18 = ctau2 * WC(0, LO) * WC(7, LO).conjugate();
2300 c19 = ctau2 * WC(0, LO) * C9mod(sh).conjugate();
2301 c27 = CF * WC(1, LO) * WC(6, LO).conjugate();
2302 c28 = CF * WC(1, LO) * WC(7, LO).conjugate();
2303 c29 = CF * WC(1, LO) * C9mod(sh).conjugate();
2304
2305 Brems_a = 2. * (c78 * tau78(sh) + c89 * tau89(sh)).real() + c88 * tau88(sh);
2306
2307 Brems_b = (c11 + c12 + c22) * tau22fit(sh) + 2. * (c17 + c27).real() * tau27fit_Re(sh) -
2308 2. * (c17 + c27).imag() * tau27fit_Im(sh) + 2. * (c18 + c28).real() * tau28fit_Re(sh) -
2309 2. * (c18 + c28).imag() * tau28fit_Im(sh) + 2. * (c19 + c29).real() * tau29fit_Re(sh) -
2310 2. * (c19 + c29).imag() * tau29fit_Im(sh);
2311
2312 return (aletilde * aletilde * alstilde * (Brems_a + Brems_b));
2313}
2314
2315double BXqll::PhiA_brems(double sh)
2316{
2317 gslpp::complex C10mod = (WC(9, QCD1 + 3*QED1) / aletilde).conjugate();
2318 // 2 als / pi = 8 alstilde, and H_A(s) = dAFB/ds * 4 / 3
2319 double factor = aletilde * aletilde * alstilde * 8. * (1. - sh) * (1. - sh) * 4. / 3.;
2320
2321 return (factor * ((WC(7, LO) * C10mod).real() * t810(sh) +
2322 ((WC(1, LO) - WC(0, LO) / 6.) * C10mod).real() * t210fit(sh)) / 3.);
2323}
2324
2325gslpp::complex BXqll::C9mod(double sh)
2326{
2327 gslpp::complex T9 = 4./3. * WC(0, LO) + WC(1, LO) + 6. * WC(2, LO) + 60. * WC(4, LO);
2328 gslpp::complex U9 = -7./2. * WC(2, LO) - 2./3. * WC(3, LO) - 38. * WC(4, LO) - 32./3. * WC(5, LO);
2329 gslpp::complex W9 = -1./2. * WC(2, LO) - 2./3. * WC(3, LO) - 8. * WC(4, LO) - 32./3. * WC(5, LO);
2330 gslpp::complex A9 = -32./27. * WC(0, LO) - 8./9. * WC(1, LO) - 16./9. * WC(2, LO) +
2331 32./27. * WC(3, LO) - 112./9. * WC(4, LO) + 512./27. * WC(5, LO);
2332 A9 *= log(1. / muh);
2333 A9 += (WC(8, QCD0 + 3*QED1) + WC(8, QCD1 + 3*QED1)) / aletilde;
2334 A9 += 4./3. * WC(2, LO) + 64./9. * WC(4, LO) + 64./27. * WC(5, LO);
2335
2336 return (A9 + T9 * h_z(z, sh) + U9 * h_z(1., sh) + W9 * h_z(0., sh));
2337}
2338
2339gslpp::complex BXqll::h_z(double zed, double sh)
2340{
2341 gslpp::complex i = gslpp::complex::i();
2342 gslpp::complex h_z;
2343
2344 if (zed == 0.)
2345 {
2346 h_z = 8./27.-4./9.*(log(sh)-i*M_PI);
2347 }
2348 else
2349 {
2350 h_z = -2./9.*(2.+4./sh*zed)*sqrt(std::abs(4.*zed-sh)/sh);
2351 if(sh > 4.*zed)
2352 h_z *= log((sqrt(sh)+sqrt(sh-4.*zed))/(sqrt(sh)-sqrt(sh-4.*zed)))-i*M_PI;
2353 else
2354 h_z *= 2.*atan(sqrt(sh/(4.*zed-sh)));
2355
2356 h_z += -4./9.*log(zed)+8./27.+16./9.*zed/sh;
2357 }
2358
2359 return(h_z);
2360}
2361
2362double BXqll::tau78(double sh)
2363{
2364 gslpp::complex i = gslpp::complex::i();
2365 double dmsh = 2.-sh;
2366 double qmsh = 4.-sh;
2367
2368 return (8./9./sh*(25.-2.*M_PI*M_PI-27.*sh+3.*sh*sh-sh*sh*sh+12.*(sh+sh*sh)*log(sh)+
2369 6.*(M_PI/2.-atan((2.-4.*sh+sh*sh)/dmsh/sqrt(sh)/sqrt(qmsh)))*
2370 (M_PI/2.-atan((2.-4.*sh+sh*sh)/dmsh/sqrt(sh)/sqrt(qmsh)))-
2371 24.*(dilog((sh-i*sqrt(sh)*sqrt(qmsh))/2.)).real()-12.*((1.-sh)*sqrt(sh)*sqrt(qmsh)-
2372 2.*atan(sqrt(sh)*sqrt(qmsh)/dmsh))*(atan(sqrt(qmsh/sh))-atan(sqrt(sh)*sqrt(qmsh)/dmsh))));
2373}
2374
2375double BXqll::tau89(double sh)
2376{
2377 gslpp::complex i = gslpp::complex::i();
2378 double dmsh = 2.-sh;
2379 double qmsh = 4.-sh;
2380
2381 return (2./3.*(sh*qmsh-3.-4.*log(sh)*(1.-sh-sh*sh)-8.*(dilog(sh/2.+i*sqrt(sh)*sqrt(qmsh)/2.)-
2382 dilog((sh*qmsh-2.)/2.+i*dmsh*sqrt(sh)*sqrt(qmsh)/2.)).real()+
2383 4.*(sh*sh*sqrt(qmsh/sh)+2.*atan(sqrt(sh)*sqrt(qmsh)/dmsh))*(atan(sqrt(qmsh/sh))-
2384 atan(sqrt(sh)*sqrt(qmsh)/dmsh))));
2385}
2386
2387double BXqll::tau88(double sh)
2388{
2389 gslpp::complex i = gslpp::complex::i();
2390 double umsh = 1.-sh;
2391 double qmsh = 4.-sh;
2392
2393 return (4./27./sh*(-8.*M_PI*M_PI+umsh*(77.-sh-4.*sh*sh)-24.*dilog((gslpp::complex) umsh).real()+
2394 3.*(10.-4.*sh-9.*sh*sh+8.*log(sqrt(sh)/umsh))*log(sh)+48.*(dilog((3.-sh)/2.+
2395 i*umsh*sqrt(qmsh)/2./sqrt(sh))).real()-6.*((20.*sh+10.*sh*sh-3.*sh*sh*sh)/sqrt(sh)/sqrt(qmsh)
2396 -8.*M_PI+8.*atan(sqrt(qmsh/sh)))*(atan(sqrt(qmsh/sh))-atan(sqrt(sh)*sqrt(qmsh)/(2.-sh)))));
2397}
2398
2399double BXqll::tau22fit(double sh)
2400{
2401 double c0, c1, c2, c3, c4, c5, c6, c7;
2402 double q2 = sh * Mb_pole * Mb_pole;
2403 double sh2 = sh*sh;
2404 double sh3 = sh*sh2;
2405
2406 if (q2 <= 6.)
2407 {
2408 c0 = -85.5881372; c1 = 597.0445584; c2 = -3951.5752864; c3 = 19597.1776773;
2409 c4 = -61781.6106761; c5 = 109385.1435323; c6 = -82410.437659; c7 = -23.4496345;
2410 }
2411 else if (q2 >= 14.4)
2412 {
2413 c0 = -606.8557516; c1 = 1655.5780054, c2 = -2473.7140331; c3 = 2628.5002517;
2414 c4 = -1766.6491578; c5 = 674.9172919; c6 = -111.776621; c7 = -230.9832054;
2415 }
2416 else
2417 throw std::runtime_error("BXqll::tau22fit: region of q^2 not implemented");
2418
2419 return (c0 + c1 * sh + c2 * sh2 + c3 * sh3 + c4 * sh2*sh2 + c5 * sh2*sh3 + c6 * sh3*sh3 + c7 * log(sh));
2420}
2421
2422double BXqll::tau27fit_Re(double sh)
2423{
2424 double c0, c1, c2, c3, c4, c5, c6, c7;
2425 double q2 = sh * Mb_pole * Mb_pole;
2426 double sh2 = sh*sh;
2427 double sh3 = sh*sh2;
2428
2429 if (q2 <= 6.)
2430 {
2431 c0 = -110.2528109; c1 = 778.1558706; c2 = -5200.7693192; c3 = 25864.3213911;
2432 c4 = -81610.5151368; c5 = 144546.094428; c6 = -108930.3093624; c7 = -31.0505015;
2433 }
2434 else if (q2 >= 14.4)
2435 {
2436 c0 = 104.4988411; c1 = -236.2534692; c2 = 265.8226647; c3 = -224.7002307;
2437 c4 = 126.6513245; c5 = -42.3215314; c6 = 6.3024015; c7 = 45.8968367;
2438 }
2439 else
2440 throw std::runtime_error("BXqll::tau27fit_Re: region of q^2 not implemented");
2441
2442 return (c0 + c1 * sh + c2 * sh2 + c3 * sh3 + c4 * sh2*sh2 + c5 * sh2*sh3 + c6 * sh3*sh3 + c7 * log(sh));
2443}
2444
2445double BXqll::tau27fit_Im(double sh)
2446{
2447 double c0, c1, c2, c3, c4, c5, c6, c7;
2448 double q2 = sh * Mb_pole * Mb_pole;
2449 double sh2 = sh*sh;
2450 double sh3 = sh*sh2;
2451
2452 if (q2 <= 6.)
2453 {
2454 c0 = -109.7117334; c1 = 770.3496716; c2 = -5092.660051; c3 = 25238.4761864;
2455 c4 = -79546.8013633; c5 = 140825.6087763; c6 = -106089.2269733; c7 = -30.1707767;
2456 }
2457 else if (q2 >= 14.4)
2458 {
2459 c0 = -3232.2938841; c1 = 8748.3523192; c2 = -12941.1559426; c3 = 13657.4529099;
2460 c4 = -9136.8371605; c5 = 3479.428199; c6 = -574.9465146; c7 = -1238.52469;
2461 }
2462 else
2463 throw std::runtime_error("BXqll::tau27fit_Im: region of q^2 not implemented");
2464
2465 return (c0 + c1 * sh + c2 * sh2 + c3 * sh3 + c4 * sh2*sh2 + c5 * sh2*sh3 + c6 * sh3*sh3 + c7 * log(sh));
2466}
2467
2468double BXqll::tau28fit_Re(double sh)
2469{
2470 double c0, c1, c2, c3, c4, c5, c6, c7;
2471 double q2 = sh * Mb_pole * Mb_pole;
2472 double sh2 = sh*sh;
2473 double sh3 = sh*sh2;
2474
2475 if (q2 <= 6.)
2476 {
2477 c0 = 31.4802423; c1 = -241.0846881; c2 = 1656.2061824; c3 = -8298.8096208;
2478 c4 = 26278.5739013; c5 = -46638.1094866; c6 = 35189.3746188; c7 = 7.7075461;
2479 }
2480 else if (q2 >= 14.4)
2481 {
2482 c0 = 4.6900829; c1 = -11.9199517; c2 = 15.9522495; c3 = -15.3178755;
2483 c4 = 9.3880956; c5 = -3.2996681; c6 = 0.5070672; c7 = 1.8726389;
2484 }
2485 else
2486 throw std::runtime_error("BXqll::tau28fit_Re: region of q^2 not implemented");
2487
2488 return (c0 + c1 * sh + c2 * sh2 + c3 * sh3 + c4 * sh2*sh2 + c5 * sh2*sh3 + c6 * sh3*sh3 + c7 * log(sh));
2489}
2490
2491double BXqll::tau28fit_Im(double sh)
2492{
2493 double c0, c1, c2, c3, c4, c5, c6, c7;
2494 double q2 = sh * Mb_pole * Mb_pole;
2495 double sh2 = sh*sh;
2496 double sh3 = sh*sh2;
2497
2498 if (q2 <= 6.)
2499 {
2500 c0 = 33.7994469; c1 = -246.7751255; c2 = 1660.6205985; c3 = -8272.13013;
2501 c4 = 26130.6171608; c5 = -46315.2850273; c6 = 34917.5988012; c7 = 8.6664919;
2502 }
2503 else if (q2 >= 14.4)
2504 {
2505 c0 = -36.7487481; c1 = 99.7027369; c2 = -147.7606619; c3 = 155.8859306;
2506 c4 = -104.1170313; c5 = 39.5573846; c6 = -6.5196117; c7 = -14.0404657;
2507 }
2508 else
2509 throw std::runtime_error("BXqll::tau28fit_Im: region of q^2 not implemented");
2510
2511 return (c0 + c1 * sh + c2 * sh2 + c3 * sh3 + c4 * sh2*sh2 + c5 * sh2*sh3 + c6 * sh3*sh3 + c7 * log(sh));
2512}
2513
2514double BXqll::tau29fit_Re(double sh)
2515{
2516 double c0, c1, c2, c3, c4, c5, c6, c7;
2517 double q2 = sh * Mb_pole * Mb_pole;
2518 double sh2 = sh*sh;
2519 double sh3 = sh*sh2;
2520
2521 if (q2 <= 6.)
2522 {
2523 c0 = 0.6273561; c1 = -0.6332694; c2 = -2.9264173; c3 = 8.9668433;
2524 c4 = -17.5473592; c5 = 25.1466317; c6 = -19.5875632; c7 = 0.0002596;
2525 }
2526 else if (q2 >= 14.4)
2527 {
2528 c0 = 5.6720573; c1 = -8.2946135; c2 = 2.2577422; c3 = 1.1204784;
2529 c4 = -0.9164538; c5 = 0.1768918; c6 = -0.0161023; c7 = 3.2956642;
2530 }
2531 else
2532 throw std::runtime_error("BXqll::tau29fit_Re: region of q^2 not implemented");
2533
2534 return (c0 + c1 * sh + c2 * sh2 + c3 * sh3 + c4 * sh2*sh2 + c5 * sh2*sh3 + c6 * sh3*sh3 + c7 * log(sh));
2535}
2536
2537double BXqll::tau29fit_Im(double sh)
2538{
2539 double c0, c1, c2, c3, c4, c5, c6, c7;
2540 double q2 = sh * Mb_pole * Mb_pole;
2541 double sh2 = sh*sh;
2542 double sh3 = sh*sh2;
2543
2544 if (q2 <= 6.)
2545 {
2546 c0 = 0.3394312; c1 = -0.2742808; c2 = -0.1496318; c3 = 1.6480602;
2547 c4 = -6.3928038; c5 = 13.9968788; c6 = -11.1080801; c7 = -0.0001847;
2548 }
2549 else if (q2 >= 14.4)
2550 {
2551 c0 = -823.7736; c1 = 2231.2869412; c2 = -3306.5226002; c3 = 3497.1579714;
2552 c4 = -2344.7802726; c5 = 894.7382317; c6 = -148.1066909; c7 = -315.648791;
2553 }
2554 else
2555 throw std::runtime_error("BXqll::tau29fit_Im: region of q^2 not implemented");
2556
2557 return (c0 + c1 * sh + c2 * sh2 + c3 * sh3 + c4 * sh2*sh2 + c5 * sh2*sh3 + c6 * sh3*sh3 + c7 * log(sh));
2558}
2559
2560double BXqll::t810(double sh)
2561{
2562 gslpp::complex i = gslpp::complex::i();
2563 double umsh = 1. - sh;
2564 double shmt = sh - 3.;
2565 double ssh = sqrt(sh);
2566 double sqmsh = sqrt(4. - sh);
2567 double umsqrt = 1. - ssh;
2568 double dmsqrt = 2. - ssh;
2569
2570 double dilog_1 = gslpp_special_functions::dilog(-ssh);
2571 double dilog_2 = (i * dilog((-2. + i * sqmsh + ssh) * ssh / (i * sqmsh - ssh))).real();
2572 double dilog_3 = (i * dilog(i * sqmsh * umsh * ssh / 2. - shmt * sh / 2.)).real();
2573
2574 return (3. * (umsqrt * umsqrt * (23. - 6. * ssh - sh) + 4. * umsh * (7. + sh) * log(1. + ssh) +
2575 2. * sh * (1. + sh - log(sh)) * log(sh)) + 2. * (-3. * MPI2 * (1. + 2. * sh) -
2576 6. * shmt * sh * log(dmsqrt) - 36. * (1. + 2. * sh) * dilog_1 - ssh / sqmsh *
2577 (6. * (2. * shmt * sh * atan((2. + ssh) / sqmsh) + 2. * M_PI * log(dmsqrt) -
2578 atan(sqmsh / ssh) * (shmt * sh + 4. * log(dmsqrt)) - atan(ssh * sqmsh / (2. - sh)) *
2579 (shmt * sh - log(sh)) + 4. * dilog_2 - 2. * dilog_3)))) / 6. / umsh / umsh;
2580}
2581
2582double BXqll::t210fit(double sh)
2583{
2584 double c0, c1, c2, c3, c4, c5, c6, c7;
2585 double q2 = sh * Mb_pole * Mb_pole;
2586 double sh2 = sh*sh;
2587 double sh3 = sh*sh2;
2588
2589 if (q2 <= 6.)
2590 {
2591 c0 = 0.4853155; c1 = -0.4841452; c2 = 2.1099434; c3 = -17.9649531;
2592 c4 = 69.0927112; c5 = -154.3567398; c6 = 145.5137293; c7 = -0.001321;
2593 }
2594 else if (q2 >= 14.4)
2595 {
2596 c0 = 0.84111; c1 = -0.2282341; c2 = -1.3546028; c3 = 1.3081906;
2597 c4 = -0.8150979; c5 = 0.2963146; c6 = -0.0476803; c7 = 1.0777667;
2598 }
2599 else
2600 throw std::runtime_error("BXqll::t210fit_Re: region of q^2 not implemented");
2601
2602 return (c0 + c1 * sh + c2 * sh2 + c3 * sh3 + c4 * sh2*sh2 + c5 * sh2*sh3 + c6 * sh3*sh3 + c7 * log(sh));
2603}
FULLNNNLO
@ FULLNNNLO
Definition: OrderScheme.h:40
NNLO
@ NNLO
Definition: OrderScheme.h:36
LO
@ LO
Definition: OrderScheme.h:34
NLO
@ NLO
Definition: OrderScheme.h:35
FULLNNLO
@ FULLNNLO
Definition: OrderScheme.h:39
FULLNLO
@ FULLNLO
Definition: OrderScheme.h:38
QED1
@ QED1
Definition: OrderScheme.h:92
QED0
@ QED0
Definition: OrderScheme.h:91
QED2
@ QED2
Definition: OrderScheme.h:93
NLO_QED11
@ NLO_QED11
Definition: OrderScheme.h:59
FULLNLO_QED
@ FULLNLO_QED
Definition: OrderScheme.h:64
LO_QED
@ LO_QED
Definition: OrderScheme.h:58
NLO_QED21
@ NLO_QED21
Definition: OrderScheme.h:60
NO_QED
@ NO_QED
Definition: OrderScheme.h:57
QCD1
@ QCD1
Definition: OrderScheme.h:76
QCD0
@ QCD0
Definition: OrderScheme.h:75
QCD2
@ QCD2
Definition: OrderScheme.h:77
BXqll::tau88
double tau88(double sh)
The finite bremsstrahlung correction from .
Definition: BXqll.cpp:2387
BXqll::getH
double getH(std::string obs, double sh)
Method to obtain each observable as defined in .
Definition: BXqll.cpp:277
BXqll::omega22em_T
double omega22em_T(double sh)
Definition: BXqll.cpp:1389
BXqll::phinv_00_01
double phinv_00_01
Definition: BXqll.h:89
BXqll::z
double z
Definition: BXqll.h:87
BXqll::tau28fit_Re
double tau28fit_Re(double sh)
The fit of the real part of finite bremsstrahlung correction from .
Definition: BXqll.cpp:2468
BXqll::Mb
double Mb
Definition: BXqll.h:86
BXqll::FH
gsl_function FH
Definition: BXqll.h:110
BXqll::KS_cc
gslpp::complex KS_cc(double sh)
Kruger-Sehgal factorizable non-perturbative charm contributions following .
Definition: BXqll.cpp:1561
BXqll::omega79_L
double omega79_L(double sh)
Definition: BXqll.cpp:1298
BXqll::phinv11
double phinv11
Definition: BXqll.h:89
BXqll::QED_max
unsigned int QED_max
Definition: BXqll.h:91
BXqll::PhiTL_brems
double PhiTL_brems(double sh)
The finite bremsstrahlung corrections to dGamma/ds for from .
Definition: BXqll.cpp:2281
BXqll::h_z
gslpp::complex h_z(double zed, double sh)
Auxiliary function from .
Definition: BXqll.cpp:2339
BXqll::M_10
std::vector< gslpp::vector< gslpp::complex > > M_10
Definition: BXqll.h:95
BXqll::Phi_u
double Phi_u(orders ord)
Normalization function for from eq. (4.8) of 1503.04849.
Definition: BXqll.cpp:1610
BXqll::Mc
double Mc
Definition: BXqll.h:86
BXqll::H_A
double H_A(double sh)
Angular observable as defined in .
Definition: BXqll.cpp:694
BXqll::phi00
double phi00
Definition: BXqll.h:88
BXqll::alstilde
double alstilde
Definition: BXqll.h:85
BXqll::S79_T
double S79_T(double sh, orders order)
Definition: BXqll.cpp:889
BXqll::F29
gslpp::complex F29(double sh)
The correction from .
Definition: BXqll.cpp:153
BXqll::omega77em_T
double omega77em_T(double sh)
Auxiliary functions from .
Definition: BXqll.cpp:1365
BXqll::eij_L
gslpp::complex eij_L(unsigned int i, unsigned int j, double sh)
Definition: BXqll.cpp:1177
BXqll::phinv20
double phinv20
Definition: BXqll.h:89
BXqll::f_Huber
gslpp::complex f_Huber(double sh, double gamma_9, double rho_c, double rho_b, double rho_0, double rho_num)
Auxiliary function from .
Definition: BXqll.cpp:1522
BXqll::BXqll
BXqll(const StandardModel &SM_i, QCD::quark quark_i, QCD::lepton lep_i)
Constructor.
Definition: BXqll.cpp:21
BXqll::F_27re
double F_27re(double muh, double z, double sh, int maxpow=20)
Definition: BXqll.cpp:215
BXqll::integrateH
double integrateH(std::string obs, double q_min, double q_max)
The integral of each observable as defined in .
Definition: BXqll.cpp:260
BXqll::omega79_T
double omega79_T(double sh)
Definition: BXqll.cpp:1257
BXqll::F_19re
double F_19re(double muh, double z, double sh, int maxpow=20)
Definition: BXqll.cpp:205
BXqll::S910_A
double S910_A(double sh, orders order)
Definition: BXqll.cpp:1021
BXqll::DeltaF_29im
double DeltaF_29im(double muh, double z, double sh, int maxpow=20)
Definition: BXqll.cpp:250
BXqll::H_T
double H_T(double sh)
Angular observable as defined in .
Definition: BXqll.cpp:296
BXqll::alsmuc
double alsmuc
Definition: BXqll.h:85
BXqll::kappa
double kappa
Definition: BXqll.h:85
BXqll::omega27em_T
gslpp::complex omega27em_T(double sh)
Definition: BXqll.cpp:1401
BXqll::M_9
std::vector< gslpp::vector< gslpp::complex > > M_9
Definition: BXqll.h:94
BXqll::S77_L
double S77_L(double sh, orders order)
Auxiliary functions from .
Definition: BXqll.cpp:936
BXqll::S77_T
double S77_T(double sh, orders order)
Auxiliary functions from .
Definition: BXqll.cpp:870
BXqll::cij_A
gslpp::complex cij_A(unsigned int i, unsigned int j, double sh)
Definition: BXqll.cpp:1120
BXqll::Mlep
double Mlep
Definition: BXqll.h:86
BXqll::H_L
double H_L(double sh)
Angular observable as defined in .
Definition: BXqll.cpp:495
BXqll::F87
gslpp::complex F87(double sh)
The correction from .
Definition: BXqll.cpp:174
BXqll::g_Huber
gslpp::complex g_Huber(double y)
Auxiliary function from .
Definition: BXqll.cpp:1545
BXqll::computeMi
void computeMi(double sh)
Vectors of auxiliary functions from Table 6 of .
Definition: BXqll.cpp:803
BXqll::lambda_1
double lambda_1
Definition: BXqll.h:87
BXqll::S99_L
double S99_L(double sh, orders order)
Definition: BXqll.cpp:974
BXqll::ale
double ale
Definition: BXqll.h:85
BXqll::phi2
double phi2
Definition: BXqll.h:88
BXqll::Hij_L
std::vector< gslpp::matrix< gslpp::complex > > Hij_L
Definition: BXqll.h:98
BXqll::tau27fit_Re
double tau27fit_Re(double sh)
The fit of the real part of finite bremsstrahlung correction from .
Definition: BXqll.cpp:2422
BXqll::C_ratio
double C_ratio
Definition: BXqll.h:90
BXqll::tau27fit_Im
double tau27fit_Im(double sh)
The fit of the imaginary part of finite bremsstrahlung correction from .
Definition: BXqll.cpp:2445
BXqll::alsmu
double alsmu
Definition: BXqll.h:85
BXqll::QCD_max
unsigned int QCD_max
Definition: BXqll.h:91
BXqll::aletilde
double aletilde
Definition: BXqll.h:85
BXqll::omega77_L
double omega77_L(double sh)
Auxiliary functions from .
Definition: BXqll.cpp:1284
BXqll::Test_WC_DF1
void Test_WC_DF1()
Temporary method to test Wilson coefficients with C10_OS1 matching and HeffDF1 evolution.
Definition: BXqll.cpp:1712
BXqll::F17
gslpp::complex F17(double sh)
The correction from .
Definition: BXqll.cpp:160
BXqll::eij_T
gslpp::complex eij_T(unsigned int i, unsigned int j, double sh)
Log-enhanced electromagnetic corrections as defined in .
Definition: BXqll.cpp:1137
BXqll::mu_c
double mu_c
Definition: BXqll.h:86
BXqll::tau28fit_Im
double tau28fit_Im(double sh)
The fit of the imaginary part of finite bremsstrahlung correction from .
Definition: BXqll.cpp:2491
BXqll::S99_T
double S99_T(double sh, orders order)
Definition: BXqll.cpp:908
BXqll::Mtau
double Mtau
Definition: BXqll.h:86
BXqll::updateParameters
void updateParameters()
The update parameter method for BXqll.
Definition: BXqll.cpp:48
BXqll::phinv00
double phinv00
Definition: BXqll.h:89
BXqll::F_27im
double F_27im(double muh, double z, double sh, int maxpow=20)
Definition: BXqll.cpp:220
BXqll::t810
double t810(double sh)
The finite bremsstrahlung correction from .
Definition: BXqll.cpp:2560
BXqll::cij_L
gslpp::complex cij_L(unsigned int i, unsigned int j, double sh, unsigned int ord)
Definition: BXqll.cpp:1082
BXqll::omega29em_L
gslpp::complex omega29em_L(double sh)
Definition: BXqll.cpp:1476
BXqll::w_H
gsl_integration_cquad_workspace * w_H
Definition: BXqll.h:111
BXqll::C9mod
gslpp::complex C9mod(double sh)
The modified coefficient from .
Definition: BXqll.cpp:2325
BXqll::omega79em_T
double omega79em_T(double sh)
Definition: BXqll.cpp:1373
BXqll::phinv01
double phinv01
Definition: BXqll.h:89
BXqll::lep
QCD::lepton lep
Definition: BXqll.h:83
BXqll::phi1
double phi1
Definition: BXqll.h:88
BXqll::omega79em_L
double omega79em_L(double sh)
Definition: BXqll.cpp:1435
BXqll::omega22em_L
double omega22em_L(double sh)
Definition: BXqll.cpp:1451
BXqll::cij_T
gslpp::complex cij_T(unsigned int i, unsigned int j, double sh, unsigned int ord)
contributions as defined in
Definition: BXqll.cpp:1044
BXqll::omega77em_L
double omega77em_L(double sh)
Auxiliary functions from .
Definition: BXqll.cpp:1427
BXqll::omega210em_A
gslpp::complex omega210em_A(double sh)
Definition: BXqll.cpp:1505
BXqll::CCH_multiplication
double CCH_multiplication(std::vector< gslpp::matrix< gslpp::complex > > &Hij)
Auxiliary function that performs the multiplication of Wilson coefficients and matrix elements.
Definition: BXqll.cpp:1668
BXqll::lambda_2
double lambda_2
Definition: BXqll.h:87
BXqll::MW
double MW
Definition: BXqll.h:86
BXqll::omega710_A
double omega710_A(double sh)
Auxiliary functions from .
Definition: BXqll.cpp:1325
BXqll::tau78
double tau78(double sh)
The finite bremsstrahlung correction from .
Definition: BXqll.cpp:2362
BXqll::KS_aux
gslpp::complex KS_aux(double sh, double m, double Gamma, double Br_ll, double Br_had)
Auxiliary function for the Kruger-Sehgal charm contributions.
Definition: BXqll.cpp:1582
BXqll::F_29im
double F_29im(double muh, double z, double sh, int maxpow=20)
Definition: BXqll.cpp:230
BXqll::F19
gslpp::complex F19(double sh)
The correction from .
Definition: BXqll.cpp:146
BXqll::phinv21
double phinv21
Definition: BXqll.h:89
BXqll::pre
double pre
Definition: BXqll.h:90
BXqll::myF_1
F_1 myF_1
Definition: BXqll.h:80
BXqll::Hij_T
std::vector< gslpp::matrix< gslpp::complex > > Hij_T
Definition: BXqll.h:97
BXqll::DeltaF_19re
double DeltaF_19re(double muh, double z, double sh, int maxpow=20)
Definition: BXqll.cpp:235
BXqll::F_17re
double F_17re(double muh, double z, double sh, int maxpow=20)
Definition: BXqll.cpp:195
BXqll::myHeff
HeffDF1 myHeff
Definition: BXqll.h:82
BXqll::tau22fit
double tau22fit(double sh)
The fit of the finite bremsstrahlung correction from .
Definition: BXqll.cpp:2399
BXqll::WC
gslpp::matrix< gslpp::complex > WC
Definition: BXqll.h:106
BXqll::tau89
double tau89(double sh)
The finite bremsstrahlung correction from .
Definition: BXqll.cpp:2375
BXqll::S1010_L
double S1010_L(double sh, orders order)
Definition: BXqll.cpp:997
BXqll::omega99em_T
double omega99em_T(double sh)
Definition: BXqll.cpp:1381
BXqll::S79_L
double S79_L(double sh, orders order)
Definition: BXqll.cpp:955
BXqll::phi01
double phi01
Definition: BXqll.h:88
BXqll::abslambdat_over_Vcb
double abslambdat_over_Vcb
Definition: BXqll.h:87
BXqll::mu_b
double mu_b
Definition: BXqll.h:86
BXqll::BXqllParameters
std::vector< std::string > BXqllParameters
Definition: BXqll.h:101
BXqll::Lbl
double Lbl
Definition: BXqll.h:87
BXqll::initializeBXqllParameters
std::vector< std::string > initializeBXqllParameters()
A method for initializing the parameters necessary for BXqll.
Definition: BXqll.cpp:41
BXqll::int_qed
unsigned int int_qed(orders_qed order_qed)
Auxiliary function that matches orders_qed to an integer.
Definition: BXqll.cpp:1662
BXqll::mySM
const StandardModel & mySM
Definition: BXqll.h:79
BXqll::omega77_T
double omega77_T(double sh)
Auxiliary functions from .
Definition: BXqll.cpp:1243
BXqll::DeltaF_19im
double DeltaF_19im(double muh, double z, double sh, int maxpow=20)
Definition: BXqll.cpp:240
BXqll::omega27em_L
gslpp::complex omega27em_L(double sh)
Definition: BXqll.cpp:1463
BXqll::tau29fit_Re
double tau29fit_Re(double sh)
The fit of the real part of finite bremsstrahlung correction from .
Definition: BXqll.cpp:2514
BXqll::phi20
double phi20
Definition: BXqll.h:88
BXqll::BR_BXcenu
double BR_BXcenu
Definition: BXqll.h:90
BXqll::omega910em_A
double omega910em_A(double sh)
Definition: BXqll.cpp:1497
BXqll::CF
double CF
Definition: BXqll.h:85
BXqll::Mc_pole
double Mc_pole
Definition: BXqll.h:86
BXqll::F27
gslpp::complex F27(double sh)
The correction from .
Definition: BXqll.cpp:167
BXqll::F_19im
double F_19im(double muh, double z, double sh, int maxpow=20)
Definition: BXqll.cpp:210
BXqll::F89
double F89(double sh)
The correction from .
Definition: BXqll.cpp:185
BXqll::M_7
std::vector< gslpp::vector< gslpp::complex > > M_7
Definition: BXqll.h:93
BXqll::F_17im
double F_17im(double muh, double z, double sh, int maxpow=20)
Definition: BXqll.cpp:200
BXqll::tau29fit_Im
double tau29fit_Im(double sh)
The fit of the imaginary part of finite bremsstrahlung correction from .
Definition: BXqll.cpp:2537
BXqll::S710_A
double S710_A(double sh, orders order)
Auxiliary functions from .
Definition: BXqll.cpp:1002
BXqll::F_BIR
gslpp::complex F_BIR(double r)
Auxiliary function from .
Definition: BXqll.cpp:1597
BXqll::dilog
friend double gslpp_special_functions::dilog(double x)
BXqll::~BXqll
virtual ~BXqll()
Destructor.
Definition: BXqll.cpp:37
BXqll::F_29re
double F_29re(double muh, double z, double sh, int maxpow=20)
Definition: BXqll.cpp:225
BXqll::getR_LOWQ2
double getR_LOWQ2(double sh)
dGamma/ds for in the low dilepton invariant mass region.
Definition: BXqll.cpp:134
BXqll::eij_A
gslpp::complex eij_A(unsigned int i, unsigned int j, double sh)
Definition: BXqll.cpp:1217
BXqll::aveH
double aveH
Definition: BXqll.h:108
BXqll::PhiA_brems
double PhiA_brems(double sh)
The finite bremsstrahlung corrections to dAFB/ds for from .
Definition: BXqll.cpp:2315
BXqll::omega910_A
double omega910_A(double sh)
Definition: BXqll.cpp:1351
BXqll::old_handler
gsl_error_handler_t * old_handler
Definition: BXqll.h:112
BXqll::allcoeff
Expanded< gslpp::vector< gslpp::complex > > allcoeff
Definition: BXqll.h:104
BXqll::f9pen_Huber
gslpp::complex f9pen_Huber(double sh)
Auxiliary function from .
Definition: BXqll.cpp:1534
BXqll::muh
double muh
Definition: BXqll.h:87
BXqll::Hij_A
std::vector< gslpp::matrix< gslpp::complex > > Hij_A
Definition: BXqll.h:99
BXqll::t210fit
double t210fit(double sh)
The fit of finite bremsstrahlung correction from .
Definition: BXqll.cpp:2582
BXqll::omega99em_L
double omega99em_L(double sh)
Definition: BXqll.cpp:1443
BXqll::DeltaF_29re
double DeltaF_29re(double muh, double z, double sh, int maxpow=20)
Definition: BXqll.cpp:245
BXqll::Mb_pole
double Mb_pole
Definition: BXqll.h:86
BXqll::phi00_2
double phi00_2
Definition: BXqll.h:88
BXqll::myF_2
F_2 myF_2
Definition: BXqll.h:81
BXqll::quark
QCD::quark quark
Definition: BXqll.h:84
BXqll::GF
double GF
Definition: BXqll.h:85
BXqll::omega99_L
double omega99_L(double sh)
Definition: BXqll.cpp:1311
BXqll::omega99_T
double omega99_T(double sh)
Definition: BXqll.cpp:1270
BXqll::phinv10
double phinv10
Definition: BXqll.h:89
BXqll::omega710em_A
double omega710em_A(double sh)
Auxiliary functions from .
Definition: BXqll.cpp:1489
BXqll::omega29em_T
gslpp::complex omega29em_T(double sh)
Definition: BXqll.cpp:1414
BXqll::Phi_u_inv
double Phi_u_inv(unsigned int ord_qcd, unsigned int ord_qed)
Inverse of the normalization function for from eq. (4.8) of 1503.04849.
Definition: BXqll.cpp:1644
BXqll::errH
double errH
Definition: BXqll.h:109
BXqll::Ms
double Ms
Definition: BXqll.h:86
BXqll::S1010_T
double S1010_T(double sh, orders order)
Definition: BXqll.cpp:931
CKM::computelamt_s
const gslpp::complex computelamt_s() const
The product of the CKM elements .
Definition: CKM.cpp:174
CKM::getV_cb
const gslpp::complex getV_cb() const
A member for returning the value of the CKM element .
Definition: CKM.h:247
F_1::F_17im
double F_17im(double muh, double z, double sh, int maxpow=20)
Definition: F_1.cpp:1889
F_1::DeltaF_19re
double DeltaF_19re(double muh, double z, double sh, int maxpow=20)
Definition: F_1.cpp:5296
F_1::F_19im
double F_19im(double muh, double z, double sh, int maxpow=20)
Definition: F_1.cpp:4628
F_1::F_19re
double F_19re(double muh, double z, double sh, int maxpow=20)
Definition: F_1.cpp:2475
F_1::F_17re
double F_17re(double muh, double z, double sh, int maxpow=20)
Definition: F_1.cpp:21
F_1::DeltaF_19im
double DeltaF_19im(double muh, double z, double sh, int maxpow=20)
Definition: F_1.cpp:5438
F_2::F_27re
double F_27re(double muh, double z, double sh, int maxpow=20)
Definition: F_2.cpp:21
F_2::F_27im
double F_27im(double muh, double z, double sh, int maxpow=20)
Definition: F_2.cpp:1850
F_2::DeltaF_29re
double DeltaF_29re(double muh, double z, double sh, int maxpow=20)
Definition: F_2.cpp:5177
F_2::DeltaF_29im
double DeltaF_29im(double muh, double z, double sh, int maxpow=20)
Definition: F_2.cpp:5314
F_2::F_29re
double F_29re(double muh, double z, double sh, int maxpow=20)
Definition: F_2.cpp:2435
F_2::F_29im
double F_29im(double muh, double z, double sh, int maxpow=20)
Definition: F_2.cpp:4517
HeffDF1::ComputeCoeff
virtual Expanded< gslpp::vector< gslpp::complex > > ComputeCoeff(double mu, schemes scheme=NDR)
Definition: HeffDF1.cpp:110
HeffDF1::LowScaleCoeff
gslpp::vector< gslpp::complex > LowScaleCoeff(qcd_orders order_qcd, qed_orders order_qed)
Definition: HeffDF1.cpp:23
Particle::getMass
const double & getMass() const
A get method to access the particle mass.
Definition: Particle.h:61
QCD::getMuc
const double getMuc() const
A get method to access the threshold between four- and three-flavour theory in GeV.
Definition: QCD.h:591
QCD::getCF
const double getCF() const
A get method to access the Casimir factor of QCD.
Definition: QCD.h:618
QCD::getOptionalParameter
const double getOptionalParameter(std::string name) const
A method to get parameters that are specific to only one set of observables.
Definition: QCD.h:450
QCD::quark
quark
An enum type for quarks.
Definition: QCD.h:323
QCD::BOTTOM
@ BOTTOM
Definition: QCD.h:329
QCD::STRANGE
@ STRANGE
Definition: QCD.h:327
QCD::CHARM
@ CHARM
Definition: QCD.h:326
QCD::Beta0
const double Beta0(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:601
QCD::lepton
lepton
An enum type for leptons.
Definition: QCD.h:310
QCD::TAU
@ TAU
Definition: QCD.h:316
QCD::getQuarks
const Particle & getQuarks(const QCD::quark q) const
A get method to access a quark as an object of the type Particle.
Definition: QCD.h:536
QCD::getMub
const double getMub() const
A get method to access the threshold between five- and four-flavour theory in GeV.
Definition: QCD.h:582
QCD::Mbar2Mp
const double Mbar2Mp(const double mbar, const quark q, const orders order=FULLNNLO) const
Converts the mass to the pole mass.
Definition: QCD.cpp:1552
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:509
StandardModel::getMuw
const double getMuw() const
A get method to retrieve the matching scale around the weak scale.
Definition: StandardModel.h:951
StandardModel::getLeptons
const Particle & getLeptons(const QCD::lepton p) const
A get method to retrieve the member object of a lepton.
Definition: StandardModel.h:744
StandardModel::Ale
const double Ale(double mu, orders order, bool Nf_thr=true) const
The running electromagnetic coupling in the scheme.
Definition: StandardModel/src/StandardModel.cpp:777
StandardModel::getCKM
const CKM & getCKM() const
A get method to retrieve the member object of type CKM.
Definition: StandardModel.h:920
StandardModel::Mw
virtual const double Mw() const
The SM prediction for the -boson mass in the on-shell scheme, .
Definition: StandardModel/src/StandardModel.cpp:1019
StandardModel::getGF
const double getGF() const
A get method to retrieve the Fermi constant .
Definition: StandardModel.h:780
StandardModel::Als
const double Als(const double mu, const orders order, const bool Nf_thr, const bool qed_flag) const
The running QCD coupling in the scheme including QED corrections.
Definition: StandardModel.h:1088
StandardModel::alphaMz
virtual const double alphaMz() const
The electromagnetic coupling at the -mass scale, .
Definition: StandardModel/src/StandardModel.cpp:938
orders
orders
An enum type for orders in QCD.
Definition: OrderScheme.h:33
qed_orders
qed_orders
An enum type for qed_orders in electroweak.
Definition: OrderScheme.h:90
orders_qed
orders_qed
An enum type for orders in electroweak.
Definition: OrderScheme.h:56
qcd_orders
qcd_orders
An enum type for qcd_orders in QCD.
Definition: OrderScheme.h:74