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