a Code for the Combination of Indirect and Direct Constraints on High Energy Physics Models Logo
bsgamma.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2015 HEPfit Collaboration
3  *
4  *
5  * For the licensing terms see doc/COPYING.
6  */
7 
8 /*
9  * phi1.4body partially hardcoded, currently switched off by setting FOUR_BODY to false
10  * EW matching missing
11  */
12 
13 #include "StandardModel.h"
14 #include "bsgamma.h"
15 #include "std_make_vector.h"
16 #include "gslpp_function_adapter.h"
17 #include <gsl/gsl_sf_dilog.h>
18 #include <gsl/gsl_sf_zeta.h>
19 #include <gsl/gsl_sf_clausen.h>
20 #include <boost/bind.hpp>
21 
22 Bsgamma::Bsgamma(const StandardModel& SM_i, QCD::quark quark_i, int obsFlag)
23 : ThObservable(SM_i),
24 Intbc_cache(2, 0.)
25 {
26  if (obsFlag > 0 and obsFlag < 3) obs = obsFlag;
27  else throw std::runtime_error("obsFlag in bsgamma can only be 1 (BR) or 2 (ACP)");
28 
29  quark = quark_i;
30 
31  SUM = false;
32  EWflag = true;
33  FOUR_BODY = false;
34  WET_NP_btos = false;
35  SMEFT_NP_btos = false;
36 
37  setParametersForObservable(make_vector<std::string>() << "mukin" << "BRsem" << "Mbkin" << "Mcatmuc" << "mupi2"
38  << "rhoD3" << "muG2" << "rhoLS3" << "BLNPcorr" << "mu_b_bsgamma" << "mu_c_bsgamma");
39 
40  Intb1Cached = 0;
41  Intb2Cached = 0;
42  Intb3Cached = 0;
43  Intb4Cached = 0;
44  Intbb1Cached = 0;
45  Intbb2Cached = 0;
46  Intbb4Cached = 0;
47  Intbc1Cached = 0;
48  Intbc2Cached = 0;
49  Intc1Cached = 0;
50  Intc2Cached = 0;
51  Intc3Cached = 0;
52  Intcc1Cached = 0;
53  IntPhi772rCached = 0;
54 
55  w_INT = gsl_integration_cquad_workspace_alloc (100);
56 }
57 
58 Bsgamma::Bsgamma(const StandardModel& SM_i, int obsFlag)
59 : ThObservable(SM_i),
60 Intbc_cache(2, 0.)
61 {
62  if (obsFlag > 0 and obsFlag < 3) obs = obsFlag;
63  else throw std::runtime_error("obsFlag in bsgamma can only be 1 (BR) or 2 (ACP)");
64 
65  SUM = true;
66  EWflag = true;
67  FOUR_BODY = false;
68 
69  Intb1Cached = 0;
70  Intb2Cached = 0;
71  Intb3Cached = 0;
72  Intb4Cached = 0;
73  Intbb1Cached = 0;
74  Intbb2Cached = 0;
75  Intbb4Cached = 0;
76  Intbc1Cached = 0;
77  Intbc2Cached = 0;
78  Intc1Cached = 0;
79  Intc2Cached = 0;
80  Intc3Cached = 0;
81  Intcc1Cached = 0;
82  IntPhi772rCached = 0;
83 
84  w_INT = gsl_integration_cquad_workspace_alloc (100);
85 }
86 
87 void Bsgamma::checkCache()
88 {
89  if (Mb_kin == Intb_cache)
90  Intb_updated = 1;
91  else {
92  Intb_cache = Mb_kin;
93  Intb_updated = 0;
94  }
95 
96  if (Mb_kin == Intbc_cache(0) && Mc == Intbc_cache(1))
97  Intbc_updated = 1;
98  else {
99  Intbc_cache(0) = Mb_kin;
100  Intbc_cache(1) = Mc;
101  Intbc_updated = 0;
102  }
103 }
104 
105 double Bsgamma::delta(double E0)
106 {
107  return 1. - 2.*E0/Mb_kin;
108 }
109 
110 double Bsgamma::rho(double E0)
111 {
112  double d=delta(E0);
113  double d4=d*d*d*d;
114 
115  return d + d4/6. + log(1. - d);
116 }
117 
118 double Bsgamma::omega(double E0)
119 {
120  double d=delta(E0);
121  double d2=d*d;
122  double d3=d2*d;
123  double d4=d3*d;
124 
125  return 3./2. * d2 - 2. * d3 + d4;
126 }
127 
128 double Bsgamma::T1(double E0,double t)
129 {
130  double d=delta(E0);
131  double d2=d*d;
132  double d3=d2*d;
133  double d4=d3*d;
134 
135  double Li2 = gsl_sf_dilog(d);
136 
137  return 109./18. * d + 17./18. * d2 - 191./108. * d3 + 23./16. * d4
138  + 79./18. * log(1. - d) - 5./3. * Li2
139  - (5./3. * rho(E0) + 2./9. * omega(E0)) * log(t*d)
140  /* + rho(E0)/9.*log(ms^5/(mu^4*md)) */;
141 }
142 
143 double Bsgamma::T2(double E0,double t)
144 {
145  double d=delta(E0);
146  double d2=d*d;
147  double d3=d2*d;
148  double d4=d3*d;
149 
150  double Li2 = gsl_sf_dilog(d);
151 
152  return 187./108. * d + 7./18. * d2 - 395./648. * d3 + 1181./2592. * d4
153  + 133./108. * log(1. - d) - Li2/2.
154  - (rho(E0)/2. + 2./27. * omega(E0)) * log(t*d)
155  /* + rho(E0)/9.*log(ms/mu) */;
156 }
157 
158 double Bsgamma::T3(double E0,double t)
159 {
160  double d=delta(E0);
161  double d2=d*d;
162  double d3=d2*d;
163  double d4=d3*d;
164 
165  double Li2 = gsl_sf_dilog(d);
166 
167  return 35./162. * d + 1./72. * d2 - 89./1944. * d3 + 341./7776. * d4
168  + 13./81. * log(1. - d) - Li2/18.
169  - (rho(E0)/18. + omega(E0)/162.) * log(t*d);
170 }
171 
172 double Bsgamma::P0_4body(double E0, double t)
173 {
174  gslpp::complex la_u =-CKMu;
175 
176  double A1sq =C1_0.abs2()*CKMusq;
177  double A2sq =C2_0.abs2()*CKMusq;
178 
179  double C13re = (C1_0.real()*C3_0.real() + C1_0.imag()*C3_0.imag());
180  double C14re = (C1_0.real()*C4_0.real() + C1_0.imag()*C4_0.imag());
181  double C15re = (C1_0.real()*C5_0.real() + C1_0.imag()*C5_0.imag());
182  double C16re = (C1_0.real()*C6_0.real() + C1_0.imag()*C6_0.imag());
183 
184  double C13im = (C1_0.real()*C3_0.imag() + C1_0.imag()*C3_0.real());
185  double C14im = (C1_0.real()*C4_0.imag() + C1_0.imag()*C4_0.real());
186  double C15im = (C1_0.real()*C5_0.imag() + C1_0.imag()*C5_0.real());
187  double C16im = (C1_0.real()*C6_0.imag() + C1_0.imag()*C6_0.real());
188 
189  double C23re = (C2_0.real()*C3_0.real() + C2_0.imag()*C3_0.imag());
190  double C24re = (C2_0.real()*C4_0.real() + C2_0.imag()*C4_0.imag());
191  double C25re = (C2_0.real()*C5_0.real() + C2_0.imag()*C5_0.imag());
192  double C26re = (C2_0.real()*C6_0.real() + C2_0.imag()*C6_0.imag());
193 
194  double C23im = (C2_0.real()*C3_0.imag() + C2_0.imag()*C3_0.real());
195  double C24im = (C2_0.real()*C4_0.imag() + C2_0.imag()*C4_0.real());
196  double C25im = (C2_0.real()*C5_0.imag() + C2_0.imag()*C5_0.real());
197  double C26im = (C2_0.real()*C6_0.imag() + C2_0.imag()*C6_0.real());
198 
199  double C13 = (C13re*la_u.real() - C13im*la_u.imag());
200  double C14 = (C14re*la_u.real() - C14im*la_u.imag());
201  double C15 = (C15re*la_u.real() - C15im*la_u.imag());
202  double C16 = (C16re*la_u.real() - C16im*la_u.imag());
203  double C23 = (C23re*la_u.real() - C23im*la_u.imag());
204  double C24 = (C24re*la_u.real() - C24im*la_u.imag());
205  double C25 = (C25re*la_u.real() - C25im*la_u.imag());
206  double C26 = (C26re*la_u.real() - C26im*la_u.imag());
207  double C33 = C3_0.abs2();
208  double C34 = (C3_0.real()*C4_0.real() + C3_0.imag()*C4_0.imag());
209  double C35 = (C3_0.real()*C5_0.real() + C3_0.imag()*C5_0.imag());
210  double C36 = (C3_0.real()*C6_0.real() + C3_0.imag()*C6_0.imag());
211  double C44 = C4_0.abs2();
212  double C45 = (C4_0.real()*C5_0.real() + C4_0.imag()*C5_0.imag());
213  double C46 = (C4_0.real()*C6_0.real() + C4_0.imag()*C6_0.imag());
214  double C55 = C5_0.abs2();
215  double C56 = (C5_0.real()*C6_0.real() + C5_0.imag()*C6_0.imag());
216  double C66 = C6_0.abs2();
217 
218  return (C33 + 20. * C35 + 2./9. * C44 + 40./9. * C46 + 136. * C55 + 272./9. * C66) * T1(E0,t) +
219 
220  (2./9. * A1sq + A2sq
221  + 8./9. * C13 - 4./27. * C14 + 128./9. * C15 - 64./27. * C16
222  + 2./3. * C23 + 8./9. * C24 + 32./3. * C25 + 128./9. * C26) * T2(E0,t) +
223 
224  (C33 + 8./3. * C34 + 32. * C35 + 128./3. * C36 - 2./9. * C44 + 128./3. * C45
225  - 64./9. * C46 + 256. * C55 + 2048./3 * C56 - 512./9. * C66) * T3(E0,t);
226 }
227 
228 double Bsgamma::zeta()
229 {
230  return Mc*Mc/Mb_kin/Mb_kin;
231 }
232 
233 gslpp::complex Bsgamma::a(double z)
234 {
235  double zeta3 = gsl_sf_zeta_int(3);
236 
237  double z2=z*z;
238  double z3=z2*z;
239  double z4=z3*z;
240  double z5=z4*z;
241  double z6=z5*z;
242 
243  double L=log(z);
244  double L2=L*L;
245  double L3=L2*L;
246 
247  double pi2=M_PI*M_PI;
248 
249  if (z == 1.) return 4.0859 + 4./9. * M_PI * gslpp::complex::i();
250  else return 16./9. * ( ( 5./2. - pi2/3. - 3.*zeta3
251  + ( 5./2. - 3./4.*pi2 )*L + L2/4. + L3/12. )*z
252  + ( 7./4. + 2./3.*pi2 - pi2*L/2. - L2/4. + L3/12. )*z2
253  + ( -7./6. - pi2/4. + 2*L - 3./4.*L2 )*z3
254  + ( 457./216. - 5./18*pi2 - L/72. - 5./6.*L2 )*z4
255  + ( 35101./8640. - 35./72.*pi2 - 185./144.*L - 35./24.*L2 )*z5
256  + ( 67801./8000. - 21./20.*pi2 - 3303./800.*L - 63./20.*L2 )*z6 +
257  gslpp::complex::i()*M_PI*( ( 2. - pi2/6. + L/2. + L2/2. )*z
258  + ( 1./2. - pi2/6. - L + L2/2. )*z2
259  + z3 + 5./9.*z4 + 49./72.*z5 + 231./200.*z6) );
260 }
261 
262 gslpp::complex Bsgamma::b(double z)
263 {
264  double z2=z*z;
265  double z3=z2*z;
266  double z4=z3*z;
267  double z5=z4*z;
268  double z6=z5*z;
269 
270  double L=log(z);
271  double L2=L*L;
272 
273  double pi2=M_PI*M_PI;
274 
275  if (z == 1.) return 0.0316 + 4./81. * M_PI * gslpp::complex::i();
276  else return -8./9. * ( ( -3. + pi2/6. - L )*z - 2./3.*pi2*pow(z,3./2.)
277  + ( 1./2. + pi2 -2.*L - L2/2. )*z2
278  + ( -25./12. - pi2/9. - 19./18.*L + 2.*L2 )*z3
279  + ( -1376./225. + 137./30.*L + 2.*L2 + 2./3.*pi2 )*z4
280  + ( -131317./11760. + 887./84.*L + 5.*L2 + 5./3.*pi2 )*z5
281  + ( -2807617./97200. + 16597./540.*L + 14.*L2 + 14./3.*pi2 )*z6 +
282  gslpp::complex::i()*M_PI*( -z + ( 1 - 2.*L )*z2
283  + ( -10./9. + 4./3.*L )*z3 + z4 + 2./3.*z5 + 7./9.*z6) );
284 }
285 
286 gslpp::complex Bsgamma::r1(int i, double z)
287 {
288  double Xb = -0.16844083981858157;
289 
290  switch(i){
291  case 1:
292  return 833./729. - (a(z) + b(z))/3. + 40./243.*gslpp::complex::i()*M_PI;
293  case 2:
294  return - 1666./243. + 2.*(a(z) + b(z)) - 80./81.*gslpp::complex::i()*M_PI;
295  case 3:
296  return 2392./243. + 8.*M_PI/3./sqrt(3.) + 32./9.*Xb - a(1.) + 2.*b(1.) + 56./81.*gslpp::complex::i()*M_PI;
297  case 4:
298  return -761./729. - 4.*M_PI/9./sqrt(3.) - 16./27.*Xb + a(1.)/6. + 5.*b(1.)/3. + 2.*b(z) - 148./243.*gslpp::complex::i()*M_PI;
299  case 5:
300  return 56680./243. + 32.*M_PI/3./sqrt(3.) + 128./9.*Xb - 16.*a(1.) + 32.*b(1.) + 896./81.*gslpp::complex::i()*M_PI;
301  case 6:
302  return 5710./729. - 16.*M_PI/9./sqrt(3.) - 64./27.*Xb - 10./3.*a(1.) + 44./3.*b(1.) + 12.*a(z) + 20.*b(z)
303  - 2296./243.*gslpp::complex::i()*M_PI;
304  case 8:
305  return 44./9. - 8./27.*M_PI*M_PI + 8./9.*gslpp::complex::i()*M_PI;
306  default:
307  std::stringstream out;
308  out << i;
309  throw std::runtime_error("Bsgamma::r1(): index " + out.str() + " not implemented");
310  }
311 }
312 
313 gslpp::complex Bsgamma::r1_ew(int i, double z)
314 {
315  double Xb = -0.16844083981858157;
316  double PI2 = M_PI*M_PI;
317  gslpp::complex iPI = gslpp::complex::i()*M_PI;
318 
319  switch(i){
320  case 1:
321  return 3332./2187. - 4.*(a(z) + b(z))/9. + 160./729.*iPI;
322  case 2:
323  return 833./729. - (a(z) + b(z))/3. + 40./243.*iPI;
324  case 3:
325  return 748./729. + 2.*M_PI/9./sqrt(3.) - 2./81.*PI2 + 8./27.*Xb
326  - a(1.)/12. + 7.*b(1.)/6. - 2.*b(z) + 26./243.*iPI;
327  case 4:
328  return 2680./2187. + 8.*M_PI/27./sqrt(3.) - 8./243.*PI2 + 32./81.*Xb
329  - a(1.)/9. + 2.*b(1.)/9. + 56./729.*iPI;
330  case 5:
331  return 78301./729. + 8.*M_PI/9./sqrt(3.) - 40./81.*PI2 + 32./27.*Xb
332  - 13.*a(1.)/3. + 38.*b(1.)/3. - 12.*a(z) - 20.*b(z) + 3908./243.*iPI;
333  case 6:
334  return 62440./2187. + 32.*M_PI/27./sqrt(3.) - 160./243.*PI2 + 128./81.*Xb
335  - 16.*a(1.)/9. + 32.*b(1.)/9. + 896./729.*iPI;
336  case 7:
337  return -25./27. - 2./9.*iPI;
338  default:
339  std::stringstream out;
340  out << i;
341  throw std::runtime_error("Bsgamma::r1_ew(): index " + out.str() + " not implemented");
342  }
343 }
344 
345 gslpp::complex Bsgamma::Gamma_t(double t)
346 {
347  if (t<4) return -2. * atan( sqrt(t/(4.-t)) ) * atan( sqrt(t/(4.-t)) );
348  else return -M_PI*M_PI/2. + 2.*log( ( sqrt(t) + sqrt(t-4.) ) / 2. )*log( ( sqrt(t) + sqrt(t-4.) ) / 2. )
349  - 2.*gslpp::complex::i()*M_PI*log( ( sqrt(t) + sqrt(t-4.) ) / 2. );
350 }
351 
352 
353 gslpp::complex Bsgamma::kappa(double Mq,double t)
354 {
355  double s = t * Mb_kin*Mb_kin/Mq/Mq;
356  return 1./2. + Gamma_t(s)/s;
357 }
358 
359 double Bsgamma::Int_b1(double E0)
360 {
361  if (Intb1Cached == 0) {
362  double t1 = (1. - delta(E0));
363 
364  INT = convertToGslFunction(boost::bind(&Bsgamma::getKb_re_1mt, &(*this), _1));
365  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
366  double mt = avaINT;
367 
368  INT = convertToGslFunction(boost::bind(&Bsgamma::getKb_re_1mt2, &(*this), _1));
369  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
370  double mt2 = avaINT;
371 
372  CacheIntb1 = delta(E0)*mt + mt2;
373  Intb1Cached = 1;
374  }
375 
376  return CacheIntb1;
377 }
378 
379 double Bsgamma::Int_b2(double E0)
380 {
381  if (Intb2Cached == 0) {
382  double t1 = (1. - delta(E0));
383 
384  INT = convertToGslFunction(boost::bind(&Bsgamma::getKb_re_t_1mt, &(*this), _1));
385  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
386  double mt = avaINT;
387 
388  INT = convertToGslFunction(boost::bind(&Bsgamma::getKb_re_t_1mt2, &(*this), _1));
389  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
390  double mt2 = avaINT;
391 
392  CacheIntb2 = delta(E0)*mt + mt2;
393  Intb2Cached = 1;
394  }
395 
396  return CacheIntb2;
397 }
398 
399 double Bsgamma::Int_b3(double E0)
400 {
401  if (Intb3Cached == 0) {
402  double t1 = (1. - delta(E0));
403 
404  INT = convertToGslFunction(boost::bind(&Bsgamma::getKb_re_t, &(*this), _1));
405  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
406  double t = avaINT;
407 
408  INT = convertToGslFunction(boost::bind(&Bsgamma::getKb_re_t_1mt, &(*this), _1));
409  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
410  double mt = avaINT;
411 
412  CacheIntb3 = delta(E0)*t + mt;
413  Intb3Cached = 1;
414  }
415 
416  return CacheIntb3;
417 }
418 
419 double Bsgamma::Int_b4(double E0)
420 {
421  if (Intb4Cached == 0) {
422  double t1 = (1. - delta(E0));
423 
424  INT = convertToGslFunction(boost::bind(&Bsgamma::getKb_re_t2_1mt, &(*this), _1));
425  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
426  double mt = avaINT;
427 
428  INT = convertToGslFunction(boost::bind(&Bsgamma::getKb_re_t2_1mt2, &(*this), _1));
429  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
430  double mt2 = avaINT;
431 
432  CacheIntb4 = delta(E0)*mt + mt2;
433  Intb4Cached = 1;
434  }
435 
436  return CacheIntb4;
437 }
438 
439 double Bsgamma::Int_bb1(double E0)
440 {
441  if (Intbb1Cached == 0) {
442  double t1 = (1. - delta(E0));
443 
444  INT = convertToGslFunction(boost::bind(&Bsgamma::getKb_abs2_1mt, &(*this), _1));
445  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
446  double mt = avaINT;
447 
448  INT = convertToGslFunction(boost::bind(&Bsgamma::getKb_abs2_1mt2, &(*this), _1));
449  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
450  double mt2 = avaINT;
451 
452  CacheIntbb1 = delta(E0)*mt + mt2;
453  Intbb1Cached = 1;
454  }
455 
456  return CacheIntbb1;
457 }
458 
459 double Bsgamma::Int_bb2(double E0)
460 {
461  if (Intbb2Cached == 0) {
462  double t1 = (1. - delta(E0));
463 
464  INT = convertToGslFunction(boost::bind(&Bsgamma::getKb_abs2_t_1mt, &(*this), _1));
465  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
466  double mt = avaINT;
467 
468  INT = convertToGslFunction(boost::bind(&Bsgamma::getKb_abs2_t_1mt2, &(*this), _1));
469  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
470  double mt2 = avaINT;
471 
472  CacheIntbb2 = delta(E0)*mt + mt2;
473  Intbb2Cached = 1;
474  }
475 
476  return CacheIntbb2;
477 }
478 
479 double Bsgamma::Int_bb4(double E0)
480 {
481  if (Intbb4Cached == 0) {
482  double t1 = (1. - delta(E0));
483 
484  INT = convertToGslFunction(boost::bind(&Bsgamma::getKb_abs2_t2_1mt, &(*this), _1));
485  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
486  double mt = avaINT;
487 
488  INT = convertToGslFunction(boost::bind(&Bsgamma::getKb_abs2_t2_1mt2, &(*this), _1));
489  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
490  double mt2 = avaINT;
491 
492  CacheIntbb4 = delta(E0)*mt + mt2;
493  Intbb4Cached = 1;
494  }
495 
496  return CacheIntbb4;
497 }
498 
499 gslpp::complex Bsgamma::Int_bc1(double E0)
500 {
501  if (Intbc1Cached == 0) {
502  double t1 = (1. - delta(E0));
503 
504  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_re_Kb_1mt, &(*this), _1));
505  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
506  gslpp::complex mt = avaINT;
507 
508  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_im_Kb_1mt, &(*this), _1));
509  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
510  mt += gslpp::complex::i() * avaINT;
511 
512  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_re_Kb_1mt2, &(*this), _1));
513  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
514  gslpp::complex mt2 = avaINT;
515 
516  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_im_Kb_1mt2, &(*this), _1));
517  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
518  mt2 += gslpp::complex::i() * avaINT;
519 
520  CacheIntbc1 = delta(E0)*mt + mt2;
521  Intbc1Cached = 1;
522  }
523 
524  return CacheIntbc1;
525 }
526 
527 gslpp::complex Bsgamma::Int_bc2(double E0)
528 {
529  if (Intbc2Cached == 0) {
530  double t1 = (1. - delta(E0));
531 
532  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_re_Kb_t_1mt, &(*this), _1));
533  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
534  gslpp::complex mt = avaINT;
535 
536  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_im_Kb_t_1mt, &(*this), _1));
537  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
538  mt += gslpp::complex::i() * avaINT;
539 
540  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_re_Kb_t_1mt2, &(*this), _1));
541  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
542  gslpp::complex mt2 = avaINT;
543 
544  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_im_Kb_t_1mt2, &(*this), _1));
545  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
546  mt2 += gslpp::complex::i() * avaINT;
547 
548  CacheIntbc2 = delta(E0)*mt + mt2;
549  Intbc2Cached = 1;
550  }
551 
552  return CacheIntbc2;
553 }
554 
555 gslpp::complex Bsgamma::Int_c1(double E0)
556 {
557  if (Intc1Cached == 0) {
558  double t1 = (1. - delta(E0));
559 
560  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_re_1mt, &(*this), _1));
561  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
562  gslpp::complex mt = avaINT;
563 
564  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_im_1mt, &(*this), _1));
565  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
566  mt += gslpp::complex::i() * avaINT;
567 
568  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_re_1mt2, &(*this), _1));
569  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
570  gslpp::complex mt2 = avaINT;
571 
572  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_im_1mt2, &(*this), _1));
573  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
574  mt2 += gslpp::complex::i() * avaINT;
575 
576  CacheIntc1 = delta(E0)*mt + mt2;
577  Intc1Cached = 1;
578  }
579 
580  return CacheIntc1;
581 }
582 
583 gslpp::complex Bsgamma::Int_c2(double E0)
584 {
585  if (Intc2Cached == 0) {
586  double t1 = (1. - delta(E0));
587 
588  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_re_t_1mt, &(*this), _1));
589  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
590  gslpp::complex mt = avaINT;
591 
592  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_im_t_1mt, &(*this), _1));
593  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
594  mt += gslpp::complex::i() * avaINT;
595 
596  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_re_t_1mt2, &(*this), _1));
597  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
598  gslpp::complex mt2 = avaINT;
599 
600  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_im_t_1mt2, &(*this), _1));
601  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
602  mt2 += gslpp::complex::i() * avaINT;
603 
604  CacheIntc2 = delta(E0)*mt + mt2;
605  Intc2Cached = 1;
606  }
607 
608  return CacheIntc2;
609 }
610 
611 gslpp::complex Bsgamma::Int_c3(double E0)
612 {
613  if (Intc3Cached == 0) {
614  double t1 = (1. - delta(E0));
615 
616  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_re_t, &(*this), _1));
617  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
618  gslpp::complex t = avaINT;
619 
620  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_im_t, &(*this), _1));
621  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
622  t += gslpp::complex::i() * avaINT;
623 
624  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_re_t_1mt, &(*this), _1));
625  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
626  gslpp::complex mt = avaINT;
627 
628  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_im_t_1mt, &(*this), _1));
629  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
630  mt += gslpp::complex::i() * avaINT;
631 
632  CacheIntc3 = delta(E0)*t + mt;
633  Intc3Cached = 1;
634  }
635 
636  return CacheIntc3;
637 }
638 
639 double Bsgamma::Int_cc(double E0)
640 {
641  if (IntccCached == 0) {
642  double t1 = (1. - delta(E0));
643 
644  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_abs2_t, &(*this), _1));
645  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
646  double mt = avaINT;
647 
648  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_abs2_t_1mt, &(*this), _1));
649  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
650  double mt2 = avaINT;
651 
652  CacheIntcc = delta(E0)*mt + 2. * mt2;
653  IntccCached = 1;
654  }
655 
656  return CacheIntcc;
657 }
658 
659 double Bsgamma::Int_cc1(double E0)
660 {
661  if (Intcc1Cached == 0) {
662  double t1 = (1. - delta(E0));
663 
664  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_abs2_1mt, &(*this), _1));
665 
666  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
667  double mt = avaINT;
668 
669  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_abs2_1mt2, &(*this), _1));
670  if (gsl_integration_cquad(&INT, t1, 1., 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
671  double mt2 = avaINT;
672 
673  CacheIntcc1 = delta(E0)*mt + mt2;
674  Intcc1Cached = 1;
675  }
676 
677  return CacheIntcc1;
678 }
679 
680 double Bsgamma::Int_cc1_part1(double E0)
681 {
682  if (Intcc1p1Cached == 0) {
683  double t1 = (1. - delta(E0));
684 
685  INT = convertToGslFunction(boost::bind(&Bsgamma::getKc_abs2_1mt, &(*this), _1));
686  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
687 
688  CacheIntcc1p1 = delta(E0)*avaINT;
689  Intcc1p1Cached = 1;
690  }
691 
692  return CacheIntcc1p1;
693 }
694 
695 double Bsgamma::ff7_dMP(double E0)
696 {
697  if (FOUR_BODY){
698  double d=delta(E0);
699  double d2=d*d;
700  double d3=d2*d;
701 
702  return 4. * d * (18. - 33.*d + 2.*d2 + 13.*d3 - 6.* d2 * (2. + d) * log(d))
703  / (81. * (d - 1.));
704  }
705  else
706  return 0.;
707 }
708 
709 double Bsgamma::ff7_sMP(double E0)
710 {
711  if (FOUR_BODY){
712  double d=delta(E0);
713  double d2=d*d;
714  double d3=d2*d;
715 
716  return (-2. * d * (72. + 39.*d - 76.*d2 - 35.*d3
717  + 6.*d*(18. + 13.*d + 2.*d2)*log(d))) / (243.*(d - 1.));
718  }
719  else
720  return 0.;
721 }
722 
723 double Bsgamma::ff8_dMP(double E0)
724 {
725  if (FOUR_BODY){
726  double d=delta(E0);
727  double d2=d*d;
728  double d3=d2*d;
729  double ld = log(d);
730  double l1d = log(1. - d);
731  double Li2 = gsl_sf_dilog(d);
732 
733  return -136./27. * d - 724./81. * d2 + 20./27. * d3
734  + (-8./9. + 16./9. * d - 8./9. * d2) * l1d* l1d
735  + (32./27. * d + 76./27. * d2 - 16./81. * d3) * ld
736  + (-104./27. - 80./9. * d + 40./9. * d2 + (32./27.
737  + 32./9. * d - 16./9. * d2) * ld) * l1d
738  + (-64./27. * d - 152./27. * d2 + 32./81. * d3
739  + (-64./27. - 64./9. * d + 32./9. * d2) * l1d) * log(Ms/Mb_kin)
740  + (32./27. + 32./9. * d - 16./9. * d2) * Li2;
741  }
742  else
743  return 0.;
744 }
745 
746 double Bsgamma::ff8_sMP(double E0)
747 {
748  if (FOUR_BODY){
749  double d=delta(E0);
750  double d2=d*d;
751  double d3=d2*d;
752  double ld = log(d);
753  double l1d = log(1. - d);
754  double Li2 = gsl_sf_dilog(d);
755 
756  return -340./243. * d - 104./81. * d2 + 16./729. * d3
757  + (-4./27. + 8./27. * d - 4./27. * d2) * l1d* l1d
758  + (8./27. * d + 4./9. * d2) * ld
759  + (-16./27. * d - 8./9. * d2) * log(Ms/Mb_kin)
760  + (-268./243. - 40./27. * d + 20./27. * d2 + (8./27.
761  + 16./27. * d - 8./27. * d2) * ld
762  + (-16./27. - 32./27. * d + 16./27. * d2) * log(Ms/Mb_kin)) * l1d
763  + (8./27. + 16./27. * d - 8./27. * d2) * Li2;
764  }
765  else
766  return 0.;
767 }
768 
769 double Bsgamma::Phi11_1(double E0)
770 {
771  return Phi22_1(E0)/36.;
772 }
773 
774 double Bsgamma::Phi12_1(double E0)
775 {
776  return -Phi22_1(E0)/3.;
777 }
778 
779 gslpp::complex Bsgamma::Phi13_1(double E0)
780 {
781  return -Phi23_1(E0)/6.;
782 }
783 
784 gslpp::complex Bsgamma::Phi14_1(double E0)
785 {
786  return -Phi24_1(E0)/6.;
787 }
788 
789 gslpp::complex Bsgamma::Phi15_1(double E0)
790 {
791  return -Phi25_1(E0)/6.;
792 }
793 
794 gslpp::complex Bsgamma::Phi16_1(double E0)
795 {
796  return -Phi26_1(E0)/6.;
797 }
798 
799 gslpp::complex Bsgamma::Phi17_1(double E0, double z)
800 {
801  return -Phi27_1(E0,z)/6.;
802 }
803 
804 gslpp::complex Bsgamma::Phi18_1(double E0, double z)
805 {
806  return Phi27_1(E0,z)/18.;
807 }
808 
809 double Bsgamma::Phi22_1(double E0)
810 {
811  return 16./27. * Int_cc1(E0);
812 }
813 
814 double Bsgamma::Phi23_1_4body(double E0)
815 {
816  if (FOUR_BODY)
817  return 0.0039849625073434735;
818  else
819  return 0.;
820 }
821 
822 gslpp::complex Bsgamma::Phi23_1(double E0)
823 {
824  return -8./27. * (Int_c1(E0) + Int_c2(E0) + 2.*Int_bc1(E0) - 2.*Int_bc2(E0))
825  - Phi23_1_4body(E0);
826 }
827 
828 double Bsgamma::Phi24_1_4body(double E0)
829 {
830  if (FOUR_BODY)
831  return 0.012330977673588935;
832  else
833  return 0.;
834 }
835 
836 gslpp::complex Bsgamma::Phi24_1(double E0)
837 {
838  return -1./6. * (Phi23_1(E0) + Phi23_1_4body(E0))
839  - Phi24_1_4body(E0);
840 }
841 
842 double Bsgamma::Phi25_1_4body(double E0)
843 {
844  if (FOUR_BODY)
845  return 0.06375940011749558;
846  else
847  return 0.;
848 }
849 
850 gslpp::complex Bsgamma::Phi25_1(double E0)
851 {
852  return -32./27. * (4.*Int_c1(E0) + Int_c2(E0) + 8.*Int_bc1(E0) - 2.*Int_bc2(E0))
853  - Phi25_1_4body(E0);
854 }
855 
856 double Bsgamma::Phi26_1_4body(double E0)
857 {
858  if (FOUR_BODY)
859  return 0.11932481422855279;
860  else
861  return 0.;
862 }
863 
864 gslpp::complex Bsgamma::Phi26_1(double E0)
865 {
866  return 16./81. * (4.*Int_c1(E0) + Int_c2(E0) - 10.*Int_bc1(E0) - 2.*Int_bc2(E0) + 36.*Int_cc1(E0))
867  - Phi26_1_4body(E0);
868 }
869 
870 gslpp::complex Bsgamma::Phi27_1(double E0, double z)
871 {
872  double d = delta(E0);
873  double d2 = d*d;
874  double z2 = z*z;
875  double Pi2 = M_PI*M_PI;
876  double st0 = sqrt(1. - 4.*z);
877  double std = sqrt( (1. - d - 4.*z) * (1. - d) );
878  double L0 = log( ( 1. + st0 ) / ( 2.*sqrt(z) ) );
879  double Ld = log( ( sqrt(1. - d) + sqrt(1. - d - 4.*z) ) / ( 2.*sqrt(z) ) );
880 
881  gslpp::complex res;
882 
883  if (d == 1) {
884  res = -2./27. + (2.*Pi2 - 7.)/9. * z + 4.*(3. - 2.*Pi2)/9. * z * z
885  + 4./3. * z * (1. - 2.*z) * st0 * L0
886  - 8./9. * z * (6.*z*z - 4.*z + 1.) * L0*L0 + 4./3. * Pi2 * z * z *z;
887  } else res = -2./27. * d * (3. - 3.*d + d2) + (2.*Pi2 - 7.)/9. * z * d * (2. - d)
888  + 4.*(3. - 2.*Pi2)/9. * z * z * d
889  + 4./3. * z * (1. - 2.*z) * ( st0 * L0 - std * Ld )
890  + 4./3. * z * d * std * Ld
891  - 8./9. * z * (6.*z*z - 4.*z + 1.) * ( L0*L0 - Ld*Ld )
892  - 8./9. * z * d * (2. - d - 4.*z) * Ld * Ld;
893 
894  if (z < (1. - d)/4.)
895  res += gslpp::complex::i() * 8./9. * M_PI * z * ( (1. - 4. * z + 6. * z2)* (L0-Ld) - 3./4. * (1. - 2. * z) * (st0-std)
896  + d * (2. - d - 4. * z) * Ld - 3./4. * d * std );
897  else
898  res += gslpp::complex::i() * 8./9. * M_PI * z * ( (1. - 4. * z + 6. * z2) * L0 - 3./4. * (1. - 2. * z) * st0 );
899 
900  return res;
901 }
902 
903 gslpp::complex Bsgamma::Phi28_1(double E0, double z)
904 {
905  return -Phi27_1(E0, z)/3.;
906 }
907 
908 double Bsgamma::Phi33_1(double E0)
909 {
910  double d=delta(E0);
911  double d2=d*d;
912  double d3=d2*d;
913  double d4=d3*d;
914 
915  return 2./27. * (Int_b1(E0) + 8.*Int_b2(E0) - 4.*Int_b4(E0)
916  + 33.*Int_bb1(E0) - 20.*Int_bb2(E0) + 4.*Int_bb4(E0))
917  + 1./18. * d * ( 1./2. - 1./2.*d2 + 1./3.*d3 - 1./15.*d4 );
918 }
919 
920 double Bsgamma::Phi34_1(double E0)
921 {
922  return -1./3.*Phi33_1(E0);
923 }
924 
925 double Bsgamma::Phi35_1(double E0)
926 {
927  double d=delta(E0);
928  double d2=d*d;
929  double d3=d2*d;
930  double d4=d3*d;
931 
932  return 32./27. * (2.*Int_b1(E0) + 4.*Int_b2(E0) - 2.*Int_b4(E0)
933  + 18.*Int_bb1(E0) - 13.*Int_bb2(E0) + 2.*Int_bb4(E0))
934  + 4./9. * d * ( 4./3. - d2 + 1./2.*d3 - 1./15.*d4 );
935 }
936 
937 gslpp::complex Bsgamma::Phi36_1(double E0)
938 {
939  double d=delta(E0);
940  double d2=d*d;
941  double d3=d2*d;
942  double d4=d3*d;
943 
944  return 8./81. * (5.*Int_b1(E0) + Int_b2(E0) + 4.*Int_b4(E0)
945  - 18.*Int_bb1(E0) + 8.*Int_bb2(E0) - 4.*Int_bb4(E0))
946  + 6. * (Phi23_1(E0) + Phi23_1_4body(E0))
947  - 2./27. * d * ( 4./3. - d2 + 1./2.*d3 - 1./15.*d4 );
948 }
949 
950 double Bsgamma::Phi37_1(double E0)
951 {
952  double d=delta(E0);
953  double d2=d*d;
954 
955  return -4./3. * Int_b3(E0) + 1./9. * d * (1. - d + 1./3.*d2) + 1./4.*ff7_sMP(E0);
956 }
957 
958 double Bsgamma::Phi38_1(double E0)
959 {
960  return -1./3. * ( Phi37_1(E0) - 1./4.*ff7_sMP(E0) ) + 1./4.*ff8_sMP(E0);
961 }
962 
963 double Bsgamma::Phi44_1(double E0)
964 {
965  return 1./36. * Phi33_1(E0);
966 }
967 
968 double Bsgamma::Phi45_1(double E0)
969 {
970  return -1./6. * Phi35_1(E0);
971 }
972 
973 gslpp::complex Bsgamma::Phi46_1(double E0)
974 {
975  return -1./6. * Phi36_1(E0);
976 }
977 
978 double Bsgamma::Phi47_1(double E0)
979 {
980  return -1./6. * ( Phi37_1(E0) - 1./4.*ff7_sMP(E0) )
981  + 1./4. * (-1./6. * ff7_sMP(E0) + ff7_dMP(E0));
982 }
983 
984 double Bsgamma::Phi48_1(double E0)
985 {
986  return -1./3. * (Phi47_1(E0) - 1./4. * (-1./6. * ff7_sMP(E0) + ff7_dMP(E0)) )
987  + 1./4. * (-1./6. * ff8_sMP(E0) + ff8_dMP(E0));
988 }
989 
990 double Bsgamma::Phi55_1(double E0)
991 {
992  double d=delta(E0);
993  double d2=d*d;
994  double d3=d2*d;
995  double d4=d3*d;
996 
997  return 128./27. * (4.*Int_b1(E0) + 2.*Int_b2(E0) - Int_b4(E0)
998  + 12.*Int_bb1(E0) - 8.*Int_bb2(E0) + Int_bb4(E0))
999  + 8./9. * d * ( 11./3. - 2.*d2 + 2./3.*d3 - 1./15.*d4 );
1000 }
1001 
1002 gslpp::complex Bsgamma::Phi56_1(double E0)
1003 {
1004  double d=delta(E0);
1005  double d2=d*d;
1006  double d3=d2*d;
1007  double d4=d3*d;
1008 
1009  return 32./81. * (20.*Int_b1(E0) + Int_b2(E0) + 4.*Int_b4(E0)
1010  + 24.*Int_bb1(E0) + 14.*Int_bb2(E0) - 4.*Int_bb4(E0))
1011  + 6. * (Phi25_1(E0) + Phi25_1_4body(E0))
1012  - 8./27. * d * ( 11./3. - 2.*d2 + 2./3.*d3 - 1./15.*d4 );
1013 }
1014 
1015 double Bsgamma::Phi57_1(double E0)
1016 {
1017  double d=delta(E0);
1018  double d2=d*d;
1019 
1020  return 16./9. * d * ( 1. - d + 1./3.*d2) + 4. * ff7_sMP(E0);
1021 }
1022 
1023 double Bsgamma::Phi58_1(double E0)
1024 {
1025  return -1./3. * (Phi57_1(E0) - 4. * ff7_sMP(E0))
1026  + 4. * ff8_sMP(E0);
1027 }
1028 
1029 gslpp::complex Bsgamma::Phi66_1(double E0)
1030 {
1031  double d=delta(E0);
1032  double d2=d*d;
1033  double d3=d2*d;
1034  double d4=d3*d;
1035 
1036  return -8./243. * (56.*Int_b1(E0) + 10.*Int_b2(E0) + 4.*Int_b4(E0)
1037  + 15.*Int_bb1(E0) - 4.*Int_bb2(E0) - 4.*Int_bb4(E0))
1038  + 32./27. * (4.*Int_c1(E0) + Int_c2(E0) - Int_bc1(E0)
1039  - 2.*Int_bc2(E0) + 9.*Int_cc1(E0))
1040  + 2./81. * d * ( 11./3. - 2.*d2 + 2./3.*d3 - 1./15.*d4 );
1041 }
1042 
1043 gslpp::complex Bsgamma::Phi67_1(double E0)
1044 {
1045  double d=delta(E0);
1046  double d2=d*d;
1047 
1048  return 8./3. * (Int_b3(E0) - 2.*Int_c3(E0)) - 8./27. * d * ( 1. - d + 1./3.*d2)
1049  + 1./4. * (-8./3. * ff7_sMP(E0) + 10. * ff7_dMP(E0));
1050 }
1051 
1052 gslpp::complex Bsgamma::Phi68_1(double E0)
1053 {
1054  return -1./3. * (Phi67_1(E0) - 1./4. * (-8./3. * ff7_sMP(E0) + 10. * ff7_dMP(E0)) )
1055  + 1./4. * (-8./3. * ff8_sMP(E0) + 10. * ff8_dMP(E0));
1056 }
1057 
1058 double Bsgamma::Phi77_1(double E0)
1059 {
1060  double d=delta(E0);
1061  double d2=d*d;
1062  double d3=d2*d;
1063 
1064  return -2./3.*pow(log(d),2.) - 7./3.*log(d) - 31./9. + 10./3.*d + d2/3. - 2./9.*d3 + d*(d - 4.)*log(d)/3.;
1065 }
1066 
1067 double Bsgamma::Phi78_1(double E0)
1068 {
1069  double d=delta(E0);
1070  double d2=d*d;
1071  double d3=d2*d;
1072 
1073  double Li2 = gsl_sf_dilog(1. - d);
1074 
1075  double pi2=M_PI*M_PI;
1076 
1077  return 8./9.*( Li2 - pi2/6. - d*log(d) + 9./4.*d - d2/4. + d3/12.);
1078 }
1079 
1080 double Bsgamma::Phi88_1(double E0)
1081 {
1082  double d=delta(E0);
1083  double d2=d*d;
1084  double d3=d2*d;
1085 
1086  double Li2 = gsl_sf_dilog(1. - d);
1087 
1088  double pi2=M_PI*M_PI;
1089 
1090  return 1./27.*( -2.*log(Mb_kin/Ms)*( d2 + 2.*d + 4.*log(1. - d) ) + 4.*Li2 - 2./3.*pi2 - d*(2. + d)*log(d)
1091  + 8.*log(1. - d) - 2./3.*d3 + 3.*d2 +7*d);
1092 }
1093 
1094 gslpp::complex Bsgamma::Kij_1(int i, int j, double E0, double mu)
1095 {
1096  if (i > 8 || j>8 || i<1 || j<1) throw std::runtime_error("Bsgamma::Kij_1(): indices (i,j) must be included in (1,..,8)");
1097 
1098  double gamma_i7[8] = {-208./243., 416./81., -176./81., -152./243., -6272./81., 4624./243., 32./3., -32./9.};
1099  gslpp::complex K_ij[8][8] = {{0.}};
1100  double Lb = log(mu/Mb_kin);
1101 
1102  K_ij[0][0] = 4.*Phi11_1(E0);
1103  K_ij[0][1] = 2.*Phi12_1(E0);
1104  K_ij[0][2] = 2.*Phi13_1(E0);
1105  K_ij[0][3] = 2.*Phi14_1(E0);
1106  K_ij[0][4] = 2.*Phi15_1(E0);
1107  K_ij[0][5] = 2.*Phi16_1(E0);
1108  K_ij[0][6] = r1(1,zeta()) - gamma_i7[0]*Lb + 2.*Phi17_1(E0, zeta());
1109  K_ij[0][7] = 2.*Phi18_1(E0, zeta());
1110 
1111  K_ij[1][1] = 4.*Phi22_1(E0);
1112  K_ij[1][2] = 2.*Phi23_1(E0);
1113  K_ij[1][3] = 2.*Phi24_1(E0);
1114  K_ij[1][4] = 2.*Phi25_1(E0);
1115  K_ij[1][5] = 2.*Phi26_1(E0);
1116  K_ij[1][6] = r1(2,zeta()) - gamma_i7[1]*Lb + 2.*Phi27_1(E0, zeta());
1117  K_ij[1][7] = 2.*Phi28_1(E0, zeta());
1118 
1119  K_ij[2][2] = 4.*Phi33_1(E0);
1120  K_ij[2][3] = 2.*Phi34_1(E0);
1121  K_ij[2][4] = 2.*Phi35_1(E0);
1122  K_ij[2][5] = 2.*Phi36_1(E0);
1123  K_ij[2][6] = r1(3,zeta()) - gamma_i7[2]*Lb + 2.*Phi37_1(E0);
1124  K_ij[2][7] = 2.*Phi38_1(E0);
1125 
1126  K_ij[3][3] = 4.*Phi44_1(E0);
1127  K_ij[3][4] = 2.*Phi45_1(E0);
1128  K_ij[3][5] = 2.*Phi46_1(E0);
1129  K_ij[3][6] = r1(4,zeta()) - gamma_i7[3]*Lb + 2.*Phi47_1(E0);
1130  K_ij[3][7] = 2.*Phi48_1(E0);
1131 
1132  K_ij[4][4] = 4.*Phi55_1(E0);
1133  K_ij[4][5] = 2.*Phi56_1(E0);
1134  K_ij[4][6] = r1(5,zeta()) - gamma_i7[4]*Lb + 2.*Phi57_1(E0);
1135  K_ij[4][7] = 2.*Phi58_1(E0);
1136 
1137  K_ij[5][5] = 4.*Phi66_1(E0);
1138  K_ij[5][6] = r1(6,zeta()) - gamma_i7[5]*Lb + 2.*Phi67_1(E0);
1139  K_ij[5][7] = 2.*Phi68_1(E0);
1140 
1141  K_ij[6][6] = -182./9. + 8./9.*M_PI*M_PI - gamma_i7[6]*2.*Lb + 4.*Phi77_1(E0);
1142  K_ij[6][7] = r1(8,zeta()) - gamma_i7[7]*Lb + 2.*Phi78_1(E0);
1143 
1144  K_ij[7][7] = 4.*Phi88_1(E0);
1145 
1146  if (j >= i ) return K_ij[i-1][j-1];
1147  else return K_ij[j-1][i-1].conjugate();
1148 }
1149 
1150 double Bsgamma::Rer22(double z)
1151 {
1152  double L = log(z);
1153  double L2 = L*L;
1154  double L3 = L2*L;
1155  double L4 = L3*L;
1156  double z32 = sqrt(z)*z;
1157  double z2 = z*z;
1158  double z52 = z32*z;
1159  double z3 = z2*z;
1160  double z72 = z52*z;
1161  double z4 = z3*z;
1162  double Pi2 = M_PI*M_PI;
1163  double zeta3 = gsl_sf_zeta_int(3);
1164 
1165  return 67454./6561. - 124./729. * Pi2
1166  - 4./1215. * (11280. - 1520. * Pi2 - 171. * Pi2*Pi2 - 5760. * zeta3
1167  + 6840. * L - 1440. * Pi2*L - 2520. * zeta3*L
1168  + 120. * L2 + 100. * L3 - 30. * L4) * z
1169  - 64./243. * Pi2*( 43. - 12. * log(2.) - 3. * L) * z32
1170  - 2./1215. * (11475. - 380. * Pi2 + 96. * Pi2*Pi2
1171  + 7200. * zeta3 - 1110. * L - 1560. * Pi2*L + 1440. * zeta3*L
1172  + 990. * L2 + 260. * L3 - 60. * L4) * z2
1173  + 2240./243. * Pi2 * z52
1174  - 2./2187. * (62471. - 2424. * Pi2 - 33264. * zeta3 - 19494. * L
1175  - 504. * Pi2*L - 5184. * L2 + 2160. * L3) * z3
1176  - 2464./6075. * Pi2 * z72
1177  + ( - 15103841./546750. + 7912./3645. * Pi2 + 2368./81. * zeta3
1178  + 147038./6075. * L + 352./243. * Pi2*L + 88./243. * L2
1179  - 512./243. * L3 ) * z4;
1180 }
1181 
1182 double Bsgamma::Phi22_2beta0(double E0, double mu)
1183 {
1184  double Lb = 2*log(mu/Mb_kin);
1185  double d = delta(E0);
1186  double d2 = d*d;
1187  double mcmb = Mc/Mb_kin;
1188  double mcmb2 = mcmb*mcmb;
1189  double mcmb3 = mcmb2*mcmb;
1190 
1191  return SM.Beta0(5) * (Phi22_1(E0)*Lb
1192  + 0.013698269459646965 + 0.3356948452887703 * d
1193  - 0.086677232161681 * d2
1194  + ( 0.3575455009710223 + 1.8248223618702617 * d
1195  - 0.374324331239819 * d2 ) * mcmb
1196  + (-2.3059130759599302 - 5.799640881350228 * d
1197  - 6.226247001127346 * d2 ) * mcmb2
1198  + ( 3.4485885608332834 - 0.5479757965141787 * d
1199  + 17.272487170738795 * d2 ) * mcmb3);
1200 }
1201 
1202 double Bsgamma::Phi28_2beta0(double E0, double mu)
1203 {
1204  double Lb = 2*log(mu/Mb_kin);
1205  double d = delta(E0);
1206  double d2 = d*d;
1207  double mcmb = Mc/Mb_kin;
1208  double mcmb2 = mcmb*mcmb;
1209  double mcmb3 = mcmb2*mcmb;
1210  double mcmb4 = mcmb3*mcmb;
1211  double mcmb5 = mcmb4*mcmb;
1212 
1213  return SM.Beta0(5) * (Phi28_1(E0, zeta()).real()*Lb
1214  + 0.026054745293391798 + 0.1678721564514209 * d
1215  - 0.19700988587274693 * d2
1216  + ( -0.03801105485376407 + 0.601712887338462 * d
1217  - 0.7557529126506585 * d2 ) * mcmb
1218  + ( 2.7551159092192132 - 10.034450524236696 * d
1219  + 11.271837772655209 * d2 ) * mcmb2
1220  + ( -27.045289848315868 + 68.46851531490181 * d
1221  - 72.50921751760909 * d2 ) * mcmb3
1222  + ( 85.86574743951778 - 289.3441408351491 * d
1223  + 297.6777008484198 * d2 ) * mcmb4
1224  + ( -91.5260435658921 + 399.81982774456964 * d
1225  - 399.85440571662446 * d2 ) * mcmb5);
1226 }
1227 
1228 double Bsgamma::Phi77_2beta0(double E0, double mu)
1229 {
1230  double Lb = 2*log(mu/Mb_kin);
1231  double d = delta(E0);
1232  double d2 = d*d;
1233  double d3 = d2*d;
1234  double Ld = log(d);
1235  double zeta3 = gsl_sf_zeta_int(3);
1236  double Li2 = gsl_sf_dilog(1. - d);
1237  Polylogarithms Poly;
1238  double Li3 = Poly.Li3(d);
1239 
1240  return SM.Beta0(5) * (Phi77_1(E0)*Lb
1241  + ( -3. + 4./3. * d - 1./3. * d2 - 4./3. * Ld ) * Li2
1242  + ( 13./18. + 2. * d - 1./2. * d2 - 4./3. * log(1. - d)
1243  + 2./3. * Ld ) * Ld*Ld
1244  - 8./3. * (Li3 - zeta3)
1245  + ( 4./9. * M_PI*M_PI - 85./18. - 47./9. * d
1246  + 19./18. * d2 + 2./9. * d3 ) * Ld
1247  - 49./6. + 80./9. * d + 1./18. * d2 - 7./9. * d3);
1248 }
1249 
1250 double Bsgamma::Phi88_2beta0(double E0, double mu)
1251 {
1252  double Lb = 2*log(mu/Mb_kin);
1253  double d = delta(E0);
1254  double d2 = d*d;
1255  double d3 = d2*d;
1256  double Ld = log(d);
1257  double L1d = log(1. - d);
1258  double Li2 = gsl_sf_dilog(1. - d);
1259  Polylogarithms Poly;
1260  double Li3 = Poly.Li3(d);
1261 
1262  return SM.Beta0(5) * (Phi88_1(E0)*Lb
1263  + 4./27. * ( - 2. * ( Li2 - 1./6. * M_PI*M_PI + 3. * L1d
1264  - 1./4. * d * (2. + d) * Ld + 8./3. * d + 5./6. * d2
1265  - 1./18. * d3 ) * log(Mb_kin/Ms) - 2. * Li3
1266  + ( 5. - 2. * Ld ) * ( Li2 - 1./6. * M_PI*M_PI )
1267  + ( 1./2. * d + 1./4. * d2 - L1d ) * Ld * Ld
1268  - 1./12. * M_PI*M_PI * d * (2. + d)
1269  + ( 151./18. - 1./3. * M_PI*M_PI ) * L1d
1270  + ( - 53./12. * d - 19./12. * d2 + 2./9. * d3 ) * Ld
1271  + 787./72. * d + 227./72. * d2 - 41./72. * d3 ));
1272 }
1273 
1274 double Bsgamma::dY1(double E0)
1275 {
1276  double z0 = 1. - delta(E0);
1277  double Li2 = gsl_sf_dilog(z0);
1278 
1279  return + 2./9. * z0*(z0*z0 + 24.)- 8./3. * (z0 - 1.)*log(1. - z0) - 8./3. * Li2;
1280 }
1281 
1282 double Bsgamma::Y1(double E0, double mu)
1283 {
1284  double Lb = log(mu/Mb_kin);
1285 
1286  return 4./9.*(29. - 2.*M_PI*M_PI) + 16./3.*Lb - dY1(E0);
1287 }
1288 
1289 double Bsgamma::Y2CF(double E0, double mu)
1290 {
1291  double Lb = log(mu/Mb_kin);
1292  double z0 = 1. - delta(E0);
1293  double z02 = z0*z0;
1294  double z03 = z02*z0;
1295  double z04 = z03*z0;
1296  double z05 = z04*z0;
1297  double z06 = z05*z0;
1298  double z07 = z06*z0;
1299  double z08 = z07*z0;
1300  double z09 = z08*z0;
1301  double z010 = z09*z0;
1302  double z011 = z010*z0;
1303  double Lz = log(1. - z0);
1304  double Li2 = gsl_sf_dilog(z0);
1305 
1306  return -21.9087 - 112.464 * Lb - 42.6667 * Lb*Lb - 77.7675 * z0 + 10.6667 * Lb*z0
1307  + 68.5848 * z02 - 5.33333 * Lb * z02 - 4.42133 * z03 + 6.22222 * Lb * z03
1308  - 4.0317 * z04 + 6.64376 * z05 - 11.647 * z06 + 15.8697 * z07
1309  - 14.8006 * z08 + 8.85514 * z09 - 2.9929 * z010 + 0.433254 * z011
1310  + (-77.7675 + 85.8251 * z0 - 28.6855 * z02
1311  + Lb * (-21.3333 - 21.3333 * z0 + 5.33333 * z02)) * Lz
1312  + (-12.2881 - 10.6667 * Lb + 6.12213 * z0 + 0.27227 * z02) * Lz*Lz
1313  + (-2.88573 + 5.77146 * z0 - 2.88573 * z02) * Lz*Lz*Lz - 32. * Lb*Li2;
1314 }
1315 
1316 double Bsgamma::Y2CA(double E0, double mu)
1317 {
1318  double Lb = log(mu/Mb_kin);
1319  double z0 = 1. - delta(E0);
1320  double z02 = z0*z0;
1321  double z03 = z02*z0;
1322  double z04 = z03*z0;
1323  double z05 = z04*z0;
1324  double z06 = z05*z0;
1325  double z07 = z06*z0;
1326  double z08 = z07*z0;
1327  double z09 = z08*z0;
1328  double z010 = z09*z0;
1329  double z011 = z010*z0;
1330  double Lz = log(1. - z0);
1331  double Li2 = gsl_sf_dilog(z0);
1332 
1333  return 22.8959 + 76.5729 * Lb + 30.2222 * Lb*Lb + 2.94616 * z0 - 60.4444 * Lb*z0
1334  - 13.2522 * z02 - 6.96907 * z03 - 2.51852 * Lb*z03 + 0.117907 * z04
1335  - 2.02988 * z05 + 2.90402 * z06 - 3.53904 * z07 + 2.55728 * z08
1336  - 0.941549 * z09 + 0.0173599 * z010 + 0.0598012 * z011
1337  + (2.94616 - 30.2222 * Lb- 12.3947 * z0 + 30.2222 * Lb*z0 + 9.44855 * z02) * Lz
1338  - 6.61587 * (1. - z0) * (1. - z0) * Lz*Lz + 30.2222 * Lb*Li2
1339  + (0.69995 - 1.3999 * z0 + 0.69995 * z02) * Lz*Lz*Lz;
1340 }
1341 
1342 double Bsgamma::Y2NL(double E0, double mu)
1343 {
1344  double z0 = 1. - delta(E0);
1345  double Lz = log(1. - z0);
1346  double Lb = log(mu/Mb_kin);
1347  double zeta3 = gsl_sf_zeta_int(3);
1348  double Li2 = gsl_sf_dilog(z0);
1349  Polylogarithms Poly;
1350  double Li3 = Poly.Li3(z0);
1351  double Li3min = Poly.Li3(1. - z0);
1352 
1353  return -16./81. * (328. - 13.*M_PI*M_PI) - 64./27. * (18. - M_PI*M_PI)*Lb
1354  -64./9. * Lb*Lb + 64./3. * zeta3
1355  +4./27. * z0*(7.*z0*z0 - 17.*z0 + 238.) + 8./3. * Lb*dY1(E0)
1356  -8./27. * (z0*z0*z0 - 6.*z0*z0 + 80.*z0 - 75. + 6.*M_PI*M_PI)*Lz
1357  +16./3. * (z0 - 1.)*Lz*Lz + 16./3. * log(z0)*Lz*Lz
1358  +32./27. * (3.*z0 - 8.)*Li2 + 32./3. * Lz*Li2
1359  -32./9. * Li3 + 32./3. * Li3min - 32./3. * zeta3;
1360 }
1361 
1362 double Bsgamma::Y2NV_PHI1(double rho)
1363 {
1364  double y = (1. - sqrt(1. - 4.*rho)) / (1. + sqrt(1. - 4.*rho));
1365 
1366  if (rho < 1./4.)
1367  return log(y)*log(y) - M_PI*M_PI;
1368  else
1369  return - acos( 1. - 1./(2. * rho) ) * acos( 1. - 1./(2. * rho) );
1370 }
1371 
1372 
1373 double Bsgamma::Y2NV_PHI2(double rho)
1374 {
1375  double y = (1. - sqrt(1. - 4.*rho)) / (1. + sqrt(1. - 4.*rho));
1376 
1377  if (rho < 1./4.)
1378  return log(y) * sqrt(1. - 4.*rho);
1379  else
1380  return - acos( 1. - 1./(2. * rho) ) * sqrt(4.*rho - 1.);
1381 }
1382 
1383 
1384 double Bsgamma::Y2NV_PHI3(double rho)
1385 {
1386  double y = (1. - sqrt(1. - 4.*rho)) / (1. + sqrt(1. - 4.*rho));
1387 
1388  if (rho < 1./4.)
1389  return ( gsl_sf_dilog(-y) + log(y)*log(y)/4. + M_PI*M_PI/12. ) * sqrt(1. - 4.*rho);
1390  else
1391  return - gsl_sf_clausen(2. * asin( 1./(2. * sqrt(rho)) )) * sqrt(4.*rho - 1.);
1392 }
1393 
1394 
1395 double Bsgamma::Y2NV_PHI4(double rho)
1396 {
1397  double y = (1. - sqrt(1. - 4.*rho)) / (1. + sqrt(1. - 4.*rho));
1398  Polylogarithms Poly;
1399  ClausenFunctions Clausen;
1400 
1401  if (rho < 1./4.)
1402  return Poly.Li3(- y) + log(y)*log(y)*log(y)/12. + log(y)*M_PI*M_PI/12.;
1403  else
1404  return Clausen.Cl3(2. * asin( 1./(2. * sqrt(rho)) ));
1405 }
1406 
1407 double Bsgamma::Y2NV(double E0, double mu)
1408 {
1409  double Lb = log(mu/Mb_kin);
1410  double rho = zeta();
1411  double rho2 = rho*rho;
1412  double rho32 = rho*sqrt(rho);
1413  double Lr = log(rho);
1414  double Li2 = gsl_sf_dilog(1. - rho);
1415  double Li2sqrt = gsl_sf_dilog(1. - sqrt(rho));
1416  Polylogarithms Poly;
1417  double Li3 = Poly.Li3(1. - rho);
1418  double Li3ov = Poly.Li3(1. - 1./rho);
1419 
1420  return -16./81. * (157. - 279.*rho - M_PI*M_PI*(5. + 9.*rho2 - 42.*rho32))
1421  -64./27. * (18. - M_PI*M_PI)*Lb - 64./9. * Lb*Lb
1422  +16./27. * (22. - M_PI*M_PI + 10.*rho)*Lr + 16./27. * (8. + 9.*rho2)*Lr*Lr
1423  -16./27. * Lr*Lr*Lr - 8./9. * (1. - 6.*rho2)*Y2NV_PHI1(rho)
1424  -8./27. * (19. - 46.*rho)*Y2NV_PHI2(rho) - 32./27. * (13. + 14.*rho)*Y2NV_PHI3(rho)
1425  -64./9. * Y2NV_PHI4(rho) - 32./9. * Lr*Li2
1426  +32./27. * (5. + 9.*rho2 + 14.*rho32)*Li2 - 1792./27. * rho32*Li2sqrt
1427  +64./9. * Li3 + 64./9. * Li3ov + 4./3.*(2.*Lb - Lr)*dY1(E0);
1428 }
1429 
1430 double Bsgamma::Y2NH(double E0, double mu)
1431 {
1432  double Lb = log(mu/Mb_kin);
1433  double zeta3 = gsl_sf_zeta_int(3);
1434  double Cl2 = gsl_sf_clausen(M_PI/3.);
1435 
1436  return 8./81. * (244. - 27.*sqrt(3.)*M_PI - 61.*M_PI*M_PI)
1437  - 64./27.*(18. - M_PI*M_PI)*Lb - 64./9.*Lb*Lb - 64./27.*zeta3
1438  + 32.*sqrt(3.)*Cl2 + 8./3.*Lb*dY1(E0);
1439 }
1440 
1441 double Bsgamma::Y2(double E0, double mu)
1442 {
1443  double CF = 4./3.;
1444  double CA = 3.;
1445  double TR = 1./2.;
1446  double NL = 3.;
1447  double NV = 1.;
1448  double NH = 1.;
1449 
1450  return CF*Y2CF(E0,mu) + CA*Y2CA(E0,mu)
1451  + TR * ( NL*Y2NL(E0,mu) + NV*Y2NV(E0,mu) + NH*Y2NH(E0,mu) );
1452 }
1453 
1454 double Bsgamma::f_NLO_1(double z)
1455 {
1456  return r1(2,z).real() + 2.*Phi27_1(0.,z).real();
1457 }
1458 
1459 double Bsgamma::zdz_f_NLO(double z, double E0)
1460 {
1461  double d = delta(E0);
1462  double sqrt1d = sqrt(1. - d);
1463  double sqrt4z = sqrt(1. - 4. * z);
1464  double sqrt1d4z = sqrt(1. - d - 4. * z);
1465  double sqrtz = sqrt(z);
1466  double sqrt1ovz = sqrt(1./z);
1467  double sqrt4m1ovz = sqrt(-4. + 1./z);
1468  double SumSqrt = sqrt1d + sqrt1d4z;
1469  double ProdSqrt = sqrt1d * sqrt1d4z;
1470  double ProdSqrtz = sqrtz * sqrt1d4z;
1471  double LogSumSqrt = log(SumSqrt/(2. * sqrtz));
1472  double LogSqrt4z = log((1. + sqrt4z)/(2. * sqrtz));
1473  double LogSqrtov = log((sqrt4m1ovz + sqrt1ovz)/2.);
1474 
1475  double z2=z*z;
1476  double z3=z2*z;
1477  double z4=z3*z;
1478  double z5=z4*z;
1479  double Lz = log(z);
1480 
1481  double Pi2 = M_PI*M_PI;
1482  double zeta3 = gsl_sf_zeta_int(3);
1483 
1484  double zdz_f_NLO_E0;
1485 
1486  if (E0 == 0. ){
1487  zdz_f_NLO_E0 = 2./27. * (3. * (-7. + 2. * M_PI*M_PI)
1488  + 2. * (36. - 24. * M_PI*M_PI) * z
1489  + 108. * M_PI*M_PI * z2
1490  - ( 36. * (-pow(1./z,3./2.)/2. - 1./(2. * sqrt4m1ovz * z2))
1491  * sqrt4m1ovz * (-1. + 2. * z))/((sqrt4m1ovz + sqrt1ovz) * pow(1./z,3./2.))
1492  - (72. * sqrt4m1ovz * LogSqrtov)/pow(1./z,3./2.)
1493  - (54. * sqrt4m1ovz * (-1. + 2. * z) * LogSqrtov)/sqrt1ovz
1494  + (18. * sqrt1ovz * (-1. + 2. * z) * LogSqrtov)/sqrt4m1ovz
1495  - (48. * (-pow(1./z,3./2.)/2. - 1./(2. * sqrt4m1ovz * z2))
1496  * z * (1. - 4. * z + 6. * z2) * LogSqrtov)/(sqrt4m1ovz + sqrt1ovz)
1497  - 24. * z * (-4. + 12. * z) * LogSqrtov * LogSqrtov
1498  - 24. * (1. - 4. * z + 6. * z2) * LogSqrtov * LogSqrtov);
1499  } else zdz_f_NLO_E0 = 2. * ((2. - d) * d * (-7. + 2. * Pi2) / 9.
1500  + 8./9. * d * (3. - 2. * Pi2) * z
1501  + (8. * d * (-SumSqrt/( 4. * pow(z,3./2.) ) - 1./ProdSqrtz)
1502  * ProdSqrt * pow(z,3./2.))/(3. * SumSqrt)
1503  + 4./3. * d * ProdSqrt * LogSumSqrt
1504  - (8. * (1. - d) * d * z * LogSumSqrt)/(3. * ProdSqrt)
1505  - (32. * d * (-SumSqrt/( 4. * pow(z,3./2.) ) - 1./ProdSqrtz)
1506  * (2. - d - 4. * z) * pow(z,3./2.) * LogSumSqrt)/(9. * SumSqrt)
1507  - 8./9. * d * (2. - d - 4. * z) * LogSumSqrt * LogSumSqrt
1508  + 32./9. * d * z * LogSumSqrt * LogSumSqrt
1509  + 4./3. * (1. - 2. * z) * z *
1510  ((2. * (-( (1. + sqrt4z)/(4. * pow(z,3./2.)) )
1511  - 1./(sqrt4z * sqrtz)) * sqrt4z * sqrtz)/(1. + sqrt4z)
1512  - (2. * ( -(SumSqrt/(4. * pow(z,3./2.)) ) - 1./ProdSqrtz)
1513  * ProdSqrt * sqrtz)/(SumSqrt) - 2. * LogSqrt4z/sqrt4z
1514  + 2. * (1. - d) * LogSumSqrt/ProdSqrt)
1515  + 4./3. * (1. - 2. * z) * (sqrt4z * LogSqrt4z - ProdSqrt * LogSumSqrt)
1516  - 8./3. * z * (sqrt4z * LogSqrt4z - ProdSqrt * LogSumSqrt)
1517  - 8./9. * z * (1 - 4. * z + 6. * z2) *
1518  (( 4. * (-( (1. + sqrt4z)/(4. * pow(z,3./2.)) )
1519  - 1./(sqrt4z * sqrtz)) * sqrtz * LogSqrt4z)/(1. + sqrt4z)
1520  - (4. * (-SumSqrt/(4. * pow(z,3./2.)) - 1./ProdSqrtz) * sqrtz * LogSumSqrt)/SumSqrt)
1521  - 8./9. * (LogSqrt4z*LogSqrt4z - LogSumSqrt*LogSumSqrt)
1522  * ( z * (-4. + 12. * z) + (1. - 4. * z + 6. * z2) )) ;
1523 
1524  return z * (zdz_f_NLO_E0
1525 
1526  - 16./9. * (-4. + Pi2/6. - Pi2 * sqrtz
1527  - Lz - z2 * (19./18. - 4. * Lz) + z * (-2. - Lz)
1528  + z3 * (137./30. + 4. * Lz)
1529  + z4 * (887./84. + 10. * Lz)
1530  + z5 * (16597./540. + 28. * Lz)
1531  - 3. * z2 * (25./12. + Pi2/9. + 19. * Lz/18. - 2. * Lz * Lz)
1532  + 2. * z * (1./2. + Pi2 - 2. * Lz - Lz * Lz/2)
1533  + 4. * z3 * (-1376./225. + 2. * Pi2/3 + 137. * Lz/30. + 2. * Lz * Lz)
1534  + 5. * z4 * (-131317./11760. + 5. * Pi2/3. + 887. * Lz/84. + 5. * Lz * Lz)
1535  + 6. * z5 * (-2807617./97200. + 14. * Pi2/3. + 16597. * Lz/540. + 14. * Lz * Lz))
1536 
1537  + 32./9. * (5./2. - Pi2/3. + (5./2. - 3. * Pi2/4.) * Lz
1538  + Lz * Lz/4 + Lz * Lz * Lz/12.
1539  + z5 * (-3303./800. - 63. * Lz/10.)
1540  + z4 * (-185./144. - 35. * Lz/12.)
1541  + z3 * (-1./72. - 5. * Lz/3.)
1542  + z2 * (2. - 3. * Lz/2.)
1543  + 6. * z5 * (67801./8000. - 21. * Pi2/20.
1544  - 3303. * Lz/800. - 63. * Lz * Lz/20.)
1545  + 5. * z4 * (35101./8640. - 35. * Pi2/72.
1546  - 185. * Lz/144. - 35. * Lz * Lz/24.)
1547  + 4. * z3 * (457./216. - 5. * Pi2/18. - Lz/72. - 5. * Lz * Lz/6.)
1548  + 3. * z2 * (-7./6. - Pi2/4. + 2. * Lz - 3. * Lz * Lz/4.)
1549  + z * (-Pi2/2. - Lz/2. + Lz * Lz/4.)
1550  + 5./2. - 3. * Pi2/4. + Lz/2. + Lz * Lz/4.
1551  + 2. * z * (7./4. + 2. * Pi2/3. - Pi2 * Lz/2. - Lz * Lz/4.
1552  + Lz * Lz * Lz/12.) - 3. * zeta3));
1553 }
1554 
1555 double Bsgamma::mddel_f_NLO(double z, double E0)
1556 {
1557  double d = delta(E0);
1558  double d2 = d*d;
1559  double sqrt1d = sqrt(1. - d);
1560  double sqrt1d4z = sqrt(1. - d - 4. * z);
1561  double sqrtz = sqrt(z);
1562  double LogSqrt = log((sqrt1d + sqrt1d4z)/(2. * sqrtz));
1563  double SumSqrt = sqrt1d + sqrt1d4z;
1564  double ProdSqrt = sqrt1d * sqrt1d4z;
1565 
1566  return 2. * (1. - d) * ( -2./27. * d * (-3. + 2. * d)
1567  - 2./27. * (3. - 3. * d + d2)
1568  + 1./9. * (2. - d) * (-7. + 2. * M_PI * M_PI) * z
1569  - 1./9. * d * (-7. + 2. * M_PI * M_PI) * z
1570  + 4. * d * (-1./(2. * sqrt1d) - 1./(2. * sqrt1d4z))
1571  * ProdSqrt * z / (3. * SumSqrt)
1572  + 4./9. * (3. - 2. * M_PI * M_PI) * z * z
1573  + 4./3. * ProdSqrt * z * LogSqrt
1574  - 16. * d * (-1./(2. * sqrt1d) - 1./(2. * sqrt1d4z))
1575  * (2. - d - 4. * z) * z * LogSqrt / (9. * SumSqrt)
1576  + 2. * d * z * (-2. + 2. * d + 4. * z) * LogSqrt / (3. * ProdSqrt)
1577  + 16. * (-1./(2. * sqrt1d) - 1./(2. * sqrt1d4z)) * z
1578  * (1 - 4. * z + 6. * z * z) * LogSqrt / (9. * SumSqrt)
1579  + 8./9. * d * z * LogSqrt * LogSqrt
1580  - 8./9. * (2. - d - 4. * z) * z * LogSqrt * LogSqrt
1581  + 4./3. * (1. - 2. * z) * z
1582  * ( (1./(2. * sqrt1d) + 1./(2. * sqrt1d4z)) * ProdSqrt/SumSqrt
1583  - (-2. + 2. * d + 4. * z) * LogSqrt/(2. * ProdSqrt)));
1584 }
1585 
1586 double Bsgamma::h27_2(double z, double E0)
1587 {
1588  double d = delta(E0);
1589  double d2 = d*d;
1590 
1591  if (E0 == 0.){
1592  return ( 41./27. - 2./9. * M_PI*M_PI
1593  - 2.24 * sqrt(z) - 7.04 * z + 23.72 * pow(z,3./2.)
1594  + ( -9.86 * z + 31.28 * z * z ) * log(z));
1595  } else return - 0.1755402735503456 - 1.4553730660088837 * d
1596  + 1.1192806367180177 * d2
1597  + ( 0.7259818237183779 - 7.230418135384073 * d
1598  + 5.977206932166958 * d2 ) * sqrt(z)
1599  + ( 13.786205094458156 + 113.71026116073105 * d
1600  - 100.3588074342665 * d2 ) * z
1601  + ( -145.05588751363894 - 307.05884309429547 * d
1602  + 388.54181686721904 * d2 ) * pow(z,3./2.)
1603  + ( 475.2039505292043 + 312.9832308573048 * d
1604  - 775.8088176670707 * d2 ) * z * z
1605  + ( -509.7299390734172 - 126.08888075477071 * d
1606  + 646.2084041395774 * d2 ) * pow(z,5./2.);
1607 }
1608 
1609 double Bsgamma::f_q(double z, double E0)
1610 {
1611  return Rer22(z) - 4./3. * h27_2(z,E0);
1612 }
1613 
1614 double Bsgamma::f_b(double z)
1615 {
1616  return -1.836 + 2.608 * z + 0.8271 * z * z - 2.441 * z * log(z);
1617 }
1618 
1619 double Bsgamma::f_c(double z)
1620 {
1621  return 9.099 + 13.20 * z - 19.68 * z * z + 25.71 * z * log(z);
1622 }
1623 
1624 double Bsgamma::F_1(double z)
1625 {
1626  return - 23.74697061848885 + 35./12. * f_q(z,0.)
1627  + (2129./936. - 9./52. * M_PI*M_PI) * f_NLO_1(z)
1628  - 0.8444138663102 * zdz_f_NLO(z,0.);
1629 }
1630 
1631 double Bsgamma::F_2(double z)
1632 {
1633  return - 3.006537367876035 - 592./81. * f_q(z,0.)
1634  - 10.344289655256379 * f_NLO_1(z)
1635  - 9.550817514525745 * zdz_f_NLO(z,0.);
1636 }
1637 
1638 double Bsgamma::delddel_Phi22_1(double E0)
1639 {
1640  double d = delta(E0);
1641 
1642  return 4. * (1. - d)/d * 16./27. * Int_cc1_part1(E0);
1643 }
1644 
1645 double Bsgamma::zdz_Phi22_1(double E0)
1646 {
1647  return 64./27. * ( Int_cc1(E0) - Int_cc(E0));
1648 }
1649 
1650 double Bsgamma::delddel_Phi28_1(double z, double E0)
1651 {
1652  double d = delta(E0);
1653  double Sq = sqrt( (1. - d) * (1. - d - 4.*z) );
1654  double Log = log( ( sqrt(1. - d) + sqrt(1. - d - 4.*z) ) / 2. / sqrt(z) );
1655  double Log2 = Log*Log;
1656 
1657  return 4. / (27. * Sq * (1 - d - 4. * z)) * Sq * Sq *
1658  (-8. * (-1. + d) * Log * z * (-1. + d + 4. * z) +
1659  Sq * (1. + d*d + (4. + 8. * Log2 - 2. * M_PI*M_PI) * z +
1660  4. * (-4. * Log2 + M_PI*M_PI) * z * z -
1661  2. * d * (1. + (2. + 4. * Log2 - M_PI*M_PI) * z)));
1662 }
1663 
1664 double Bsgamma::zdz_Phi28_1(double z, double E0)
1665 {
1666  double d = delta(E0);
1667  double Sq = sqrt( (1. - d) * (1. - d - 4.*z) );
1668  double Log1 = log( ( sqrt(1. - d) + sqrt(1. - d - 4.*z) ) / 2. / sqrt(z) );
1669  double Log2 = log( ( 1. + sqrt(1. - 4.*z) ) / 2. / sqrt(z) );
1670 
1671  return 2./27. * z * (d*d * (-7. - 8. * (-1. + Log1) * Log1 + 2. * M_PI*M_PI)
1672  - 20. * Log2 * sqrt(1. - 4. * z) + 72. * Log2 * sqrt(1. - 4. * z) * z
1673  + (48. * Log1 * Sq * z * z) / (-1. + d + 4. * z)
1674  + (8. * z * ( -(3. + 4. * Log1) * Sq*Sq + 2. * (3. + 8. * Log1) * Sq * z
1675  - 24. * Log1 * z * z )) / (-1. + d - Sq + 4. * z)
1676  - 2. * d * (-7. - 3. * Sq + Log1 * (8. + 6. * Sq) + M_PI*M_PI * (2. - 8. * z)
1677  + 12. * z + 8. * Log1*Log1 * (-1. + 4. * z))
1678  + 2. * (-3. * (-1. + Sq + 2. * (1. + Sq) * z)
1679  + Log1 * (4. + Sq * (6. - 52. * z) + 24. * z * z)
1680  + 4. * ( Log2 * Log2 - Log1 * Log1) * (1. + 2. * z * (-4. + 9. * z))));
1681 }
1682 
1683 double Bsgamma::delddel_Phi88_1(double E0)
1684 {
1685  double d = delta(E0);
1686  double Ld = log(d);
1687 
1688  return 4./27. * (1. - d) * (5. - 8./(1. - d) + 5. * d - 2. * d * d -
1689  2. * (2. - 4./(1. - d) + 2. * d) * log(Mb_kin/Ms) + (4. * Ld)/(1. - d)
1690  - d * Ld - (2. + d) * Ld);
1691 }
1692 
1693 double Bsgamma::f(double r)
1694 {
1695  double r2 = r*r;
1696  double r3 = r2*r;
1697  double r4 = r3*r;
1698  double Lr = log(r);
1699  double zeta3 = gsl_sf_zeta_int(3);
1700  Polylogarithms Poly;
1701 
1702  if (r==1.)
1703  return 7126./81. - 356./27. * M_PI*M_PI - 16./3. * zeta3;
1704  else
1705  return - 16./3. * Poly.Li3(r2)
1706  + 8. * r * (35./9. * r2 + 9.) * ( gsl_sf_dilog(r) - atanh(r) * Lr - 1./4. * M_PI*M_PI)
1707  + 2. * (8./3. * Lr - 6. * r4 - 35./9. * r3 - 9. * r - 62./9. ) * gsl_sf_dilog(r2)
1708  - 8. * (3. * r4 + 31./9.) * log(1. - r2) * Lr + 32./9. * Lr*Lr*Lr
1709  + 8. * (3. * r4 + 25./9.) * Lr*Lr + 64./9. * r2 * Lr
1710  + 4. * M_PI*M_PI * r4 + 172./9. * r2 + 5578./81.;
1711 }
1712 
1713 double Bsgamma::Delta(double r)
1714 {
1715  double r2 = r*r;
1716  double r3 = r2*r;
1717  double Lmr = log(1. - r);
1718  double Lpr = log(1. + r);
1719  double Lr = log(r);
1720  double Lr2 = Lr*Lr;
1721 
1722  if (r==1.)
1723  return -3./8. + 1./8. * M_PI*M_PI;
1724  else
1725  return 1./4. * (1. - r) * (1. - r3) * ( gsl_sf_dilog(r) + Lr * Lmr - 1./2. * Lr2 - 1./3. * M_PI*M_PI )
1726  - 1./4. * (1. + r) * (1. + r3) * ( gsl_sf_dilog(-1./r) - Lr * Lpr + Lr2 )
1727  + 1./4. * Lr2 - 1./4. * r2 * Lr - 3./8. * r2 + 1./24. * M_PI*M_PI;
1728 }
1729 
1730 double Bsgamma::f_u(double r)
1731 {
1732  double r2 = r*r;
1733  double r3 = r2*r;
1734  double r4 = r3*r;
1735  double r5 = r4*r;
1736  double r6 = r5*r;
1737  double r7 = r6*r;
1738  double Lr = log(r);
1739  double Lr2 = Lr*Lr;
1740  double Pi2 = M_PI*M_PI;
1741  double zeta3 = gsl_sf_zeta(3.);
1742 
1743  if (r==1.)
1744  return 6335./288. - 1./2. * M_PI*M_PI - 16. * zeta3;
1745  else
1746  return -5./6. * Pi2 * r + ( 14. + 16./9. * Pi2) * r2
1747  + (64./9. * Lr + 128./9. * log(2.) - 95./54.) * Pi2 * r3
1748  + (-16./3. * Lr2 + 365./9. * Lr + 4. * Pi2 * Lr
1749  - 4375./54. - 25./9. * Pi2 + 32. * zeta3) * r4
1750  - 224./45. * Pi2 * r5 + (-128./27. * Lr2 + 16./15. * Lr
1751  + 15608./2025. + 128./81. * Pi2) * r6 - 16./7. * Pi2 * r7;
1752 }
1753 
1754 double Bsgamma::omega77(double z)
1755 {
1756  double z2 = z*z;
1757  double z3 = z2*z;
1758  double z4 = z3*z;
1759  double z5 = z4*z;
1760  double z6 = z5*z;
1761  double z7 = z6*z;
1762  double z8 = z7*z;
1763  double omz = 1. - z;
1764  double omz2 = omz * omz;
1765  double omz3 = omz2 * omz;
1766  double Lz = log(z);
1767  double Lomz = log (1. - z);
1768  double Lomz2 = Lomz*Lomz;
1769  double Lomz3 = Lomz2*Lomz;
1770  double Ltmz = log (2. - z);
1771  double Ltmz2 = Ltmz*Ltmz;
1772  double Li2omz = gsl_sf_dilog(1. - z);
1773  double Li2zmo = gsl_sf_dilog(z - 1.);
1774  Polylogarithms Poly;
1775  double Pi2 = M_PI*M_PI;
1776 
1777  return 4./9. * (z3 - 4. * z2 + 4. * z + 1.)/omz * ( 2. * Poly.Li3( 1./(2.-z) )
1778  - Poly.Li3( z/(2.-z) ) + Poly.Li3( z/(z-2.) )
1779  + Ltmz * ( Lomz2 - 1./3. * Ltmz2 + 1./6. * Pi2 ) )
1780  + 4./9. * (z3 + 36. * z - 43.)/omz * Poly.Li3(z)
1781  + 8./9. * (z3 - 2. * z2 + 19. * z - 22.)/omz * Poly.Li3(1.-z)
1782  - 16./9. * omz2 * Poly.Li3(z-1.)
1783  - 4./9. * (z3 + 35. * z - 44.)/omz * Li2omz * Lomz
1784  - 4./9. * (z3 - 2. * z2 + 2. * z - 3.)/omz * Li2zmo * Lomz
1785  - 4./27. * (23. * z6 - 106. * z5 + 145. * z4 + 3. * z3
1786  - 180. * z2 + 147. * z - 36.)/(z * omz3) * (Li2omz + Lomz * Lz)
1787  + 2./27. * (z8 - 6. * z7 + 9. * z6 + 27. * z5 - 140. * z4 + 219. * z3
1788  - 124. * z2 + 28. * z - 6.)/(z * omz3) * (Li2zmo + Lomz * Ltmz)
1789  - 8./9. * (z2 + 8. * z - 11.)/omz * Lomz2 * Lz
1790  - 2./9. * (z4 - 3. * z3 - 5. * z2 + 15. * z + 8.)/(z * omz) * Lomz3
1791  - (z6 - 4. * z5 - 46. * z4 + 101. * z3 - 461. * z2 + 1057. * z - 72.)/(27. * z * omz) * Lomz2
1792  + 2./27. * (z3 - 2. * z2 + 4. * z - 5.)/omz * Pi2 * Lomz
1793  + (2. * z5 - 29. * z4 - 113. * z3 + 153. * z2 - 827. * z - 162.)/(27. * z * omz) * Lomz
1794  - (3. * z3 - 8. * z2 + 144. * z - 157.)/(9. * omz) * gsl_sf_zeta(3.)
1795  + (z6 - 4. * z5 + 48. * z4 - 106. * z3 - 58. * z2 + 158. * z - 75.)/(81. * z * omz) * Pi2
1796  + (2. * z4 - 92. * z3 + 88. * z2 - 713. * z - 18.)/(27. * omz);
1797 }
1798 
1799 double Bsgamma::Int_Phi77_2rem(double E0)
1800 {
1801  if (IntPhi772rCached == 0) {
1802  double t1 = (1. - delta(E0));
1803 
1804  INT = convertToGslFunction(boost::bind(&Bsgamma::omega77, &(*this), _1));
1805  if (gsl_integration_cquad(&INT, 0., t1, 1.e-2, 1.e-1, w_INT, &avaINT, &errINT, NULL) != 0) return std::numeric_limits<double>::quiet_NaN();
1806 
1807  CacheIntPhi772r = avaINT;
1808  IntPhi772rCached = 1;
1809  }
1810 
1811  return CacheIntPhi772r;
1812 }
1813 
1814 double Bsgamma::Phi77_2rem(double E0)
1815 {
1816  double xm = 8./9. * M_PI * alsUps;
1817  double d = delta(E0);
1818  double d2 = d*d;
1819  double d3 = d2*d;
1820  double d4 = d3*d;
1821 
1822  return Int_Phi77_2rem(E0) - xm/3./ d *
1823  ((4. - 6. * d2 + 2. * d3) * log(d) + 7. - 13. * d + 3. * d2 + 5. * d3 - 2. * d4);
1824 }
1825 
1826 double Bsgamma::K77_2_z1(double E0, double mu)
1827 {
1828  double K77_1 = Kij_1(7,7,E0,mu).real();
1829  double Pi2 = M_PI*M_PI;
1830  double xm = 8./9. * M_PI * alsUps;
1831  double Lb = 2.*log(mu_b/Mb_kin);
1832 
1833  return ( K77_1 - 4. * Phi77_1(E0) ) * K77_1 - 1178948./729. + 18593./729. * Pi2
1834  - 628./405. * Pi2*Pi2 + 428./27. * Pi2 * log(2.) + 61294./81. * gsl_sf_zeta(3.)
1835  - 880./9. * Lb * Lb + ( 440./27. * Pi2 - 14698./27. ) * Lb
1836  + 64./3. * xm + 4. * (Phi77_2beta0(E0,mu) + Phi77_2rem(E0));
1837 }
1838 
1839 double Bsgamma::Kij_2(int i, int j, double E0, double mu_b, double mu_c)
1840 {
1841  if (i<1 || i > 2)
1842  if (i < 7)
1843  throw std::runtime_error("Bsgamma::Kij_2(): index i must be included in (1,2,7,8)");
1844  if (j<1 || j > 2)
1845  if (j < 7)
1846  throw std::runtime_error("Bsgamma::Kij_2(): index j must be included in (1,2,7,8)");
1847 
1848  int temp;
1849 
1850  if (i > j) {temp=i; i=j; j=temp;}
1851 
1852  double K_ij[8][8] = {{0.}};
1853  double z = zeta();
1854  double d = delta(E0);
1855  double r = sqrt(z);
1856  double Lb = 2.*log(mu_b/Mb_kin);
1857  double Lb2 = Lb*Lb;
1858  double Lc = 2.*log(mu_c/Mc);
1859  double Lcb = log(Mc/Mb_kin);
1860  double xm = 8./9. * M_PI * alsUps;
1861 
1862  double A1 = 22.604961348474838;
1863  double A2 = 75.60281607240395;
1864 
1865  K_ij[1][1] = (r1(2,zeta()).real() - 208./81.*Lb) * (r1(2,zeta()).real() - 208./81.*Lb)
1866  + r1(2,zeta()).imag() * r1(2,zeta()).imag()
1867  + 4.*Phi22_2beta0(E0,mu_b) * SM.Beta0(3)/SM.Beta0(5) + 16./3. * Phi22_1(E0) * (Lcb - Lb)
1868  + xm * (delddel_Phi22_1(E0) - 2. * zdz_Phi22_1(E0));
1869  K_ij[1][6] = A2 + F_2(z) - 27./2. * f_q(z,E0) + f_b(z) + f_c(z)
1870  + 4./3. * Phi27_1(E0,z).real() * log(z) + (8. * Lc - 2. * xm) * zdz_f_NLO(z,E0)
1871  + xm * mddel_f_NLO(z,E0) + 416./81. * xm
1872  + (10./3. * Kij_1(2,7,E0,mu_b).real() - 2./3. * Kij_1(4,7,E0,mu_b).real()
1873  - 208./81. * Kij_1(7,7,E0,mu_b).real() - 35./27. * Kij_1(7,8,E0,mu_b).real()
1874  - 254./81.) * Lb - 5948./729. * Lb2;
1875  K_ij[1][7] = (r1(2,zeta()).real() - 208./81.*Lb) * (r1(8,zeta()).real() + 16./9.*Lb)
1876  + r1(2,zeta()).imag() * r1(8,zeta()).imag()
1877  + 2.*Phi28_2beta0(E0,mu_b) * SM.Beta0(3)/SM.Beta0(5) + 8./3. * Phi28_1(E0,z).real() * (Lcb - Lb)
1878  + xm * (delddel_Phi28_1(z,E0) - 2. * zdz_Phi28_1(z,E0));
1879 
1880  K_ij[0][0] = 1./36. * K_ij[1][1];
1881  K_ij[0][1] = -1./6. * K_ij[1][1];
1882  K_ij[0][6] = - 1./6. * K_ij[1][6] + A1 + F_1(z)
1883  + (- 3./2. * Kij_1(2,7,E0,mu_b).real() - 3./4. * Kij_1(7,8,E0,mu_b).real()
1884  + 94./81.) * Lb - 34./27. * Lb2;
1885  K_ij[0][7] = -1./6. * K_ij[1][7];
1886 
1887  K_ij[6][6] = K77_2_z1(E0,mu_b) + ( 1972./81. - 16./27. * M_PI*M_PI + 8./3. * Phi77_1(E0)) * log(zeta())
1888  + 2./3. * (f(r) - f(1.)) - 128./3. * (Delta(r) - Delta(1.))
1889  - 16. * (f_u(r) - f_u(1.));
1890  K_ij[6][7] = 2./3. * Y2(E0,mu_b) + (16./9.*M_PI*M_PI - 164./9. - 32./6. * Lb) * Y1(E0,mu_b)
1891  - 32./81. * alsUps * M_PI * (3. + 7.*d - 3.*d*d + d*d*d - 4.*d*log(d) );
1892 
1893  K_ij[7][7] = (r1(8,zeta()).real() + 16./9.*Lb) * (r1(8,zeta()).real() + 16./9.*Lb)
1894  + r1(8,zeta()).imag() * r1(8,zeta()).imag()
1895  + 4.*Phi88_2beta0(E0,mu_b) * SM.Beta0(3)/SM.Beta0(5) + 16./3. * Phi88_1(E0) * (Lcb - Lb)
1896  + xm * (delddel_Phi88_1(E0));
1897 
1898  return K_ij[i-1][j-1];
1899 }
1900 
1901 void Bsgamma::computeCoeff(double mu)
1902 {
1903 
1904  /*allcoeff = SM.getMyFlavour()->ComputeCoeffsgamma(160.);
1905 
1906  C1_0 = (*(allcoeff[LO]))(0);
1907  C2_0 = (*(allcoeff[LO]))(1);
1908  C3_0 = (*(allcoeff[LO]))(2);
1909  C4_0 = (*(allcoeff[LO]))(3);
1910  C5_0 = (*(allcoeff[LO]))(4);
1911  C6_0 = (*(allcoeff[LO]))(5);
1912  C7_0 = (*(allcoeff[LO]))(6);
1913  C8_0 = (*(allcoeff[LO]))(7);
1914 
1915  C1_1 = (*(allcoeff[NLO]))(0)/Alstilde;
1916  C2_1 = (*(allcoeff[NLO]))(1)/Alstilde;
1917  C3_1 = (*(allcoeff[NLO]))(2)/Alstilde;
1918  C4_1 = (*(allcoeff[NLO]))(3)/Alstilde;
1919  C5_1 = (*(allcoeff[NLO]))(4)/Alstilde;
1920  C6_1 = (*(allcoeff[NLO]))(5)/Alstilde;
1921  C7_1 = (*(allcoeff[NLO]))(6)/Alstilde;
1922  C8_1 = (*(allcoeff[NLO]))(7)/Alstilde;
1923 
1924  std::cout << "C_0(MuW): (" << C1_0.real() << "," << C2_0.real() << ","
1925  << C3_0.real() << "," << C4_0.real() << "," << C5_0.real() << ","
1926  << C6_0.real() << "," << C7_0.real() << "," << C8_0.real() << ")" << std::endl;
1927  std::cout << "C_1(MuW): (" << C1_1.real() << "," << C2_1.real() << ","
1928  << C3_1.real() << "," << C4_1.real() << "," << C5_1.real() << ","
1929  << C6_1.real() << "," << C7_1.real() << "," << C8_1.real() << ")" << std::endl << std::endl;*/
1930 
1931  allcoeff = SM.getFlavour().ComputeCoeffsgamma(mu);
1932  allcoeffprime = SM.getFlavour().ComputeCoeffprimesgamma(mu);
1933 
1934  C1_0 = (*(allcoeff[LO]))(0);
1935  C2_0 = (*(allcoeff[LO]))(1);
1936  C3_0 = (*(allcoeff[LO]))(2);
1937  C4_0 = (*(allcoeff[LO]))(3);
1938  C5_0 = (*(allcoeff[LO]))(4);
1939  C6_0 = (*(allcoeff[LO]))(5);
1940  C7_0 = (*(allcoeff[LO]))(6) + C_7_NP;
1941  C8_0 = (*(allcoeff[LO]))(7);
1942 
1943  C1_1 = (*(allcoeff[NLO]))(0)/Alstilde;
1944  C2_1 = (*(allcoeff[NLO]))(1)/Alstilde;
1945  C3_1 = (*(allcoeff[NLO]))(2)/Alstilde;
1946  C4_1 = (*(allcoeff[NLO]))(3)/Alstilde;
1947  C5_1 = (*(allcoeff[NLO]))(4)/Alstilde;
1948  C6_1 = (*(allcoeff[NLO]))(5)/Alstilde;
1949  C7_1 = (*(allcoeff[NLO]))(6)/Alstilde;
1950  C8_1 = (*(allcoeff[NLO]))(7)/Alstilde;
1951 
1952  C7_2 = (*(allcoeff[NNLO]))(6)/Alstilde/Alstilde;
1953 
1954  C7p_0 = (*(allcoeffprime[LO]))(6) + Ms/Mb*((*(allcoeff[LO]))(6)) + C_7p_NP;
1955  C7p_1 = ((*(allcoeffprime[NLO]))(6) + Ms/Mb*((*(allcoeff[NLO]))(6)))/Alstilde; /*Implement the other WCs*/
1956 
1957  /*std::cout << "C_0(mu): (" << C1_0.real() << "," << C2_0.real() << ","
1958  << C3_0.real() << "," << C4_0.real() << "," << C5_0.real() << ","
1959  << C6_0.real() << "," << C7_0.real() << "," << C8_0.real() << ")" << std::endl;
1960  std::cout << "C_1(mu): (" << C1_1.real() << "," << C2_1.real() << ","
1961  << C3_1.real() << "," << C4_1.real() << "," << C5_1.real() << ","
1962  << C6_1.real() << "," << C7_1.real() << "," << C8_1.real() << ")" << std::endl << std::endl;
1963  std::cout << "C_2^7(mu): " << C7_2.real() << std::endl << std::endl;*/
1964 
1965  C7_1ew = 4.868;
1966 
1967 }
1968 
1969 double Bsgamma::P0(double E0)
1970 {
1971  return C7_0.abs2() + C7p_0.abs2() + P0_4body(E0,Mb_kin*Mb_kin/Ms/Ms);
1972 }
1973 
1974 double Bsgamma::P11()
1975 {
1976  return 2.*( C7_0.real()*C7_1.real() + C7_0.imag()*C7_1.imag()
1977  + C7p_0.real()*C7p_1.real() + C7p_0.imag()*C7p_1.imag() ); /*CHECK SIGN*/
1978 }
1979 
1980 double Bsgamma::P21(double E0, double mu)
1981 {
1982  int i,j;
1983  gslpp::complex C0[8]={C1_0,C2_0,C3_0,C4_0,C5_0,C6_0,C7_0,C8_0};
1984  gslpp::complex C0p[8]={C7p_0}; /*IMPLEMENT OTHER WC*/
1985  double p21=0.;
1986 
1987  for(i=0;i<8;i++)
1988  {
1989  for(j=0;j<8;j++)
1990  {
1991  p21 += ( C0[i].real()*C0[j].real() + C0[i].imag()*C0[j].imag() ) * Kij_1(i+1,j+1,E0,mu).real();
1992  }
1993  }
1994 
1995  for(i=6;i<7;i++) /*CHECK ALGORITHM*/
1996  {
1997  for(j=6;j<7;j++)
1998  {
1999  p21 += (C0p[i].real()*C0p[j].real() + C0p[i].imag()*C0p[j].imag()) * Kij_1(i+1,j+1,E0,mu).real();
2000  }
2001  }
2002 
2003  return p21;
2004 }
2005 
2006 double Bsgamma::P21_CPodd(double E0, double mu)
2007 {
2008  int i,j;
2009  gslpp::complex C0[8]={C1_0,C2_0,C3_0,C4_0,C5_0,C6_0,C7_0,C8_0};
2010  gslpp::complex C0p[8]={C7p_0}; /*IMPLEMENT OTHER WC*/
2011  double p21=0.;
2012 
2013  for(i=0;i<8;i++)
2014  {
2015  for(j=0;j<8;j++)
2016  {
2017  p21 += - ( C0[i].real()*C0[j].imag() - C0[i].imag()*C0[j].real() ) * Kij_1(i+1,j+1,E0,mu).imag();
2018  }
2019  }
2020 
2021  for(i=6;i<7;i++) /*CHECK ALGORITHM*/
2022  {
2023  for(j=6;j<7;j++)
2024  {
2025  p21 += - ( C0p[i].real()*C0p[j].imag() - C0p[i].imag()*C0p[j].real() ) * Kij_1(i+1,j+1,E0,mu).imag();
2026  }
2027  }
2028 
2029  return p21;
2030 }
2031 
2032 double Bsgamma::P12()
2033 {
2034 
2035  return C7_1.abs2() + C7p_1.abs2() + 2.*(C7_0*C7_2).real(); /*CHECK SIGN*/
2036 }
2037 
2038 double Bsgamma::P22(double E0, double mu_b, double mu_c)
2039 {
2040 
2041  int i,j, temp_i,temp_j;
2042  gslpp::complex C0[4]={C1_0,C2_0,C7_0,C8_0};
2043  double p22=0.;
2044 
2045  for(i=0;i<4;i++)
2046  {
2047  for(j=0;j<4;j++)
2048  {
2049  if (i > 1) {
2050  temp_i=i+4;
2051  } else temp_i=i;
2052  if (j > 1) {
2053  temp_j=j+4;
2054  } else temp_j=j;
2055  p22 += (C0[i]*C0[j]).real() * Kij_2(temp_i+1,temp_j+1,E0,mu_b,mu_c);
2056  }
2057  }
2058 
2059  return p22;
2060 }
2061 
2062 double Bsgamma::P32(double E0, double mu)
2063 {
2064 
2065  int i,j;
2066  gslpp::complex C0[8]={C1_0,C2_0,C3_0,C4_0,C5_0,C6_0,C7_0,C8_0};
2067  gslpp::complex C1[8]={C1_1,C2_1,C3_1,C4_1,C5_1,C6_1,C7_1,C8_1};
2068  double p32=0.;
2069 
2070  for(i=0;i<8;i++)
2071  {
2072  for(j=0;j<8;j++)
2073  {
2074  p32 += 2.*(C0[i]*C1[j]).real() * Kij_1(i+1,j+1,E0,mu).real();
2075  }
2076  }
2077 
2078  return p32;
2079 }
2080 
2081 double Bsgamma::EW_NLO(double mu)
2082 {
2083 
2084  if(EWflag) {
2085  double ew_nlo = 0.;
2086  double ga_eff_ew_7[7] = {-832./729., -208./243., -20./243., -176./729., -22712./243., -6272./729., 16./9.};
2087  double Lb = log(mu/Mb_kin);
2088  double Lz = 2. * log(Mz/mu);
2089  gslpp::complex C[7] = {C1_0, C2_0, C3_0, C4_0, C5_0, C6_0, C7_0};
2090  gslpp::complex r[7] = {0.};
2091 
2092  r[0] = r1_ew(1,zeta()) - ga_eff_ew_7[0] * Lb;
2093  r[1] = r1_ew(2,zeta()) - ga_eff_ew_7[1] * Lb;
2094  r[2] = r1_ew(3,zeta()) - ga_eff_ew_7[2] * Lb;
2095  r[3] = r1_ew(4,zeta()) - ga_eff_ew_7[3] * Lb;
2096  r[4] = r1_ew(5,zeta()) - ga_eff_ew_7[4] * Lb;
2097  r[5] = r1_ew(6,zeta()) - ga_eff_ew_7[5] * Lb;
2098  r[6] = r1_ew(7,zeta()) - ga_eff_ew_7[6] * Lb;
2099 
2100  for(int i=0;i<7;i++){
2101  ew_nlo += 2. * (C7_0.real()*C[i].real() + C7_0.imag()*C[i].imag()) * r[i].real()
2102  - 2. * (C7_0.real()*C[i].imag() - C7_0.imag()*C[i].real()) * r[i].imag();
2103  }
2104 
2105  ew_nlo += 2. * (C7_0.real() * (-2. * C7_0 * Lz + C7_1ew).real()
2106  + C7_0.imag() * (-2. * C7_0 * Lz + C7_1ew).imag());
2107 
2108  return ew_nlo;
2109  }
2110 
2111  else return 0.;
2112 }
2113 
2114 double Bsgamma::Vub_NLO_2body()
2115 {
2116  double z = zeta();
2117 
2118  return 4. * Alstilde *
2119  ((C7_0.real()*( C2_0-C1_0/6. ).real() + C7_0.imag()*( C2_0-C1_0/6. ).imag()) * CKMu.real() +
2120  (C7_0.real()*( C2_0-C1_0/6. ).imag() - C7_0.imag()*( C2_0-C1_0/6. ).real()) * CKMu.imag())
2121  * ( a(z)+b(z) ).real();
2122 }
2123 
2124 double Bsgamma::Vub_NLO_2body_CPodd()
2125 {
2126  double z = zeta();
2127 
2128  return - 4. * Alstilde *
2129  ((C7_0.real()*( C2_0-C1_0/6. ).real() + C7_0.imag()*( C2_0-C1_0/6. ).imag()) * CKMu.imag() +
2130  (C7_0.real()*( C2_0-C1_0/6. ).imag() - C7_0.imag()*( C2_0-C1_0/6. ).real()) * CKMu.real() )
2131  * ( a(z)+b(z) ).imag();
2132 }
2133 
2134 double Bsgamma::Vub_NLO_3body_A(double E0)
2135 {
2136  double d = delta(E0);
2137 
2138  return 64./27. * Alstilde * ( C2_0 - C1_0/6. ).abs2() *
2139  ( CKMu.real() * ( 2. * Int_cc1(E0) - Int_c1(E0).real() )
2140  + CKMusq * ( Int_cc1(E0) - Int_c1(E0).real() + 1./8. * d * ( 1. - 1./3. * d*d ) ));
2141 }
2142 
2143 double Bsgamma::Vub_NLO_3body_A_CPodd(double E0)
2144 {
2145  return - 64./27. * Alstilde * ( C2_0 - C1_0/6. ).abs2() * CKMu.imag() * Int_c1(E0).imag();
2146 }
2147 
2148 double Bsgamma::Vub_NLO_3body_B(double E0)
2149 {
2150  double d = delta(E0);
2151  gslpp::complex wc1 = C7_0 - C8_0/3.;
2152  gslpp::complex wc2 = C2_0 - C1_0/6.;
2153  gslpp::complex me = Phi27_1(E0,zeta()) + 2./9. * d * ( 1. - d + 1./3. * d*d );
2154 
2155  return 4. * Alstilde *
2156  (( wc1.real() * wc2.real() + wc1.imag() * wc2.imag()) * CKMu.real() +
2157  ( wc1.real() * wc2.imag() - wc1.imag() * wc2.real()) * CKMu.imag() )
2158  * me.real();
2159 }
2160 
2161 double Bsgamma::Vub_NLO_3body_B_CPodd(double E0)
2162 {
2163  double d = delta(E0);
2164  gslpp::complex wc1 = C7_0 - C8_0/3.;
2165  gslpp::complex wc2 = C2_0 - C1_0/6.;
2166  gslpp::complex me = Phi27_1(E0,zeta()) + 2./9. * d * ( 1. - d + 1./3. * d*d );
2167 
2168  return - 4. * Alstilde *
2169  (( wc1.real() * wc2.real() + wc1.imag() * wc2.imag()) * CKMu.imag() +
2170  ( wc1.real() * wc2.imag() - wc1.imag() * wc2.real()) * CKMu.real() )
2171  * me.imag();
2172 }
2173 
2174 double Bsgamma::Vub_NLO_4body(double E0)
2175 {
2176  if (FOUR_BODY) {
2177  double d = delta(E0);
2178  double d2 = d*d;
2179  double d3 = d2*d;
2180  double Ld = log(d);
2181  double Lumd = log(1. - d);
2182  double Lq = log(Ms/Mb_kin);
2183 
2184  double uphib427 = ( 2. * d * (-63. + 30. * d + 35. * d2 - 2. * d3
2185  + 3. * d * (-18. - 7. * d + d2) * Ld) ) / ( 243. * (d - 1.) );
2186  double uphib428 = ( 108. * (d - 1.) * (d - 1.) * Lumd*Lumd
2187  - 12. * Lumd * (- 25. - 18. * Lq - 18. * d * (5. + 4. * Lq)
2188  + 9. * d2 * (5. + 4. * Lq) + (9. + 36. * d - 18. * d2) * Ld)
2189  + d * (24. * (17. + 9. * Lq) + 27. * d * (43. + 26. * Lq)
2190  - d2 * (127. + 72. * Lq) + 9. * (-12. - 39. * d + 4. * d2) * Ld)
2191  + 108. * (-1. - 4. * d + 2. * d2) * gsl_sf_dilog(d) ) / 729.;
2192 
2193  return 4. * Alstilde * (
2194  ( C2_0 - C1_0/6. ).abs2() * CKMu.real() * 0.005025213076791178
2195  + (( C2_0 - C1_0/6. ).real() * (C7_0.real() * uphib427 + C8_0.real() * uphib428)
2196  + ( C2_0 - C1_0/6. ).imag() * (C7_0.imag() * uphib427 + C8_0.imag() * uphib428) ) * CKMu.real()
2197  + (( C2_0 - C1_0/6. ).real() * (C7_0.imag() * uphib427 + C8_0.imag() * uphib428)
2198  - ( C2_0 - C1_0/6. ).imag() * (C7_0.real() * uphib427 + C8_0.real() * uphib428) ) * CKMu.imag() );
2199  }
2200 
2201  else return 0.;
2202 }
2203 
2204 double Bsgamma::Vub_NLO_4body_CPodd(double E0)
2205 {
2206  if (FOUR_BODY) {
2207  return 4. * Alstilde * ( C2_0 - C1_0/6. ).abs2() * CKMu.imag() * 0.013978889449487913;
2208  }
2209 
2210  else return 0.;
2211 }
2212 
2213 double Bsgamma::Vub_NLO(double E0)
2214 {
2215  return Vub_NLO_2body() + Vub_NLO_3body_A(E0) + Vub_NLO_3body_B(E0) + Vub_NLO_4body(E0);
2216 }
2217 
2218 double Bsgamma::Vub_NLO_CPodd(double E0)
2219 {
2220  return Vub_NLO_2body_CPodd() + Vub_NLO_3body_A_CPodd(E0) + Vub_NLO_3body_B_CPodd(E0) + Vub_NLO_4body_CPodd(E0);
2221 }
2222 
2223 double Bsgamma::Vub_NNLO(double E0)
2224 {
2225  double r12 = (( C2_1 - C1_1/6. )/( C2_0 - C1_0/6. )).real();
2226  double r78 = (( C7_1 - C8_1/3. )/( C7_0 - C8_0/3. )).real();
2227  double r7 = (C7_1/C7_0).real();
2228 
2229  return Alstilde * ( (r12 + r7) * Vub_NLO_2body()
2230  + 2. * r12 * Vub_NLO_3body_A(E0) + (r12 + r78) * Vub_NLO_3body_B(E0));
2231 }
2232 
2233 double Bsgamma::P(double E0, double mu_b, double mu_c, orders order)
2234 {
2235 
2236  switch(order) {
2237  case NNLO:
2238  /*std::cout << "p0 w/ tree, VubLO: " << P0(E0) << std::endl;
2239  std::cout << "p11: " << P11() << std::endl;
2240  std::cout << "p21: " << P21(E0,mu_b) << std::endl;
2241  std::cout << "p12: " << P12() << std::endl;
2242  std::cout << "p22: " << P22(E0,mu_b,mu_c) << std::endl;
2243  std::cout << "p32: " << P32(E0,mu_b) << std::endl;
2244  std::cout << "Vub_NLO: " << Vub_NLO(E0) << std::endl;
2245  std::cout << "Vub_NNLO: " << Vub_NNLO(E0) << std::endl;
2246  std::cout << "EW_NLO: " << EW_NLO(mu_b) << std::endl;*/
2247  return P0(E0)
2248  + Alstilde * (P11() + P21(E0,mu_b)) + Vub_NLO(E0) + AleatMztilde * EW_NLO(mu_b)
2249  + Alstilde * Alstilde * (P12() + P22(E0,mu_b,mu_c) + P32(E0,mu_b)) + Vub_NNLO(E0);
2250  break;
2251  case NLO:
2252  return P0(E0)
2253  + Alstilde * (P11() + P21(E0,mu_b)) + Vub_NLO(E0) + AleatMztilde * EW_NLO(mu_b);
2254  break;
2255  case LO:
2256  return P0(E0);
2257  break;
2258  default:
2259  std::stringstream out;
2260  out << order;
2261  throw std::runtime_error("Bsgamma::P(): order " + out.str() + " not implemented");
2262  }
2263 }
2264 
2265 double Bsgamma::N_27()
2266 {
2267  double mcnorm = 1.131; // value fixed according to arXiv:1003.5012, in order to employ the remaining corrections given in that work
2268  double lambda2 = mu_G2/3.;
2269 
2270  return -1./18. * (C7_0 * ( 2.*C2_0 - C1_0/3. )).real() * lambda2/mcnorm/mcnorm;
2271 }
2272 
2273 double Bsgamma::N_77(double E0, double mu)
2274 {
2275  double z = 1. - delta(E0);
2276  double z2 = z*z;
2277  double z3 = z2*z;
2278  double z4 = z3*z;
2279  double umz2 = (1.-z)*(1.-z);
2280  double Lz = log(1. - z);
2281  double Lz2 = Lz*Lz;
2282  double Lb = 2. * log(mu/Mb_kin);
2283 
2284  double corrLambda2_rad;
2285  double corrLambda2_sem;
2286  double corrLambda2_mix;
2287  double corrLambda2;
2288  double corrLambda3;
2289 
2290  double alsb = SM.Alstilde5(Mb_kin);
2291  double Lambda_pert = 64./9. * alsb * mu_kin *
2292  (1. + 4. * alsb * (9./2. * (log(Mb_kin/2./mu_kin) + 8./3.)
2293  - 3. * (M_PI*M_PI/6. - 13./12.)) );
2294  double mu_pi2_pert = 3./4. * mu_kin * Lambda_pert - 48. * alsb*alsb * mu_kin*mu_kin;
2295  double rho_D3_pert = 1./2. * mu_kin*mu_kin * Lambda_pert - 128./3. * alsb*alsb * mu_kin*mu_kin*mu_kin;
2296 
2297  double lambda1 = -mu_pi2 + mu_pi2_pert;
2298  double lambda2 = mu_G2/3.;
2299  double rho1 = rho_D3 - rho_D3_pert;
2300 
2301  double f1EGN = 16./9. * ( 4. - M_PI*M_PI) - 8./3. * Lz2 -
2302  ( 4. * z * ( 30. - 63. * z + 31. * z2 + 5. * z3))/(9. * umz2) -
2303  ( 4. * (30. - 72. * z + 51. * z2 - 2. * z3 - 3. * z4))/(9. * umz2) * Lz;
2304  double f2EGN = -2./9. * ( 87. + 32. * M_PI*M_PI) - 32./3. * Lz2 +
2305  2. * ( 162. - 244. * z + 113. * z2 - 7. * z3)/(3. * (1. - z)) * Lz +
2306  2. * z * ( 54 - 49. * z + 15. * z2)/(1. - z);
2307 
2308  corrLambda2_rad = lambda1 * ( f1EGN/8. - 4./3. * (Lb + 1.) )
2309  + lambda2 * (f2EGN/8. + 12. * (Lb + 1.) );
2310  corrLambda2_sem = -3. * 4.98 * lambda2 + (25. - 4. * M_PI*M_PI)/12.*lambda1;
2311  corrLambda2_mix = 1./8. * (9. * lambda2 - lambda1) * Kij_1(7,7,E0,mu).real();
2312 
2313  corrLambda2 = corrLambda2_rad - corrLambda2_sem + corrLambda2_mix;
2314 
2315  corrLambda3 = (-88./6. + 16.*log(2.))* rho1 /Mb_kin/Mb_kin/Mb_kin;
2316 
2317  return (C7_0.abs2() + C7p_0.abs2()) * (4. * Alstilde / Mb_kin / Mb_kin * corrLambda2 + corrLambda3);
2318 }
2319 
2320 double Bsgamma::N(double E0, double mu)
2321 {
2322  return N_27() + N_77(E0,mu) + BLNPcorr * P0(E0);
2323 }
2324 
2325 double Bsgamma::C_sem()
2326 {
2327  double z=zeta();
2328  return (1. - 8. * z + 8. * z*z*z - z*z*z*z - 12. * z*z * log(z)) * ( 0.903
2329  - 0.588 * (SM.Alstilde5(4.6)*4*M_PI - 0.22) + 0.0650 * (Mb_kin - 4.55)
2330  - 0.1080 * (Mc - 1.05) - 0.0122 * mu_G2 - 0.199 * rho_D3 + 0.004 * rho_LS3);
2331 }
2332 
2333 void Bsgamma::updateParameters()
2334 {
2335  mu_kin=SM.getOptionalParameter("mukin");
2336  BRsl=SM.getOptionalParameter("BRsem")/100.;
2337  Mb_kin=SM.getOptionalParameter("Mbkin");
2338  Mc=SM.getOptionalParameter("Mcatmuc");
2339  mu_pi2=SM.getOptionalParameter("mupi2");
2340  rho_D3=SM.getOptionalParameter("rhoD3");
2341  mu_G2=SM.getOptionalParameter("muG2");
2342  rho_LS3=SM.getOptionalParameter("rhoLS3");
2343  mu_b=SM.getOptionalParameter("mu_b_bsgamma");
2344  mu_c=SM.getOptionalParameter("mu_c_bsgamma");
2345  C=C_sem();
2346 
2347  ale=SM.getAle();
2348  alsUps=8./M_PI * mu_kin/Mb_kin * ( 1. + 3./8. * mu_kin/Mb_kin );
2349  Alstilde = SM.Alstilde5(mu_b);
2350  AleatMztilde=SM.ale_OS(SM.getMz())/4./M_PI;
2351  Ms=SM.getQuarks(QCD::STRANGE).getMass();
2352  Mb=SM.getQuarks(QCD::BOTTOM).getMass();
2353  Mz=SM.getMz();
2354  V_ub=SM.getCKM().getV_ub();
2355  V_cb=SM.getCKM().getV_cb();
2356  V_tb=SM.getCKM().getV_tb();
2357 
2358  if(WET_NP_btos){
2359  C_7_NP = SM.getOptionalParameter("C7_NP");
2360  C_7p_NP = SM.getOptionalParameter("C7p_NP");
2361  }
2362  else if(SMEFT_NP_btos){
2363  gslpp::complex SMEFT_factor = (M_PI/SM.getAle())*(SM.v()*1.e-3)*(SM.v()*1.e-3)/SM.getCKM().computelamt_s();
2364  C_7_NP = SM.getOptionalParameter("CdB")-SM.getOptionalParameter("CdW");
2365  C_7_NP *= SMEFT_factor*SM.getAle()*8.*M_PI*SM.v()/Mb;
2366  C_7p_NP = SM.getOptionalParameter("CpdB")-SM.getOptionalParameter("CpdW");
2367  C_7p_NP *= SMEFT_factor*SM.getAle()*8.*M_PI*SM.v()/Mb;
2368  }
2369  else{
2370  C_7_NP = 0.;
2371  C_7p_NP = 0.;
2372  }
2373 
2374  if (SUM) {
2375  CKMratio=(V_tb/V_cb).abs2()*(1. - V_tb.abs2());
2376  CKMu=-V_ub.abs2()/(1. - V_tb.abs2());
2377  CKMusq = (V_ub/V_tb).abs2() * (1. - V_ub.abs2())/(1. - V_tb.abs2());
2378  }
2379  else
2380  switch (quark) {
2381  case StandardModel::STRANGE:
2382  CKMratio=(SM.getCKM().computelamt_s()/V_cb).abs2();
2383  CKMu=SM.getCKM().computelamu_s().conjugate() / SM.getCKM().computelamt_s().conjugate(); // -0.00802793 + 0.0180942*gslpp::complex::i(); //
2384  CKMusq = CKMu.abs2();
2385  break;
2386  case StandardModel::DOWN:
2387  CKMratio=(SM.getCKM().computelamt_d()/V_cb).abs2();
2388  CKMu=SM.getCKM().computelamu_d().conjugate() / SM.getCKM().computelamt_d().conjugate(); // 0.00745398 - 0.40416*gslpp::complex::i(); //
2389  CKMusq = CKMu.abs2();
2390  break;
2391  default:
2392  std::stringstream out;
2393  out << quark;
2394  throw std::runtime_error("bqgamma: quark " + out.str() + " not implemented");
2395  }
2396 
2397  BLNPcorr=SM.getOptionalParameter("BLNPcorr");
2398 
2399  checkCache();
2400 
2401  if (Intb_updated == 0) {
2402  Intb1Cached = 0;
2403  Intb2Cached = 0;
2404  Intb3Cached = 0;
2405  Intb4Cached = 0;
2406  Intbb1Cached = 0;
2407  Intbb2Cached = 0;
2408  Intbb4Cached = 0;
2409  IntPhi772rCached = 0;
2410  }
2411  if (Intbc_updated == 0) {
2412  Intbc1Cached = 0;
2413  Intbc2Cached = 0;
2414  Intc1Cached = 0;
2415  Intc1imCached = 0;
2416  Intc2Cached = 0;
2417  Intc3Cached = 0;
2418  IntccCached = 0;
2419  Intcc1Cached = 0;
2420  Intcc1p1Cached = 0;
2421  }
2422 
2423  computeCoeff(mu_b);
2424 
2425  overall = BRsl * CKMratio * 6. * ale / (M_PI * C);
2426 
2427 }
2428 
2429 double Bsgamma::computeThValue()
2430 {
2431  double E0 = getBinMin();
2432 
2433  updateParameters();
2434 
2435  if (obs == 1)
2436  return overall * ( P(E0, mu_b, mu_c, NNLO) + N(E0, mu_b) );
2437  if (obs == 2)
2438  return (Alstilde * P21_CPodd(E0, mu_b) + Vub_NLO_CPodd(E0) ) / (P(E0, mu_b, mu_c, NNLO) + N(E0, mu_b) );
2439 
2440  throw std::runtime_error("Bsgamma::computeThValue(): Observable type not defined. Can be only 1 or 2");
2441 }
Bsgamma::BRsl
double BRsl
Definition: bsgamma.h:1757
Bsgamma::getKc_re_t_1mt2
double getKc_re_t_1mt2(double t)
The function .
Definition: bsgamma.h:395
Bsgamma::getKc_im_Kb_1mt2
double getKc_im_Kb_1mt2(double t)
The function .
Definition: bsgamma.h:637
Bsgamma::Intcc1Cached
unsigned int Intcc1Cached
Definition: bsgamma.h:1827
Bsgamma::allcoeff
gslpp::vector< gslpp::complex > ** allcoeff
Definition: bsgamma.h:1777
std_make_vector.h
Bsgamma::C7_2
gslpp::complex C7_2
Definition: bsgamma.h:1800
Bsgamma::a
gslpp::complex a(double z)
The funcion as defined in .
Definition: bsgamma.cpp:233
Bsgamma::overall
double overall
Definition: bsgamma.h:1765
Bsgamma::Int_b4
double Int_b4(double E0)
Integral of the functions getKb_re_t2_1mt() and getKb_re_t2_1mt2().
Definition: bsgamma.cpp:419
Bsgamma::Int_b3
double Int_b3(double E0)
Integral of the functions getKb_re_t() and getKb_re_t_1mt().
Definition: bsgamma.cpp:399
Bsgamma::Int_c1
gslpp::complex Int_c1(double E0)
Integral of the functions getKc_re_1mt(), getKc_im_1mt() and getKc_re_1mt2(), getKc_im_1mt2().
Definition: bsgamma.cpp:555
Bsgamma::Phi12_1
double Phi12_1(double E0)
The function from .
Definition: bsgamma.cpp:774
Bsgamma::Phi34_1
double Phi34_1(double E0)
The function obtained using the prescription of .
Definition: bsgamma.cpp:920
Bsgamma::Intbc_updated
unsigned int Intbc_updated
Definition: bsgamma.h:1849
Bsgamma::Phi26_1_4body
double Phi26_1_4body(double E0)
The function obtained from .
Definition: bsgamma.cpp:856
Bsgamma::C2_1
gslpp::complex C2_1
Definition: bsgamma.h:1790
Bsgamma::getKc_im_t_1mt
double getKc_im_t_1mt(double t)
The function .
Definition: bsgamma.h:384
Bsgamma::P32
double P32(double E0, double mu)
The perturbative part of the BR as defined in .
Definition: bsgamma.cpp:2062
Bsgamma::Y2NV_PHI3
double Y2NV_PHI3(double rho)
The function from arXiv:0805.3911v2.
Definition: bsgamma.cpp:1384
Bsgamma::getKc_re_t_1mt
double getKc_re_t_1mt(double t)
The function .
Definition: bsgamma.h:373
Bsgamma::CacheIntc3
gslpp::complex CacheIntc3
Definition: bsgamma.h:1842
StandardModel::v
virtual double v() const
The Higgs vacuum expectation value.
Definition: StandardModel.cpp:943
Bsgamma::getKb_re_t2_1mt2
double getKb_re_t2_1mt2(double t)
The function .
Definition: bsgamma.h:560
Bsgamma::CacheIntb3
double CacheIntb3
Definition: bsgamma.h:1833
Bsgamma::Int_Phi77_2rem
double Int_Phi77_2rem(double E0)
The integral of omega77()
Definition: bsgamma.cpp:1799
CKM::computelamu_d
gslpp::complex computelamu_d() const
The product of the CKM elements .
Definition: CKM.cpp:130
Bsgamma::getKc_im_Kb_t_1mt
double getKc_im_Kb_t_1mt(double t)
The function .
Definition: bsgamma.h:659
Bsgamma::SMEFT_NP_btos
bool SMEFT_NP_btos
Definition: bsgamma.h:1743
Bsgamma::T2
double T2(double E0, double t)
The cutoff energy function as defined in .
Definition: bsgamma.cpp:143
QCD::BOTTOM
Definition: QCD.h:329
Bsgamma::Vub_NLO_2body_CPodd
double Vub_NLO_2body_CPodd()
The CP odd part of the 2 body NLO Vub part of the as defined in , .
Definition: bsgamma.cpp:2124
Bsgamma::kappa
gslpp::complex kappa(double Mq, double t)
The function as defined in .
Definition: bsgamma.cpp:353
Li3
cd Li3(cd x)
Definition: hpl.h:1016
Li2
cd Li2(cd x)
Definition: hpl.h:1011
Bsgamma::F_1
double F_1(double z)
The interpolated function from arXiv:1503.01791.
Definition: bsgamma.cpp:1624
Bsgamma::Phi78_1
double Phi78_1(double E0)
The function from .
Definition: bsgamma.cpp:1067
CKM::getV_tb
gslpp::complex getV_tb() const
A member for returning the value of the CKM element .
Definition: CKM.h:264
Bsgamma::Phi58_1
double Phi58_1(double E0)
The function obtained using the prescription of and adding the 4-body contribution from .
Definition: bsgamma.cpp:1023
ThObservable::setParametersForObservable
void setParametersForObservable(std::vector< std::string > parametersForObservable_i)
A set method to get the parameters for the specific observable.
Definition: ThObservable.h:109
Bsgamma::getKc_im_1mt2
double getKc_im_1mt2(double t)
The function .
Definition: bsgamma.h:450
Bsgamma::Intb_cache
double Intb_cache
Definition: bsgamma.h:1851
Bsgamma::BLNPcorr
double BLNPcorr
Definition: bsgamma.h:1770
Bsgamma::Intc2Cached
unsigned int Intc2Cached
Definition: bsgamma.h:1824
Bsgamma::getKb_abs2_t2_1mt
double getKb_abs2_t2_1mt(double t)
The function .
Definition: bsgamma.h:505
Bsgamma::Phi24_1
gslpp::complex Phi24_1(double E0)
The function obtained using the prescription of and adding the 4-body contribution from .
Definition: bsgamma.cpp:836
Bsgamma::getKb_abs2_t_1mt2
double getKb_abs2_t_1mt2(double t)
The function .
Definition: bsgamma.h:494
Bsgamma::CacheIntcc1p1
double CacheIntcc1p1
Definition: bsgamma.h:1845
Bsgamma::getKc_im_Kb_1mt
double getKc_im_Kb_1mt(double t)
The function .
Definition: bsgamma.h:615
Bsgamma::C6_0
gslpp::complex C6_0
Definition: bsgamma.h:1785
Bsgamma::Phi24_1_4body
double Phi24_1_4body(double E0)
The function obtained from .
Definition: bsgamma.cpp:828
Bsgamma::P0_4body
double P0_4body(double E0, double t)
The 4-body LO contribution as defined in .
Definition: bsgamma.cpp:172
Bsgamma::Phi45_1
double Phi45_1(double E0)
The function obtained using the prescription of .
Definition: bsgamma.cpp:968
Bsgamma::Rer22
double Rer22(double z)
The function from .
Definition: bsgamma.cpp:1150
Bsgamma::Phi66_1
gslpp::complex Phi66_1(double E0)
The function obtained using the prescription of .
Definition: bsgamma.cpp:1029
Bsgamma::Vub_NLO_3body_A
double Vub_NLO_3body_A(double E0)
The first piece of the 3 body NLO Vub part of the , .
Definition: bsgamma.cpp:2134
Bsgamma::allcoeffprime
gslpp::vector< gslpp::complex > ** allcoeffprime
Definition: bsgamma.h:1778
Bsgamma::Gamma_t
gslpp::complex Gamma_t(double t)
The function as defined in .
Definition: bsgamma.cpp:345
Bsgamma::Phi56_1
gslpp::complex Phi56_1(double E0)
The function obtained using the prescription of and adding the 4-body contribution from .
Definition: bsgamma.cpp:1002
Bsgamma::SUM
bool SUM
Definition: bsgamma.h:1740
Bsgamma::getKc_abs2_t_1mt
double getKc_abs2_t_1mt(double t)
The function .
Definition: bsgamma.h:329
Bsgamma::Int_c2
gslpp::complex Int_c2(double E0)
Integral of the functions getKc_re_t_1mt(), getKc_im_t_1mt() and getKc_re_t_1mt2(),...
Definition: bsgamma.cpp:583
Bsgamma::Phi57_1
double Phi57_1(double E0)
The function obtained using the prescription of and adding the 4-body contribution from .
Definition: bsgamma.cpp:1015
Bsgamma::getKb_abs2_1mt2
double getKb_abs2_1mt2(double t)
The function .
Definition: bsgamma.h:472
Bsgamma::Phi77_2beta0
double Phi77_2beta0(double E0, double mu)
The function from ..
Definition: bsgamma.cpp:1228
Bsgamma::Mz
double Mz
Definition: bsgamma.h:1756
Bsgamma::Phi25_1_4body
double Phi25_1_4body(double E0)
The function obtained from .
Definition: bsgamma.cpp:842
Bsgamma::errINT
double errINT
Definition: bsgamma.h:1811
Bsgamma::P22
double P22(double E0, double mu_b, double mu_c)
The perturbative part of the BR as defined in .
Definition: bsgamma.cpp:2038
Bsgamma::P0
double P0(double E0)
The perturbative part of the BR as defined in .
Definition: bsgamma.cpp:1969
Bsgamma::Y2
double Y2(double E0, double mu)
The function from arXiv:0805.3911v2 and arXiv:1005.5587v1.
Definition: bsgamma.cpp:1441
Bsgamma::Mc
double Mc
Definition: bsgamma.h:1753
Bsgamma::V_cb
gslpp::complex V_cb
Definition: bsgamma.h:1761
CKM::computelamt_s
gslpp::complex computelamt_s() const
The product of the CKM elements .
Definition: CKM.cpp:136
Bsgamma::Y2NL
double Y2NL(double E0, double mu)
The function from arXiv:0805.3911v2.
Definition: bsgamma.cpp:1342
Bsgamma::Intbb4Cached
unsigned int Intbb4Cached
Definition: bsgamma.h:1819
Bsgamma::Phi13_1
gslpp::complex Phi13_1(double E0)
The function obtained using the prescription of .
Definition: bsgamma.cpp:779
Bsgamma::C3_1
gslpp::complex C3_1
Definition: bsgamma.h:1791
Bsgamma::CacheIntPhi772r
double CacheIntPhi772r
Definition: bsgamma.h:1846
Bsgamma::mddel_f_NLO
double mddel_f_NLO(double z, double E0)
The function from arXiv:1503.01791.
Definition: bsgamma.cpp:1555
Bsgamma::CKMusq
double CKMusq
Definition: bsgamma.h:1764
Bsgamma::updateParameters
void updateParameters()
The update parameter method for bsgamma.
Definition: bsgamma.cpp:2333
Bsgamma::delddel_Phi28_1
double delddel_Phi28_1(double z, double E0)
Derivative of the function Phi28_1() used to compute effects of massive quark loops on gluon lines.
Definition: bsgamma.cpp:1650
Bsgamma::zdz_Phi22_1
double zdz_Phi22_1(double E0)
Derivative of the function Phi22_1() used to compute effects of massive quark loops on gluon lines.
Definition: bsgamma.cpp:1645
Bsgamma::C1_1
gslpp::complex C1_1
Definition: bsgamma.h:1789
make_vector
Definition: std_make_vector.h:15
LO
Definition: OrderScheme.h:33
Bsgamma::V_tb
gslpp::complex V_tb
Definition: bsgamma.h:1762
Bsgamma::f_b
double f_b(double z)
The function from arXiv:1503.01791.
Definition: bsgamma.cpp:1614
Bsgamma::Y1
double Y1(double E0, double mu)
The function from arXiv:0805.3911v2.
Definition: bsgamma.cpp:1282
Polylogarithms::Li3
double Li3(const double x) const
The trilogarithm .
Definition: Polylogarithms.cpp:36
Bsgamma::CacheIntb2
double CacheIntb2
Definition: bsgamma.h:1832
Bsgamma::C4_1
gslpp::complex C4_1
Definition: bsgamma.h:1792
Bsgamma::K77_2_z1
double K77_2_z1(double E0, double mu)
The function computed in the limit .
Definition: bsgamma.cpp:1826
Bsgamma::CacheIntcc1
double CacheIntcc1
Definition: bsgamma.h:1844
StandardModel.h
Bsgamma::f_c
double f_c(double z)
The function from arXiv:1503.01791.
Definition: bsgamma.cpp:1619
Bsgamma::CacheIntb4
double CacheIntb4
Definition: bsgamma.h:1834
Bsgamma::delta
double delta(double E0)
The cutoff energy function .
Definition: bsgamma.cpp:105
Bsgamma::Phi68_1
gslpp::complex Phi68_1(double E0)
The function obtained using the prescription of and adding the 4-body contribution from .
Definition: bsgamma.cpp:1052
Bsgamma::CacheIntc2
gslpp::complex CacheIntc2
Definition: bsgamma.h:1841
gslpp::complex
A class for defining operations on and functions of complex numbers.
Definition: gslpp_complex.h:35
Bsgamma::r1
gslpp::complex r1(int i, double z)
The funcion as defined in .
Definition: bsgamma.cpp:286
Bsgamma::mu_b
double mu_b
Definition: bsgamma.h:1749
CKM::getV_cb
gslpp::complex getV_cb() const
A member for returning the value of the CKM element .
Definition: CKM.h:237
gslpp::log
complex log(const complex &z)
Definition: gslpp_complex.cpp:342
Bsgamma::Phi27_1
gslpp::complex Phi27_1(double E0, double z)
The function from .
Definition: bsgamma.cpp:870
Bsgamma::C_7_NP
gslpp::complex C_7_NP
Definition: bsgamma.h:1805
Bsgamma::getKc_im_Kb_t_1mt2
double getKc_im_Kb_t_1mt2(double t)
The function .
Definition: bsgamma.h:681
Bsgamma::N_27
double N_27()
The non perturbative part of the due to interference as defined in , .
Definition: bsgamma.cpp:2265
Bsgamma::Vub_NLO_CPodd
double Vub_NLO_CPodd(double E0)
The CP odd part of the total NLO Vub part of the , .
Definition: bsgamma.cpp:2218
Bsgamma::CacheIntbb2
double CacheIntbb2
Definition: bsgamma.h:1836
Bsgamma::INT
gsl_function INT
Definition: bsgamma.h:1808
Bsgamma::Int_cc1
double Int_cc1(double E0)
Integral of the functions getKc_abs2_1mt() and getKc_abs2_1mt^().
Definition: bsgamma.cpp:659
Bsgamma::C5_1
gslpp::complex C5_1
Definition: bsgamma.h:1793
gslpp::complex::abs2
double abs2() const
Definition: gslpp_complex.cpp:86
Bsgamma::Delta
double Delta(double r)
The function from Z. Phys. C 48, 673 (1990).
Definition: bsgamma.cpp:1713
StandardModel
A model class for the Standard Model.
Definition: StandardModel.h:477
Bsgamma::Phi28_1
gslpp::complex Phi28_1(double E0, double z)
The function from .
Definition: bsgamma.cpp:903
Bsgamma::getKc_re_1mt
double getKc_re_1mt(double t)
The function .
Definition: bsgamma.h:417
Bsgamma::Intbc1Cached
unsigned int Intbc1Cached
Definition: bsgamma.h:1820
Bsgamma::WET_NP_btos
bool WET_NP_btos
Definition: bsgamma.h:1743
Bsgamma::omega
double omega(double E0)
The cutoff energy function as defined in .
Definition: bsgamma.cpp:118
Bsgamma::getKb_abs2_1mt
double getKb_abs2_1mt(double t)
The function .
Definition: bsgamma.h:461
Bsgamma::b
gslpp::complex b(double z)
The funcion as defined in .
Definition: bsgamma.cpp:262
Bsgamma::Int_cc
double Int_cc(double E0)
Integral of the functions getKc_abs2_t() and getKc_abs2_t_1mt().
Definition: bsgamma.cpp:639
Bsgamma::Phi22_1
double Phi22_1(double E0)
The function from .
Definition: bsgamma.cpp:809
Bsgamma::Phi26_1
gslpp::complex Phi26_1(double E0)
The function obtained using the prescription of and adding the 4-body contribution from .
Definition: bsgamma.cpp:864
Bsgamma::Phi48_1
double Phi48_1(double E0)
The function obtained using the prescription of and adding the 4-body contribution from .
Definition: bsgamma.cpp:984
Bsgamma::Phi88_1
double Phi88_1(double E0)
The function from .
Definition: bsgamma.cpp:1080
Bsgamma::CacheIntcc
double CacheIntcc
Definition: bsgamma.h:1843
Bsgamma::Intbc_cache
gslpp::vector< double > Intbc_cache
Definition: bsgamma.h:1852
Bsgamma::alsUps
double alsUps
Definition: bsgamma.h:1747
StandardModel::Alstilde5
double Alstilde5(const double mu) const
The value of at any scale with the number of flavours and full EW corrections.
Definition: StandardModel.cpp:899
Bsgamma::getKb_re_t_1mt
double getKb_re_t_1mt(double t)
The function .
Definition: bsgamma.h:538
Bsgamma::P21
double P21(double E0, double mu)
The perturbative part of the BR as defined in .
Definition: bsgamma.cpp:1980
Bsgamma::dY1
double dY1(double E0)
The function from arXiv:0805.3911v2.
Definition: bsgamma.cpp:1274
QCD::Beta0
double Beta0(const double nf) const
The coefficient for a certain number of flavours .
Definition: QCD.cpp:466
lambda1
An observable class for the quartic Higgs potential coupling .
Definition: THDMquantities.h:382
bsgamma.h
Bsgamma::getKc_im_t
double getKc_im_t(double t)
The function .
Definition: bsgamma.h:362
Bsgamma::computeCoeff
void computeCoeff(double mu)
Compute the Wilson Coefficient.
Definition: bsgamma.cpp:1901
Bsgamma::r1_ew
gslpp::complex r1_ew(int i, double z)
The funcion as defined in .
Definition: bsgamma.cpp:313
CKM::computelamu_s
gslpp::complex computelamu_s() const
The product of the CKM elements .
Definition: CKM.cpp:146
Bsgamma::Kij_1
gslpp::complex Kij_1(int i, int j, double E0, double mu)
The function from .
Definition: bsgamma.cpp:1094
Bsgamma::Y2NV
double Y2NV(double E0, double mu)
The function from arXiv:0805.3911v2.
Definition: bsgamma.cpp:1407
Bsgamma::Intb3Cached
unsigned int Intb3Cached
Definition: bsgamma.h:1815
ClausenFunctions::Cl3
double Cl3(const double phi) const
The Clausen function of index 3, .
Definition: ClausenFunctions.cpp:26
Bsgamma::CacheIntbb1
double CacheIntbb1
Definition: bsgamma.h:1835
Bsgamma::Vub_NLO_2body
double Vub_NLO_2body()
The 2 body NLO Vub part of the as defined in , .
Definition: bsgamma.cpp:2114
Bsgamma::Y2NV_PHI1
double Y2NV_PHI1(double rho)
The function from arXiv:0805.3911v2.
Definition: bsgamma.cpp:1362
Bsgamma::Phi28_2beta0
double Phi28_2beta0(double E0, double mu)
The function from arXiv:1009.5685.
Definition: bsgamma.cpp:1202
Bsgamma::Kij_2
double Kij_2(int i, int j, double E0, double mu_b, double mu_c)
The function from arXiv:1503.01791.
Definition: bsgamma.cpp:1839
gslpp::complex::imag
const double & imag() const
Definition: gslpp_complex.cpp:59
Bsgamma::C5_0
gslpp::complex C5_0
Definition: bsgamma.h:1784
Bsgamma::AleatMztilde
double AleatMztilde
Definition: bsgamma.h:1746
Bsgamma::Y2NH
double Y2NH(double E0, double mu)
The function from arXiv:0805.3911v2.
Definition: bsgamma.cpp:1430
Bsgamma::Intbb1Cached
unsigned int Intbb1Cached
Definition: bsgamma.h:1817
Bsgamma::getKc_abs2_1mt2
double getKc_abs2_1mt2(double t)
The function .
Definition: bsgamma.h:340
StandardModel::ale_OS
double ale_OS(const double mu, orders order=FULLNLO) const
The running electromagnetic coupling in the on-shell scheme.
Definition: StandardModel.cpp:533
Bsgamma::quark
QCD::quark quark
Definition: bsgamma.h:1738
Bsgamma::C7_0
gslpp::complex C7_0
Definition: bsgamma.h:1786
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
Bsgamma::C_sem
double C_sem()
The ratio as defined in , but with coefficients slightly modified due to different imput parameters...
Definition: bsgamma.cpp:2325
Bsgamma::Intcc1p1Cached
unsigned int Intcc1p1Cached
Definition: bsgamma.h:1828
Bsgamma::delddel_Phi22_1
double delddel_Phi22_1(double E0)
Derivative of the function Phi22_1() used to compute effects of massive quark loops on gluon lines.
Definition: bsgamma.cpp:1638
Bsgamma::CacheIntbb4
double CacheIntbb4
Definition: bsgamma.h:1837
Bsgamma::Intbb2Cached
unsigned int Intbb2Cached
Definition: bsgamma.h:1818
Bsgamma::h27_2
double h27_2(double z, double E0)
The function from arXiv:1009.5685 and arXiv:1503.01791.
Definition: bsgamma.cpp:1586
Bsgamma::Phi44_1
double Phi44_1(double E0)
The function obtained using the prescription of .
Definition: bsgamma.cpp:963
Bsgamma::Int_bc1
gslpp::complex Int_bc1(double E0)
Integral of the functions getKc_re_Kb_1mt(), getKc_im_Kb_1mt() and getKc_re_Kb_1mt2(),...
Definition: bsgamma.cpp:499
Bsgamma::getKc_abs2_1mt
double getKc_abs2_1mt(double t)
The function .
Definition: bsgamma.h:318
Bsgamma::IntPhi772rCached
unsigned int IntPhi772rCached
Definition: bsgamma.h:1829
Bsgamma::Intc1imCached
unsigned int Intc1imCached
Definition: bsgamma.h:1823
Bsgamma::mu_pi2
double mu_pi2
Definition: bsgamma.h:1766
Bsgamma::w_INT
gsl_integration_cquad_workspace * w_INT
Definition: bsgamma.h:1809
Bsgamma::P
double P(double E0, double mu_b, double mu_c, orders order)
The perturbative part of the as defined in , .
Definition: bsgamma.cpp:2233
Bsgamma::Int_bc2
gslpp::complex Int_bc2(double E0)
Integral of the functions getKc_re_Kb_t_1mt(), getKc_im_Kb_t_1mt() and getKc_re_Kb_t_1mt2(),...
Definition: bsgamma.cpp:527
Bsgamma::getKb_re_1mt2
double getKb_re_1mt2(double t)
The function .
Definition: bsgamma.h:593
gslpp::pow
complex pow(const complex &z1, const complex &z2)
Definition: gslpp_complex.cpp:395
Bsgamma::C2_0
gslpp::complex C2_0
Definition: bsgamma.h:1781
Bsgamma::ff7_sMP
double ff7_sMP(double E0)
The 4-body NLO correction due to and s, , from .
Definition: bsgamma.cpp:709
gslpp::sqrt
complex sqrt(const complex &z)
Definition: gslpp_complex.cpp:385
Bsgamma::C
double C
Definition: bsgamma.h:1758
gslpp::complex::i
static const complex & i()
Definition: gslpp_complex.cpp:154
Bsgamma::Int_cc1_part1
double Int_cc1_part1(double E0)
Integral of the functions getKc_abs2_1mt().
Definition: bsgamma.cpp:680
Bsgamma::rho
double rho(double E0)
The cutoff energy function as defined in .
Definition: bsgamma.cpp:110
Bsgamma::mu_kin
double mu_kin
Definition: bsgamma.h:1751
Bsgamma::Mb_kin
double Mb_kin
Definition: bsgamma.h:1752
Bsgamma::getKc_abs2_t
double getKc_abs2_t(double t)
The function .
Definition: bsgamma.h:307
Bsgamma::CacheIntbc1
gslpp::complex CacheIntbc1
Definition: bsgamma.h:1838
Bsgamma::zeta
double zeta()
The squared ratio between and , .
Definition: bsgamma.cpp:228
Bsgamma::ff7_dMP
double ff7_dMP(double E0)
The 4-body NLO correction due to and d, , from .
Definition: bsgamma.cpp:695
Bsgamma::omega77
double omega77(double z)
The function, linear combination of the functions , and from hep-ph/0505097.
Definition: bsgamma.cpp:1754
Bsgamma::Phi47_1
double Phi47_1(double E0)
The function from and adding the 4-body contribution from .
Definition: bsgamma.cpp:978
Bsgamma::C8_0
gslpp::complex C8_0
Definition: bsgamma.h:1787
Bsgamma::Intb2Cached
unsigned int Intb2Cached
Definition: bsgamma.h:1814
Bsgamma::getKb_abs2_t2_1mt2
double getKb_abs2_t2_1mt2(double t)
The function .
Definition: bsgamma.h:516
NNLO
Definition: OrderScheme.h:35
Bsgamma::Y2NV_PHI2
double Y2NV_PHI2(double rho)
The function from arXiv:0805.3911v2.
Definition: bsgamma.cpp:1373
CKM::computelamt_d
gslpp::complex computelamt_d() const
The product of the CKM elements .
Definition: CKM.cpp:120
lambda2
An observable class for the quartic Higgs potential coupling .
Definition: THDMquantities.h:405
Bsgamma::computeThValue
double computeThValue()
Computes the Branching Ratio for the decay.
Definition: bsgamma.cpp:2429
gslpp_function_adapter.h
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
Bsgamma::Phi77_2rem
double Phi77_2rem(double E0)
The part of the function with no dependance, as defined in .
Definition: bsgamma.cpp:1814
Bsgamma::Intbc2Cached
unsigned int Intbc2Cached
Definition: bsgamma.h:1821
Bsgamma::Phi38_1
double Phi38_1(double E0)
The function obtained using the prescription of and adding the 4-body contribution from .
Definition: bsgamma.cpp:958
Bsgamma::Phi22_2beta0
double Phi22_2beta0(double E0, double mu)
The function from arXiv:1009.5685.
Definition: bsgamma.cpp:1182
ThObservable::SM
const StandardModel & SM
A reference to an object of StandardMode class.
Definition: ThObservable.h:121
Bsgamma::Vub_NLO_4body_CPodd
double Vub_NLO_4body_CPodd(double E0)
The CP odd part of the 4 body NLO Vub part of the obtained from , .
Definition: bsgamma.cpp:2204
Bsgamma::IntccCached
unsigned int IntccCached
Definition: bsgamma.h:1826
Bsgamma::N
double N(double E0, double mu)
The non perturbative part of the as defined in , .
Definition: bsgamma.cpp:2320
Bsgamma::C7p_1
gslpp::complex C7p_1
Definition: bsgamma.h:1803
Bsgamma::getKc_im_1mt
double getKc_im_1mt(double t)
The function .
Definition: bsgamma.h:428
QCD::quark
quark
An enum type for quarks.
Definition: QCD.h:323
Bsgamma::P11
double P11()
The perturbative part of the BR as defined in .
Definition: bsgamma.cpp:1974
Bsgamma::Alstilde
double Alstilde
Definition: bsgamma.h:1748
Bsgamma::V_ub
gslpp::complex V_ub
Definition: bsgamma.h:1760
Bsgamma::Int_b1
double Int_b1(double E0)
Integral of the functions getKb_re_1mt() and getKb_re_1mt2().
Definition: bsgamma.cpp:359
Bsgamma::Vub_NLO_3body_B
double Vub_NLO_3body_B(double E0)
The second piece of the 3 body NLO Vub part of the , .
Definition: bsgamma.cpp:2148
StandardModel::getFlavour
const Flavour & getFlavour() const
Definition: StandardModel.h:1020
Bsgamma::Phi55_1
double Phi55_1(double E0)
The function obtained using the prescription of .
Definition: bsgamma.cpp:990
Bsgamma::Phi33_1
double Phi33_1(double E0)
The function obtained using the prescription of .
Definition: bsgamma.cpp:908
Bsgamma::Phi23_1
gslpp::complex Phi23_1(double E0)
The function obtained using the prescription of and adding the 4-body contribution from .
Definition: bsgamma.cpp:822
Bsgamma::Bsgamma
Bsgamma(const StandardModel &SM_i, QCD::quark quark_i, int obsFlag)
Constructor.
Definition: bsgamma.cpp:22
Bsgamma::Phi77_1
double Phi77_1(double E0)
The function from .
Definition: bsgamma.cpp:1058
Bsgamma::C_7p_NP
gslpp::complex C_7p_NP
Definition: bsgamma.h:1806
Bsgamma::T1
double T1(double E0, double t)
The cutoff energy function as defined in .
Definition: bsgamma.cpp:128
Bsgamma::C7_1
gslpp::complex C7_1
Definition: bsgamma.h:1795
Bsgamma::Phi35_1
double Phi35_1(double E0)
The function obtained using the prescription of .
Definition: bsgamma.cpp:925
Bsgamma::C1_0
gslpp::complex C1_0
Definition: bsgamma.h:1780
Bsgamma::T3
double T3(double E0, double t)
The cutoff energy function as defined in .
Definition: bsgamma.cpp:158
Bsgamma::avaINT
double avaINT
Definition: bsgamma.h:1810
Bsgamma::Phi46_1
gslpp::complex Phi46_1(double E0)
The function obtained using the prescription of and adding the 4-body contribution from .
Definition: bsgamma.cpp:973
Bsgamma::Phi67_1
gslpp::complex Phi67_1(double E0)
The function obtained using the prescription of and adding the 4-body contribution from .
Definition: bsgamma.cpp:1043
Bsgamma::FOUR_BODY
bool FOUR_BODY
Definition: bsgamma.h:1742
Bsgamma::Intb1Cached
unsigned int Intb1Cached
Definition: bsgamma.h:1813
Bsgamma::C8_1
gslpp::complex C8_1
Definition: bsgamma.h:1796
Bsgamma::Intc1Cached
unsigned int Intc1Cached
Definition: bsgamma.h:1822
Bsgamma::getKb_abs2_t_1mt
double getKb_abs2_t_1mt(double t)
The function .
Definition: bsgamma.h:483
orders
orders
An enum type for orders in QCD.
Definition: OrderScheme.h:31
Bsgamma::Phi37_1
double Phi37_1(double E0)
The function obtained using the prescription of and adding the 4-body contribution from .
Definition: bsgamma.cpp:950
Bsgamma::Int_c3
gslpp::complex Int_c3(double E0)
Integral of the functions getKc_re_t(), getKc_im_t() and getKc_re_t_1mt(), getKc_im_t_1mt().
Definition: bsgamma.cpp:611
Bsgamma::Int_bb4
double Int_bb4(double E0)
Integral of the functions getKb_abs2_t2_1mt() and getKb_abs2_t2_1mt2().
Definition: bsgamma.cpp:479
Bsgamma::Intc3Cached
unsigned int Intc3Cached
Definition: bsgamma.h:1825
Bsgamma::f
double f(double r)
The function from hep-ph/0611123.
Definition: bsgamma.cpp:1693
Bsgamma::ff8_dMP
double ff8_dMP(double E0)
The 4-body NLO correction due to and d, , from .
Definition: bsgamma.cpp:723
Bsgamma::C7_1ew
gslpp::complex C7_1ew
Definition: bsgamma.h:1798
StandardModel::getMz
double getMz() const
A get method to access the mass of the boson .
Definition: StandardModel.h:721
Bsgamma::f_NLO_1
double f_NLO_1(double z)
The function from arXiv:1503.01791.
Definition: bsgamma.cpp:1454
Bsgamma::C3_0
gslpp::complex C3_0
Definition: bsgamma.h:1782
Bsgamma::Vub_NLO_3body_A_CPodd
double Vub_NLO_3body_A_CPodd(double E0)
The CP odd part of the first piece of the 3 body NLO Vub part of the , .
Definition: bsgamma.cpp:2143
Bsgamma::Phi23_1_4body
double Phi23_1_4body(double E0)
The function obtained from .
Definition: bsgamma.cpp:814
Bsgamma::EW_NLO
double EW_NLO(double mu)
The NLO electroweak correction to the BR as defined in .
Definition: bsgamma.cpp:2081
ThObservable
A class for a model prediction of an observable.
Definition: ThObservable.h:25
Bsgamma::Mb
double Mb
Definition: bsgamma.h:1755
Bsgamma::Intb_updated
unsigned int Intb_updated
Definition: bsgamma.h:1848
Bsgamma::Phi14_1
gslpp::complex Phi14_1(double E0)
The function obtained using the prescription of .
Definition: bsgamma.cpp:784
Bsgamma::Intb4Cached
unsigned int Intb4Cached
Definition: bsgamma.h:1816
Bsgamma::F_2
double F_2(double z)
The interpolated function from arXiv:1503.01791.
Definition: bsgamma.cpp:1631
QCD::STRANGE
Definition: QCD.h:327
ClausenFunctions
A class for the Clausen functions.
Definition: ClausenFunctions.h:22
Cl2
double Cl2(double x)
Definition: hpl.h:1026
Bsgamma::CacheIntb1
double CacheIntb1
Definition: bsgamma.h:1831
Bsgamma::mu_c
double mu_c
Definition: bsgamma.h:1750
gslpp::complex::real
const double & real() const
Definition: gslpp_complex.cpp:53
Bsgamma::ale
double ale
Definition: bsgamma.h:1745
Bsgamma::CacheIntc1
gslpp::complex CacheIntc1
Definition: bsgamma.h:1840
Bsgamma::getKb_re_1mt
double getKb_re_1mt(double t)
The function .
Definition: bsgamma.h:582
CKM::getV_ub
gslpp::complex getV_ub() const
A member for returning the value of the CKM element .
Definition: CKM.h:210
Bsgamma::CacheIntbc2
gslpp::complex CacheIntbc2
Definition: bsgamma.h:1839
Bsgamma::getKc_re_Kb_1mt2
double getKc_re_Kb_1mt2(double t)
The function .
Definition: bsgamma.h:626
ThObservable::getBinMin
double getBinMin()
A get method to get the minimum value of the bin.
Definition: ThObservable.h:82
Bsgamma::getKc_re_t
double getKc_re_t(double t)
The function .
Definition: bsgamma.h:351
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
Flavour::ComputeCoeffsgamma
gslpp::vector< gslpp::complex > ** ComputeCoeffsgamma(double mu, bool noSM=false, schemes scheme=NDR) const
Computes the Wilson coefficient for the process .
Definition: Flavour.cpp:124
Bsgamma::Int_bb1
double Int_bb1(double E0)
Integral of the functions getKb_abs2_1mt() and getKb_abs2_1mt2().
Definition: bsgamma.cpp:439
Bsgamma::CKMratio
double CKMratio
Definition: bsgamma.h:1759
Bsgamma::Phi88_2beta0
double Phi88_2beta0(double E0, double mu)
The function from arXiv:1009.5685.
Definition: bsgamma.cpp:1250
Bsgamma::C4_0
gslpp::complex C4_0
Definition: bsgamma.h:1783
Bsgamma::getKb_re_t_1mt2
double getKb_re_t_1mt2(double t)
The function .
Definition: bsgamma.h:571
Bsgamma::getKc_re_Kb_1mt
double getKc_re_Kb_1mt(double t)
The function .
Definition: bsgamma.h:604
Bsgamma::obs
int obs
Definition: bsgamma.h:1772
Flavour::ComputeCoeffprimesgamma
gslpp::vector< gslpp::complex > ** ComputeCoeffprimesgamma(double mu, schemes scheme=NDR) const
Computes the chirality flipped Wilson coefficient for the process .
Definition: Flavour.cpp:129
Bsgamma::Y2CA
double Y2CA(double E0, double mu)
The function from arXiv:1005.5587v1.
Definition: bsgamma.cpp:1316
Bsgamma::getKc_re_Kb_t_1mt
double getKc_re_Kb_t_1mt(double t)
The function .
Definition: bsgamma.h:648
NLO
Definition: OrderScheme.h:34
Bsgamma::Phi18_1
gslpp::complex Phi18_1(double E0, double z)
The function from .
Definition: bsgamma.cpp:804
Bsgamma::Phi17_1
gslpp::complex Phi17_1(double E0, double z)
The function from .
Definition: bsgamma.cpp:799
Bsgamma::ff8_sMP
double ff8_sMP(double E0)
The 4-body NLO correction due to and s, , from .
Definition: bsgamma.cpp:746
Bsgamma::Phi25_1
gslpp::complex Phi25_1(double E0)
The function obtained using the prescription of and adding the 4-body contribution from .
Definition: bsgamma.cpp:850
convertToGslFunction
gsl_function convertToGslFunction(const F &f)
Definition: gslpp_function_adapter.h:24
Bsgamma::zdz_Phi28_1
double zdz_Phi28_1(double z, double E0)
Derivative of the function Phi28_1() used to compute effects of massive quark loops on gluon lines.
Definition: bsgamma.cpp:1664
Bsgamma::P21_CPodd
double P21_CPodd(double E0, double mu)
The CP odd part of the perturbative part of the BR as defined in .
Definition: bsgamma.cpp:2006
Bsgamma::Vub_NNLO
double Vub_NNLO(double E0)
The NNLO Vub part of the as defined in xxxxxxxxx, .
Definition: bsgamma.cpp:2223
Polylogarithms
A class for the polylogarithms.
Definition: Polylogarithms.h:24
Bsgamma::f_q
double f_q(double z, double E0)
The function from arXiv:1503.01791.
Definition: bsgamma.cpp:1609
Bsgamma::Y2CF
double Y2CF(double E0, double mu)
The function from arXiv:1005.5587v1.
Definition: bsgamma.cpp:1289
Bsgamma::Int_b2
double Int_b2(double E0)
Integral of the functions getKb_re_t_1mt() and getKb_re_t_1mt2().
Definition: bsgamma.cpp:379
Bsgamma::Phi16_1
gslpp::complex Phi16_1(double E0)
The function obtained using the prescription of .
Definition: bsgamma.cpp:794
Bsgamma::rho_D3
double rho_D3
Definition: bsgamma.h:1768
Bsgamma::N_77
double N_77(double E0, double mu)
The non perturbative part of the due to interference as defined in arXiv:0911.2175,...
Definition: bsgamma.cpp:2273
Bsgamma::getKc_im_t_1mt2
double getKc_im_t_1mt2(double t)
The function .
Definition: bsgamma.h:406
Bsgamma::Vub_NLO_4body
double Vub_NLO_4body(double E0)
The 4 body NLO Vub part of the obtained from , .
Definition: bsgamma.cpp:2174
Bsgamma::Phi15_1
gslpp::complex Phi15_1(double E0)
The function obtained using the prescription of .
Definition: bsgamma.cpp:789
Bsgamma::Ms
double Ms
Definition: bsgamma.h:1754
Bsgamma::zdz_f_NLO
double zdz_f_NLO(double z, double E0)
The function from arXiv:1503.01791.
Definition: bsgamma.cpp:1459
QCD::DOWN
Definition: QCD.h:325
Bsgamma::C7p_0
gslpp::complex C7p_0
Definition: bsgamma.h:1802
Bsgamma::rho_LS3
double rho_LS3
Definition: bsgamma.h:1769
Bsgamma::CKMu
gslpp::complex CKMu
Definition: bsgamma.h:1763
Bsgamma::delddel_Phi88_1
double delddel_Phi88_1(double E0)
Derivative of the function Phi88_1() used to compute effects of massive quark loops on gluon lines.
Definition: bsgamma.cpp:1683
Bsgamma::getKb_re_t2_1mt
double getKb_re_t2_1mt(double t)
The function .
Definition: bsgamma.h:549
StandardModel::getAle
double getAle() const
A get method to retrieve the fine-structure constant .
Definition: StandardModel.h:748
Bsgamma::checkCache
void checkCache()
The caching method for bsgamma.
Definition: bsgamma.cpp:87
Bsgamma::C6_1
gslpp::complex C6_1
Definition: bsgamma.h:1794
Bsgamma::getKc_re_1mt2
double getKc_re_1mt2(double t)
The function .
Definition: bsgamma.h:439
Bsgamma::Vub_NLO
double Vub_NLO(double E0)
The total NLO Vub part of the , .
Definition: bsgamma.cpp:2213
Bsgamma::Int_bb2
double Int_bb2(double E0)
Integral of the functions getKb_abs2_t_1mt() and getKb_abs2_t_1mt2().
Definition: bsgamma.cpp:459
Bsgamma::Phi11_1
double Phi11_1(double E0)
The function from .
Definition: bsgamma.cpp:769
Bsgamma::Phi36_1
gslpp::complex Phi36_1(double E0)
The function obtained using the prescription of and adding the 4-body contribution from .
Definition: bsgamma.cpp:937
Bsgamma::P12
double P12()
The perturbative part of the BR as defined in .
Definition: bsgamma.cpp:2032
Bsgamma::Y2NV_PHI4
double Y2NV_PHI4(double rho)
The function from arXiv:0805.3911v2.
Definition: bsgamma.cpp:1395
Bsgamma::Vub_NLO_3body_B_CPodd
double Vub_NLO_3body_B_CPodd(double E0)
The CP odd part of the second piece of the 3 body NLO Vub part of the , .
Definition: bsgamma.cpp:2161
StandardModel::getCKM
CKM getCKM() const
A get method to retrieve the member object of type CKM.
Definition: StandardModel.h:879
Bsgamma::EWflag
bool EWflag
Definition: bsgamma.h:1741
Bsgamma::f_u
double f_u(double r)
The function obtained after multiplying the fitted function of arXiv:0803.0960 for and subtracting...
Definition: bsgamma.cpp:1730
Bsgamma::getKc_re_Kb_t_1mt2
double getKc_re_Kb_t_1mt2(double t)
The function .
Definition: bsgamma.h:670
Bsgamma::getKb_re_t
double getKb_re_t(double t)
The function .
Definition: bsgamma.h:527
Bsgamma::mu_G2
double mu_G2
Definition: bsgamma.h:1767