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