Hyperiso 1.0.3
Modular flavour-physics calculations, Wilson coefficients and statistical inference
Loading...
Searching...
No Matches
BXsllDecay.cpp
Go to the documentation of this file.
1#include "BXsllDecay.h"
2#include "wcoef_ids.hpp"
3
5 cache.cc_res_mass = {3.096916, 3.68609, 3.77292, 4.039 , 4.153 , 4.421 };
6 cache.cc_res_br = {5.93e-2 , 7.7e-3 , 1.1e-5 , 1.4e-5, 1.0e-5 , 1.1e-5};
7 cache.cc_res_width_tot = {9.29e-5 , 3.04e-4, 2.73e-2, 8.0e-2, 1.03e-1, 6.2e-2};
8 cache.cc_res_width_had = {8.147e-5, 2.9746e-4, 2.36e-2, 5.2e-2, 7.8e-2 , 4.3e-2};
9
11
12 cache.alpha_em = (*p)(ParamId{ParameterType::SM, "EW", {1, 1}}, DataType::VALUE);
13 // LOG_INFO("Bsll alpha_em =", cache.alpha_em);
14 double m_c = (*p)(ParamId{ParameterType::SM, "MASS", 4}, DataType::VALUE);
15 cache.m_b_1S = (*p)(ParamId{ParameterType::SM, "QCD", {5, 3}}, DataType::VALUE);
16 complex_t V_tb = (*p)(ParamId{ParameterType::SM, "VCKM", {2, 2}}, DataType::VALUE);
17 complex_t V_ts = (*p)(ParamId{ParameterType::SM, "VCKM", {2, 1}}, DataType::VALUE);
18 complex_t V_cb = (*p)(ParamId{ParameterType::SM, "VCKM", {1, 2}}, DataType::VALUE);
19 complex_t V_cs = (*p)(ParamId{ParameterType::SM, "VCKM", {1, 1}}, DataType::VALUE);
20 cache.m_c_hat = m_c / cache.m_b_1S;
21 cache.m_D_hat = (*p)(ParamId{ParameterType::FLAVOR, "FMASS", 421}, DataType::VALUE) / cache.m_b_1S;
22 double mu_b = (*p)(ParamId{ParameterType::WILSON, "B_SCALE", 1}, DataType::VALUE);
23 cache.alpha_s_mu_b = (*iobs_qcdp)(AlphasConfig(mu_b, MassType::POLE, MassType::POLE));
24 cache.z = std::pow(cache.m_c_hat, 2);
25 cache.L_b = std::log(mu_b / cache.m_b_1S);
26 cache.L_b_5GeV = std::log(mu_b / 5.0);
27
28 cache.pref_dB0_ds = 1. + 3. * (*p)(ParamId{ParameterType::DECAY, "B_Xs", 6}, DataType::VALUE) * g_lambda(cache.z) / (2. * pow(cache.m_b_1S, 2) * f(cache.z)) - (*p)(ParamId{ParameterType::DECAY, "B_Xsll", 2}, DataType::VALUE) * g_rho(cache.z) / (6. * pow(cache.m_b_1S, 3) * f(cache.z));
29 cache.pref_dB_ds = (*p)(ParamId{ParameterType::DECAY, "B_Xs", 2}, DataType::VALUE) * std::pow(std::abs(V_tb * std::conj(V_ts) / V_cb), 2) / (std::pow(2 * PI / cache.alpha_em, 2) * f(cache.z) * kappa(cache.z));
30 cache.pref_A0_0 = 1 + 3 * (*p)(ParamId{ParameterType::DECAY, "B_Xs", 6}, DataType::VALUE) * g_lambda(cache.z) / (2 * std::pow(cache.m_b_1S, 2) * f(cache.z));
31 cache.pref_A0_1 = 4 * (*p)(ParamId{ParameterType::DECAY, "B_Xsll", 1}, DataType::VALUE) / (3 * std::pow(cache.m_b_1S, 2));
32 cache.pref_delta_mb2 = 3. * (*p)(ParamId{ParameterType::DECAY, "B_Xs", 6}, DataType::VALUE) / (2. * std::pow(cache.m_b_1S * std::abs(V_tb), 2));
33 cache.pref_delta_mb3 = -(*p)(ParamId{ParameterType::DECAY, "B_Xsll", 2}, DataType::VALUE) / (pow(cache.m_b_1S, 3) * std::pow(std::abs(V_tb), 2));
34 cache.pref_delta_mc2 = 8. * (*p)(ParamId{ParameterType::DECAY, "B_Xs", 6}, DataType::VALUE) / (9. * std::pow(m_c, 2)) * std::abs(std::conj(V_cs) * V_cb / (std::conj(V_ts) * std::pow(V_tb, 3)));
35 cache.pref_delta_brems = cache.alpha_s_mu_b / (4. * PI);
36 cache.pref_delta_em = cache.alpha_em / (4. * PI);
37
38 size_t ff_order {20};
39 fill_cache(BV::f_17, 0, 1, cache.F_17_lookup, cache.L_b, cache.z, ff_order);
40 fill_cache(BV::f_27, 0, 1, cache.F_27_lookup, cache.L_b, cache.z, ff_order);
41 fill_cache(BV::f_19_1S, 0, 1, cache.F_19_lookup, cache.L_b, cache.z, ff_order);
42 fill_cache(BV::f_29_1S, 0, 1, cache.F_29_lookup, cache.L_b, cache.z, ff_order);
43
44 auto bound_func = std::bind(&BXsllDecay::delta_bremB_base, &*this, std::placeholders::_1);
45 fill_cache(bound_func, 0, 1, cache.delta_brems_lookup);
46
48
49 // auto hatify = [this] (double q2) { return q2 / pow(cache.m_b_1S, 2); };
50
51 // printf("alpha = %.5e\n", cache.alpha_em);
52 // printf("m_b_1S = %.5e\n", cache.m_b_1S);
53 // printf("z = %.5e\n", cache.z);
54 // printf("m_D = %.5e\n", cache.m_D_hat * cache.m_b_1S);
55
56 // double s = 0.0454;
57 // double w = 0.5;
58 // printf("pref dB_ds = %.5e\n", cache.pref_dB_ds);
59 // printf("pref dB0_ds = %.5e\n", cache.pref_dB0_ds);
60 // printf("pref dB_mb2 = %.5e\n", cache.pref_dB_ds * cache.pref_delta_mb2 / (*p)(ParamId{ParameterType::DECAY, "B_Xs", 2}).real());
61 // printf("pref dB_mb3 = %.5e\n", cache.pref_dB_ds * cache.pref_delta_mb3 / (*p)(ParamId{ParameterType::DECAY, "B_Xs", 2}).real());
62 // printf("pref dB_mc2 = %.5e\n", cache.pref_dB_ds * cache.pref_delta_mc2 / (*p)(ParamId{ParameterType::DECAY, "B_Xs", 2}).real());
63 // printf("pref dB_brems = %.5e\n", cache.pref_dB_ds * cache.pref_delta_brems);
64 // printf("pref dB_em = %.5e\n", cache.pref_dB_ds * cache.pref_delta_em / (*p)(ParamId{ParameterType::DECAY, "B_Xs", 2}).real());
65
66 // printf("C7_new = %.4e + %.4e i\n", real(C7_new(s, false)), imag(C7_new(s, false)));
67 // printf("C9_new = %.4e + %.4e i\n", real(C9_new(s, false)), imag(C9_new(s, false)));
68 // printf("C10_new = %.4e + %.4e i\n", real(C10_new(s, false)), imag(C10_new(s, false)));
69
70 // printf("g(0) = %.4e + %.4e i\n", real(g(0, s)), imag(g(0, s)));
71 // printf("g(1) = %.4e + %.4e i\n", real(g(1, s)), imag(g(1, s)));
72 // printf("g(m_c) = %.4e + %.4e i\n", real(g_ld(cache.m_c_hat, s)), imag(g_ld(cache.m_c_hat, s)));
73
74 // printf("f(z) = %.4e\n", f(cache.z));
75 // printf("h(z) = %.4e\n", h(cache.z));
76 // printf("g_rho(z) = %.4e\n", g_rho(cache.z));
77 // printf("g_lambda(z) = %.4e\n", g_lambda(cache.z));
78 // printf("kappa(z) = %.4e\n", kappa(cache.z));
79 // printf("f_7(s) = %.4e\n", f_7(s));
80 // printf("f_9(s) = %.4e\n", f_9(s));
81 // printf("Gm1 = %.4e + %.4e i\n", real(Gm1(s / cache.z)), imag(Gm1(s / cache.z)));
82 // printf("G0 = %.4e + %.4e i\n", real(G0(s / cache.z)), imag(G0(s / cache.z)));
83 // printf("Delta_i_23 = %.4e + %.4e i\n", real(Delta_i_23(s, cache.z, w)), imag(Delta_i_23(s, cache.z, w)));
84 // printf("Delta_i_27 = %.4e + %.4e i\n", real(Delta_i_27(s, cache.z, w)), imag(Delta_i_27(s, cache.z, w)));
85 // printf("tau_22 = %.4e + %.4e i\n", real(tau_22(s, w, Delta_i_23(s, cache.z, w), Delta_i_27(s, cache.z, w))), imag(tau_22(s, w, Delta_i_23(s, cache.z, w), Delta_i_27(s, cache.z, w))));
86 // printf("tau_27 = %.4e + %.4e i\n", real(tau_27(s, w, Delta_i_23(s, cache.z, w), Delta_i_27(s, cache.z, w))), imag(tau_27(s, w, Delta_i_23(s, cache.z, w), Delta_i_27(s, cache.z, w))));
87 // printf("tau_28 = %.4e + %.4e i\n", real(tau_28(s, w, Delta_i_23(s, cache.z, w), Delta_i_27(s, cache.z, w))), imag(tau_28(s, w, Delta_i_23(s, cache.z, w), Delta_i_27(s, cache.z, w))));
88 // printf("tau_29 = %.4e + %.4e i\n", real(tau_29(s, w, Delta_i_23(s, cache.z, w), Delta_i_27(s, cache.z, w))), imag(tau_29(s, w, Delta_i_23(s, cache.z, w), Delta_i_27(s, cache.z, w))));
89 // printf("tau_210 = %.4e + %.4e i\n", real(tau_210(s, cache.z)), imag(tau_210(s, cache.z)));
90 // printf("tau_77(s) = %.4e\n", tau_77(s));
91 // printf("tau_78(s) = %.4e\n", tau_78(s));
92 // printf("tau_79(s) = %.4e\n", tau_79(s));
93 // printf("tau_710(s) = %.4e\n", tau_710(s));
94 // printf("tau_88(s) = %.4e\n", tau_88(s));
95 // printf("tau_89(s) = %.4e\n", tau_89(s));
96 // printf("tau_810(s) = %.4e\n", tau_810(s));
97 // printf("tau_99(s) = %.4e\n", tau_99(s));
98 // printf("tau_910(s) = %.4e\n", tau_910(s));
99 // printf("sigma(s) = %.4e\n", sigma(s));
100 // printf("sigma_9(s) = %.4e\n", sigma_9(s));
101 // printf("sigma_7(s) = %.4e\n", sigma_7(s, cache.L_b));
102 // printf("F(s/z) = %.4e + %.4e i\n", real(F(s / cache.z)), imag(F(s / cache.z)));
103
104 // printf("pref dB_mb2 (s = %.4f) = %.4e\n", s, cache.pref_dB_ds * cache.pref_delta_mb2);
105 // printf("pref dB_mb3 (s = %.4f) = %.4e\n", s, cache.pref_dB_ds * cache.pref_delta_mb3);
106 // printf("pref dB_mc2 (s = %.4f) = %.4e\n", s, cache.pref_dB_ds * cache.pref_delta_mc2);
107
108 // printf("dB0_ds (s = %.4f) = %.4e\n", s, cache.pref_dB_ds * cache.pref_dB0_ds * dB0_ds(s, cache.m_l_hat));
109 // printf("dB_mb2 (s = %.4f) = %.4e\n", s, cache.pref_dB_ds * cache.pref_delta_mb2 * delta_mb2(s));
110 // printf("dB_mb3 (s = %.4f) = %.4e\n", s, cache.pref_dB_ds * cache.pref_delta_mb3 * delta_mb3(s));
111 // printf("dB_mc2 (s = %.4f) = %.4e\n", s, cache.pref_dB_ds * cache.pref_delta_mc2 * delta_mc2(s));
112 // printf("dB_brems A (s = %.4f) = %.4e\n", s, cache.pref_dB_ds * cache.pref_delta_brems * (delta_bremA(s)));
113 // printf("dB_brems B (s = %.4f) = %.4e\n", s, cache.pref_dB_ds * cache.pref_delta_brems * (delta_bremB(s)));
114 // printf("dB_em (s = %.4f) = %.4e\n", s, cache.pref_dB_ds * cache.pref_delta_em * delta_em(s, cache.L_l));
115}
116
118 if (cfg.gen != gen) {
119 cfg.gen = gen;
121 }
122}
123
124// void BXsllDecay::fill_wilson_cache() {
125// cache.C = w_proxy->getAFR(WGroup::B, w_config.order);
126// auto C_P = w_proxy->getAFR(WGroup::BPrime, w_config.order);
127// auto C_S = w_proxy->getAFR(WGroup::BScalar, w_config.order);
128// cache.C.insert(C_P.begin(), C_P.end());
129// cache.C.insert(C_S.begin(), C_S.end());
130
131// cache.C_LO = w_proxy->getAR(WGroup::B, QCDOrder::LO);
132// auto C_P_LO = w_proxy->getAR(WGroup::BPrime, QCDOrder::LO);
133// cache.C_LO.insert(C_P_LO.begin(), C_P_LO.end());
134// }
135
137 cache.C.clear();
138 cache.C_LO.clear();
139
140 auto C_B = w_proxy->getAFR(WGroup::B, w_config.order);
141 auto C_P = w_proxy->getAFR(WGroup::BPrime, w_config.order);
142 auto C_S = w_proxy->getAFR(WGroup::BScalar, w_config.order);
143
144 for (const auto& [id, val] : C_B) {
145 cache.C[id] = val;
146 }
147 for (const auto& [id, val] : C_P) {
148 cache.C[id] = val;
149 }
150 for (const auto& [id, val] : C_S) {
151 cache.C[id] = val;
152 }
153
154 auto C_B_LO = w_proxy->getAR(WGroup::B, QCDOrder::LO);
155 auto C_P_LO = w_proxy->getAR(WGroup::BPrime, QCDOrder::LO);
156
157 for (const auto& [id, val] : C_B_LO) {
158 cache.C_LO[id] = val;
159 }
160 for (const auto& [id, val] : C_P_LO) {
161 cache.C_LO[id] = val;
162 }
163}
164
166 cache.m_l_hat = (*p)(ParamId{ParameterType::SM, "MASS", 11 + 2 * (int)cfg.gen}, DataType::VALUE) / cache.m_b_1S;
167 cache.L_l = -2 * std::log(cache.m_l_hat);
168
169 for (size_t i = 0; i < 3; i++)
170 cache.rand_err[i] = (*p)(ParamId{ParameterType::DECAY, "B_Xsll", {3, 11 + 2 * (int)cfg.gen, i}}, DataType::VALUE);
171}
172
173double BXsllDecay::f(double z) {
174 return 1 - 8 * z + 8 * pow(z, 3) - pow(z, 4) - 12 * pow(z, 2) * log(z);
175}
176
177double BXsllDecay::h(double z) {
178 double z2 = z * z;
179 double z3 = z2 * z;
180 double z4 = z3 * z;
181 double lz = log(z);
182 double lw = log(1 - z);
183 double lz2 = lz * lz;
184
185 return -(1 - z2) * (25./4 - 239./3 * z + 25./4 * z2) + z * lz * (20. + 90. * z - 4./3 * z2 + 17./3 * z3)
186 + z2 * lz2 * (36. + z2) + (1. - z2) * (17./3. - 64./3. * z + 17./3. * z2) * lw - 4. * (1. + 30. * z2 + z4) * lz * lw
187 - (1. + 16. * z2 + z4) * (6. * Li2(z) - PI2) - 32. * pow(z, 1.5) * (1. + z) * (PI2 - 4. * Li2(sqrt(z)) + 4. * Li2(-sqrt(z)) - 2. * lz * log((1.-sqrt(z))/(1.+sqrt(z))));
188}
189
190double BXsllDecay::g_rho(double z) {
191 return 77.-88.*z+24.*pow(z,2.)-8.*pow(z,3)+5.*pow(z,4)+48.*log(z)+36.*pow(z,2)*log(z);
192}
193
194double BXsllDecay::g_lambda(double z) {
195 return 3.-8.*z+24.*pow(z,2)-24.*pow(z,3)+5.*pow(z,4)+12.*pow(z,2)*log(z);
196}
197
198double BXsllDecay::kappa(double z) {
199 return 1 - 2. * cache.alpha_s_mu_b * h(z) / (3. * PI * f(z));
200}
201
202double BXsllDecay::f_7(double s) {
203 return 1./6./(s-1.)/(s-1.)*(24.*(1.+13.*s-4.*s*s)*Li2(sqrt(s))+12.*(1.-17.*s+6.*s*s)*Li2(s)+6.*s*(6.-7.*s)*log(s)
204 +24.*(1.-s)*(1.-s)*log(s)*log(1.-s)+12.*(-13.+16.*s-3.*s*s)*(log(1.-sqrt(s))-log(1.-s))
205 +39.-2.*PI2+252.*s-26.*PI2*s+21.*s*s+8.*PI2*s*s-180.*sqrt(s)-132.*s*sqrt(s));
206}
207
208double BXsllDecay::f_9(double s) {
209 return -1./6./(s-1.)/(s-1.)*(48.*s*(-5.+2.*s)*Li2(sqrt(s))+24.*(-1.+7.*s-3.*s*s)*Li2(s)+6.*s*(-6.+7.*s)*log(s)
210 -24.*(1.-s)*(1.-s)*log(s)*log(1.-s)+24.*(5.-7.*s+2.*s*s)*(log(1.-sqrt(s))-log(1.-s))
211 -21.-156.*s+20.*PI2*s+9.*s*s-8.*PI2*s*s+120.*sqrt(s)+48.*s*sqrt(s));
212}
213
215 if(t>4.) return -2.*I*PI*log((sqrt(t)+sqrt(t-4.))/2.)-PI2/2.+2.*pow(log((sqrt(t)+sqrt(t-4.))/2.),2.);
216 else return 2.*PI*atan(sqrt((4.-t)/t))-PI2/2.-2.*pow(atan(sqrt((4.-t)/t)),2.);
217}
218
220 if(t>4.) return -I*PI*sqrt((t-4.)/t)-2.+2.*sqrt((t-4.)/t)*log((sqrt(t)+sqrt(t-4.))/2.);
221 else return PI*sqrt((4.-t)/t)-2.-2.*sqrt((4.-t)/t)*atan(sqrt((4.-t)/t));
222}
223
224complex_t BXsllDecay::Delta_i_23(double s, double z, double w) {
225 return -2.+4./(w-s)*(z*Gm1(s/z)-z*Gm1(w/z)-s/2.*G0(s/z)+s/2.*G0(w/z));
226}
227
228complex_t BXsllDecay::Delta_i_27(double s, double z, double w) {
229 return 2.*(G0(s/z) - G0(w/z));
230}
231
232double BXsllDecay::tau_22(double s, double w, complex_t Delta_23, complex_t Delta_27) {
233 return 8./27.*(w-s)*(1.-w)*(1.-w)/s/w/w/w*((3.*w*w+2.*s*s*(2.+w)-s*w*(5.-2.*w))*pow(abs(Delta_23),2.)
234 +(2.*s*s*(2.+w)+s*w*(1.+2.*w))*pow(abs(Delta_27),2.)
235 +4.*s*(w*(1.-w)-s*(2.+w))*real(Delta_23*conj(Delta_27)));
236}
237
238complex_t BXsllDecay::tau_27(double s, double w, complex_t Delta_23, complex_t Delta_27) {
239 return 8./3./s/w*(((1.-w)*(4.*s*s-s*w+w*w)+s*w*(4.+s-w)*log(w))*Delta_23
240 -(4.*s*s*(1.-w)+s*w*(4.+s-w)*log(w))*Delta_27);
241}
242
243complex_t BXsllDecay::tau_28(double s, double w, complex_t Delta_23, complex_t Delta_27) {
244 return 8./9./s/w/(w-s)*((pow(w-s,2.)*(2.*s-w)*(1.-w))*Delta_23
245 -(2.*s*pow(w-s,2.)*(1.-w))*Delta_27
246 +s*w*((1.+2.*s-2.*w)*Delta_23-2.*(1.+s-w)*Delta_27)*log(s/((1.+s-w)*(w*w+s*(1.-w)))));
247}
248
249complex_t BXsllDecay::tau_29(double s, double w, complex_t Delta_23, complex_t Delta_27) {
250 return 4./3./w*((2.*s*(1.-w)*(s+w)+4.*s*w*log(w))*Delta_23
251 -(2.*s*(1.-w)*(s+w)+w*(3.*s+w)*log(w))*Delta_27);
252}
253
254double BXsllDecay::tau_77(double s) {
255 return -2./9./(2.+s)*(2.*(1.-s)*(1.-s)*log(1.-s)+6.*s*(2.-2.*s-s*s)/(1.-s)/(1.-s)*log(s)+(11.-7.*s-10.*s*s)/(1.-s));
256}
257
258double BXsllDecay::tau_78(double s) {
259 return 8./9./s*(25.-2.*PI2-27.*s+3.*s*s-s*s*s+12.*(s+s*s)*log(s)+6.*pow((PI/2.-atan((2.-4.*s+s*s)/(2.-s)*sqrt(s)*sqrt(4.-s))),2.)
260 -24.*real(CLi2((s-I*sqrt(s)*sqrt(4.-s))/2.))-12.*((1.-s)*sqrt(s)*sqrt(4.-s)-atan((sqrt(s)*sqrt(4.-s))/(2.-s)))
261 *(atan(sqrt((4.-s)/s))-atan((sqrt(s)*sqrt(4.-s))/(2.-s))));
262}
263
264double BXsllDecay::tau_88(double s) {
265 return 4./27./s*(-8.*PI2+(1.-s)*(77.-s-4.*s*s)-24.*Li2(1.-s)+3.*(10.-4.*s-9.*s*s+8.*log((sqrt(s))/(1.-s)))*log(s)
266 +48.*real(CLi2((3.-s)/2.+I*(1.-s)*sqrt(4.-s)/2./sqrt(s)))-6.*((20.*s+10.*s*s-3.*s*s*s)/sqrt(s)/sqrt(4.-s)-8.*PI+8.*atan(sqrt((4.-s)/s)))
267 *(atan(sqrt((4.-s)/s))-atan((sqrt(s)*sqrt(4.-s))/(2.-s))));
268}
269
270double BXsllDecay::tau_89(double s) {
271 return 2./3.*(s*(4.-s)-3.-4.*log(s)*(1.-s-s*s)
272 -8.*real(CLi2(s/2.+I*sqrt(s)*sqrt(4.-s)/2.)-CLi2((-2.+s*(4.-s))/2.+I*((2.-s)*sqrt(s)*sqrt(4.-s))/2.))
273 +4.*(s*s*sqrt((4.-s)/s)+2.*atan(sqrt(s)*sqrt(4.-s)/(2.-s)))*(atan(sqrt((4.-s)/s))-atan((sqrt(s)*sqrt(4.-s))/(2.-s))));
274}
275
276double BXsllDecay::tau_99(double s) {
277 return -4./9./(1.+2.*s)*(2.*(1.-s)*(1.-s)*log(1.-s)+3.*s*(1.+s)*(1.-2.*s)/(1.-s)/(1.-s)*log(s)+3.*(1.-3.*s*s)/(1.-s));
278}
279
280double BXsllDecay::tau_79(double s) {
281 return -4.*(1.-s)*(1.-s)/9./s*log(1.-s)-4.*s*(3.-2.*s)*log(s)/9./(1.-s)/(1.-s)-2./9.*(5.-3.*s)/(1.-s);
282}
283
284complex_t BXsllDecay::tau_210(double s, double z) {
285 auto f = [&] (double w) {
286 return -s/(s-w)/(1.-s)/(1.-s)*((4.*(1.-s)*(1.+w)-2.*fabs(s-w*w)*(w*(3.+w)-s*(1.-w))/w/w
287 +(2.+5.*w+2.*w*w+s*(3.+4.*w))*log((s+w*w+fabs(s-w*w))/2./w)-(s-w)/s/sqrt((1.+w)*(1.+w)-4.*s)
288 *(w*(2.-w)-s*(6.-5.*w))*(log(1.+w-s*(3.-w)+(1.-s)*sqrt((1.+w)*(1.+w)-4.*s))
289 -log(s*(1.-3.*w)+w*w*(1.+w)+fabs(s-w*w)*sqrt((1.+w)*(1.+w)-4.*s))))*Delta_i_23(s, z, w)
290 -(2.*(1.-s)*(1.+2.*w)-2.*fabs(s-w*w)*(w*(2.+w)-s*(1.-w))/w/w
291 +2.*(s*(1.+2.*w)+w*(2.+w))*log((s+w*w+fabs(s-w*w))/2./w)
292 +4.*(1.-w)*(s-w)/sqrt((1.+w)*(1.+w)-4.*s)*(log(1.+w-s*(3.-w)+(1.-s)*sqrt((1.+w)*(1.+w)-4.*s))
293 -log(s*(1.-3.*w)+w*w*(1.+w)+fabs(s-w*w)*sqrt((1.+w)*(1.+w)-4.*s))))*Delta_i_27(s, z, w));
294 };
295
296 return c_integrate(f, s, 1, 1e-2);
297}
298
299double BXsllDecay::tau_710(double s) {
300 return -5./2.+1./3./(1.-3.*s)-1./3.*s*(6.-7.*s)*log(s)/(1.-s)/(1.-s)-1./9.*(3.-7.*s+4.*s*s)*log(1.-s)/s+f_7(s)/3.;
301}
302
303double BXsllDecay::tau_810(double s){
304 return 1./6./(1.-s)/(1.-s)*(3.*((1.-sqrt(s))*(1.-sqrt(s))*(23.-6.*sqrt(s)-s)+4.*(1.-s)*(7.+s)*log(1.+sqrt(s))
305 +2.*s*(1.+s-log(s))*log(s))+2.*(-3.*PI2*(1.+2.*s)+6.*(3.-s)*s*log(2.-sqrt(s))
306 -36.*(1.+2.*s)*CLi2(-sqrt(s))-6.*sqrt(s/(4.-s))*(2.*(-3.+s)*s*atan((2.+sqrt(s))/sqrt(4.-s))
307 +2.*PI*log(2.-sqrt(s))-atan(sqrt((4.-s)/s))*((-3.+s)*s+4.*log(2.-sqrt(s)))
308 -atan(sqrt(s*(4.-s))/(2.-s))*((-3.+s)*s-log(s))+4.*real(I*CLi2((-2.+I*sqrt(4.-s)+sqrt(s))*sqrt(s)/(I*sqrt(4.-s)-sqrt(s))))
309 -2.*real(I*CLi2(I/2.*sqrt(4.-s)*(1.-s)*sqrt(s)+(3.-s)*s/2.)))));
310}
311
312double BXsllDecay::tau_910(double s) {
313 return -5./2.+1./3./(1.-s)-1./3.*s*(6.-7.*s)*log(s)/(1.-s)/(1.-s)-2./9.*(3.-5.*s+2.*s*s)*log(1.-s)/s+f_9(s)/3.;
314}
315
316double BXsllDecay::sigma(double s) {
317 return -4./3.*Li2(s)-2./3.*log(s)*log(1.-s)-2./9.*PI2-log(1.-s)-2./9.*(1.-s)*log(1.-s);
318}
319
320double BXsllDecay::sigma_9(double s) {
321 return sigma(s)+1.5;
322}
323
324double BXsllDecay::sigma_7(double s, double L_mu) {
325 return sigma(s)+1./6.-8./3.*L_mu;
326}
327
329 if(r < 1) {
330 return 3./2./r*(1./sqrt(r*(1.-r))*atan(sqrt(r/(1.-r)))-1.);
331 }
332 return 3./2./r*(1./2./sqrt(r*(r-1.))*(log((1.-sqrt(1.-1./r))/(1.+sqrt(1.-1./r)))+I*PI)-1.);
333}
334
336 if (abs(s) < 0.4)
337 return {
338 23.787-120.948*s+365.373*s*s-584.206*s*s*s,
339 1.653+6.009*s-17.080*s*s+115.880*s*s*s
340 };
341
342 double d = 1 - s;
343 return {
344 -148.061*d*d+492.539*d*d*d-1163.847*pow(d,4.)+1189.528*pow(d,5.),
345 -261.287*d*d+1170.856*d*d*d-2546.948*pow(d,4.)+2540.023*pow(d,5.)
346 };
347}
348
349double BXsllDecay::Sigma_2(double s) {
350 if (abs(s) < 0.4)
351 return 11.488-36.987*s+255.330*s*s-812.388*s*s*s+1011.791*s*s*s*s;
352
353 double d = 1 - s;
354 return -221.904*d*d+900.822*d*d*d-2031.620*pow(d,4.)+1984.303*pow(d,5.);
355}
356
358 if (abs(s) < 0.4)
359 return {
360 109.311-846.039*s+2890.115*s*s-4179.072*s*s*s,
361 4.606+17.650*s-53.244*s*s+348.069*s*s*s
362 };
363
364 double d = 1 - s;
365 return {
366 -298.730*d*d+828.0675*d*d*d-2217.6355*pow(d,4.)+2241.792*pow(d,5.),
367 -528.759*d*d+2095.723*d*d*d-4681.843*pow(d,4.)+5036.677*pow(d,5.)
368 };
369}
370
371complex_t BXsllDecay::Sigma_7(double s, double z) {
372 if (abs(s) < 0.4) {
373 double a = pow(4 * z, 2);
374 return {
375 -0.259023-28.424*s+205.533*s*s-603.219*s*s*s+722.031*s*s*s*s,
376 (-12.20658-215.8208*(s-a)+412.1207*(s-a)*(s-a))*(s-a)*(s-a)*(s>a)
377 };
378 }
379
380 double d = 1 - s;
381 return {
382 77.0256*d*d-264.705*d*d*d+595.814*pow(d,4.)-610.1637*pow(d,5.),
383 135.858*d*d-618.990*d*d*d+1325.040*pow(d,4.)-1277.170*pow(d,5.)
384 };
385}
386
387double BXsllDecay::omega_22(double s, double L_l) {
388 return L_l * (Sigma_2(s)/8./(1.-s)/(1.-s)/(1.+2.*s)+real(Sigma_1(s))/9./(1.-s)/(1.-s)/(1.+2.*s)*cache.L_b_5GeV)
389 +64./81.*omega_1010(s, L_l) * cache.L_b_5GeV * cache.L_b_5GeV;
390}
391
392complex_t BXsllDecay::omega_27(double s, double L_l) {
393 return L_l*(Sigma_3(s)/96./(1.-s)/(1.-s))+8./9.*omega_79(s, L_l)*cache.L_b_5GeV;
394}
395
396complex_t BXsllDecay::omega_29(double s, double L_l) {
397 return L_l * (Sigma_1(s)/8./(1.-s)/(1.-s)/(1.+2.*s))+16./9.*omega_1010(s,L_l)*cache.L_b_5GeV;
398}
399
400complex_t BXsllDecay::omega_210(double s, double L_l, double z) {
401 return L_l*(-Sigma_7(s, z)/24./s/(1.-s)/(1.-s))+8./9.*omega_910(s,L_l)*cache.L_b_5GeV;
402}
403
404double BXsllDecay::omega_77(double s, double L_l) {
405 return L_l*(s/2./(1.-s)/(2.+s)+log(1.-s)-s*(-3.+2.*s*s)/2./(1.-s)/(1.-s)/(2.+s)*log(s));
406}
407
408double BXsllDecay::omega_79(double s, double L_l) {
409 return L_l*(-0.5/(1.-s)+log(1.-s)+(-1.+2.*s-2.*s*s)/2./(1.-s)/(1.-s)*log(s));
410}
411
412double BXsllDecay::omega_710(double s, double L_l) {
413 return L_l*((7.-16.*sqrt(s)+9.*s)/4./(1.-s)+log(1.-sqrt(s))+(1.+3.*s)/(1.-s)*log((1.+sqrt(s))/2.)-s*log(s)/(1.-s));
414}
415
416double BXsllDecay::omega_99(double s, double L_l) {
417 return L_l*(-(1.+4.*s-8.*s*s)/6./(1.-s)/(1.+2.*s)+log(1.-s)-(1.-6.*s*s+4.*s*s*s)*log(s)/2./(1.-s)/(1.-s)/(1.+2.*s))
418 -Li2(s)/9.+4./27.*PI2-(37.-3.*s-6.*s*s)/72./(1.-s)/(1.+2.*s)-((41.+76.*s)*log(1.-s))/36./(1.+2.*s)
419 +((6.-10.*s-17.*s*s+14.*s*s*s)/18./(1.-s)/(1.-s)/(1.+2.*s)+17.*log(1.-s)/18.)*log(s)-(1.-6.*s*s+4.*s*s*s)*log(s)*log(s)/2./(1.-s)/(1.-s)/(1.+2.*s);
420}
421
422double BXsllDecay::omega_910(double s, double L_l) {
423 return L_l*(-(5.-16.*sqrt(s)+11.*s)/4./(1.-s)+log(1.-sqrt(s))+(1.-5.*s)/(1.-s)*log((1.+sqrt(s))/2.)-(1.-3.*s)*log(s)/(1.-s));
424}
425
426double BXsllDecay::omega_1010(double s, double L_l) {
427 return L_l*(-(1.+4.*s-8.*s*s)/6./(1.-s)/(1.+2.*s)+log(1.-s)-(1.-6.*s*s+4.*s*s*s)*log(s)/2./(1.-s)/(1.-s)/(1.+2.*s));
428}
429
430complex_t BXsllDecay::g(double z, double s) {
431 double z2=z*z;
432
433 if(s==0.)
434 return -4./9.*log(z2)+8./27.-4./9.;
435
436 if(z==0.)
437 return 8./27.-4./9.*(log(s)-I*PI);
438
439 if(4.*z2<s) {
440 return -4./9.*log(z2)+8./27.+16./9.*z2/s -2./9.*sqrt(1.-4.*z2/s)*(2.+4.*z2/s)*(log((1.+sqrt(1.-4.*z2/s))/(1.-sqrt(1.-4.*z2/s)))-I*PI);
441 } else if(4.*z2>s) {
442 return -4./9.*log(z2)+8./27.+16./9.*z2/s -4./9.*sqrt(4.*z2/s-1.)*(2.+4.*z2/s)*atan(1./sqrt(4.*z2/s-1.));
443 }
444
445 else
446 return -4./9.*log(z2)+8./27.+16./9.*z2/s;
447}
448
449double BXsllDecay::breit_wigner(double s, double m_V, double br, double gamma_tot, double gamma_had) {
450 double m_V_hat = m_V / cache.m_b_1S;
451 double gamma_tot_hat = gamma_tot / cache.m_b_1S;
452 double gamma_had_hat = gamma_had / cache.m_b_1S;
453 // LOG_INFO("m_V_hat", m_V_hat);
454 // LOG_INFO("gamma_tot_hat", gamma_tot_hat);
455 // LOG_INFO("gamma_had_hat", gamma_had_hat);
456 // LOG_INFO("br * gamma_tot_hat * gamma_had_hat",br * gamma_tot_hat * gamma_had_hat);
457 // LOG_INFO("denom", (std::pow(s - std::pow(m_V_hat, 2), 2) + std::pow(m_V_hat * gamma_tot_hat, 2)));
458 // LOG_INFO("tout", br * gamma_tot_hat * gamma_had_hat / (std::pow(s - std::pow(m_V_hat, 2), 2) + std::pow(m_V_hat * gamma_tot_hat, 2)));
459 return br * gamma_tot_hat * gamma_had_hat / (std::pow(s - std::pow(m_V_hat, 2), 2) + std::pow(m_V_hat * gamma_tot_hat, 2));
460}
461
462double BXsllDecay::R_cc_cont(double s) {
463 return s > 0.6 ? (s > 0.69 ? 1.02 : 11.33 * s - 6.8) : 0;
464}
465
466double BXsllDecay::R_cc(double s) {
467 double R_cc_res = 0;
468
469 for (size_t k = 0; k < cache.cc_res_mass.size(); k++) {
470 R_cc_res += breit_wigner(s, cache.cc_res_mass[k], cache.cc_res_br[k], cache.cc_res_width_tot[k], cache.cc_res_width_had[k]);
471 }
472
473 return 9. * s / std::pow(cache.alpha_em, 2) * R_cc_res + R_cc_cont(s);
474}
475
476double BXsllDecay::PV_breit_wigner(double s, double m_V, double br, double gamma_tot, double gamma_had) {
477 double s_c = 4 * pow(cache.m_D_hat, 2);
478 double m_V_hat = m_V / cache.m_b_1S;
479 double m_V_hat2 = pow(m_V_hat, 2);
480 double gamma_tot_hat = gamma_tot / cache.m_b_1S;
481 double den = (pow(s - pow(m_V_hat, 2), 2) + pow(m_V_hat * gamma_tot_hat, 2));
482 double B = breit_wigner(s, m_V, br, gamma_tot, gamma_had);
483 return 9. / pow(cache.alpha_em, 2) * B * (0.5 * log(den / pow(s_c - s, 2)) + (s - m_V_hat2) / gamma_tot_hat * m_V_hat * (atan((s_c - m_V_hat2) / gamma_tot_hat * m_V_hat) - PI / 2.));
484}
485
486double BXsllDecay::PV_R_cc_cont(double s) {
487 double s_c = 4 * pow(cache.m_D_hat, 2);
488 return 1 / s * ((11.33 * s - 6.8) * log(abs((0.69 - s) / (s_c - s))) - 1.02 * log(abs(0.69 - s)) - 6.8 * log(s_c) - 2.90171798847631);
489}
490
491double BXsllDecay::PV_R_cc(double s) {
492 double PV_res = 0;
493
494 for (size_t k = 0; k < cache.cc_res_mass.size(); k++) {
495 PV_res += PV_breit_wigner(s, cache.cc_res_mass[k], cache.cc_res_br[k], cache.cc_res_width_tot[k], cache.cc_res_width_had[k]);
496 }
497
498 return PV_res + PV_R_cc_cont(s);
499}
500
501complex_t BXsllDecay::g_ld(double z, double s) {
502 // printf("g(z, 0) = %.4e + %.4e i\n", real(g(z, 0)), imag(g(z, 0)));
503 // printf("PV_R_cc(s = %.3f) = %.4e + %.4e i\n", s, real(PV_R_cc(s)), imag(PV_R_cc(s)));
504 // printf("R_cc(s = %.3f) = %.4e + %.4e i\n", s, real(R_cc(s)), imag(R_cc(s)));
505
506 return g(z, 0) + (s * PV_R_cc(s) + I * PI * R_cc(s)) / 3.;
507}
508
509complex_t BXsllDecay::C9_eff(double s, QCDOrder order, bool prime) {
510 s = std::clamp(s, 1e-6, 1. - 1e-6);
511 auto C = order == QCDOrder::LO ? cache.C_LO : cache.C;
512 auto C_ids = WCoefMapper::get_group(prime ? WGroup::BPrime : WGroup::B);
513 complex_t g_0 = g(0, s);
514 complex_t g_1 = g(1, s);
515 complex_t g_mc = g_ld(cache.m_c_hat, s);
516
517 // printf("g_0(s = %.3f) = %.4e + %.4e i\n", s, real(g_0), imag(g_0));
518 // printf("g_1(s = %.3f) = %.4e + %.4e i\n", s, real(g_1), imag(g_1));
519 // printf("g_mc(s = %.3f) = %.4e + %.4e i\n", s, real(g_mc), imag(g_mc));
520
521 return C[C_ids[8]]
522 -(-32./27.*C[C_ids[0]]-8./9.*C[C_ids[1]]-16./9.*C[C_ids[2]]+32./27.*C[C_ids[3]]-112./9.*C[C_ids[4]]+512./27.*C[C_ids[5]]) * cache.L_b
523 +4./3.*C[C_ids[2]]+64./9.*C[C_ids[4]]+64./27.*C[C_ids[5]]
524 +g_mc*(4./3.*C[C_ids[0]]+C[C_ids[1]]+6.*C[C_ids[2]]+60.*C[C_ids[4]])
525 +g_1*(-7./2.*C[C_ids[2]]-2./3.*C[C_ids[3]]-38.*C[C_ids[4]]-32./3.*C[C_ids[5]])
526 +g_0*(-1./2.*C[C_ids[2]]-2./3.*C[C_ids[3]]-8.*C[C_ids[4]]-32./3.*C[C_ids[5]]);
527}
528
530 return lerp(s, cache.F_17_lookup);
531}
532
534 return lerp(s, cache.F_27_lookup);
535}
536
538 return lerp(s, cache.F_19_lookup);
539}
540
542 return lerp(s, cache.F_29_lookup);
543}
544
545complex_t BXsllDecay::C7_new(double s, bool prime) {
546 auto C_0 = cache.C_LO;
547 auto C7_eff = cache.C[prime ? WCoef::CP7 : WCoef::C7];
548 // LOG_INFO("s_hat =", s);
549 // LOG_INFO("C7eff =", C7_eff);
550 // LOG_INFO("sigma_7 =", sigma_7(s, cache.L_b));
551 // LOG_INFO("F_17 =", F_17(s));
552 // LOG_INFO("F_27 =", F_27(s));
553 // LOG_INFO("F_87 =", f_87(s, cache.L_b));
554 return (1.+cache.alpha_s_mu_b/PI*sigma_7(s, cache.L_b))*C7_eff
555 -cache.alpha_s_mu_b/4./PI*(C_0.at(prime ? WCoef::CP1 : WCoef::C1)*F_17(s)+C_0.at(prime ? WCoef::CP2 :WCoef::C2)*F_27(s)+C_0.at(prime ? WCoef::CP8 :WCoef::C8)*BV::f_87(s, cache.L_b));
556}
557
558complex_t BXsllDecay::C9_new(double s, bool prime) {
559 if (abs(s - 1) < 1e-6) s = 1;
560 auto C_0 = cache.C_LO;
561 // printf("C9_eff = %.4e + %.4e i\n", real(C9_eff(s, this->w_config.order, prime)), imag(C9_eff(s, this->w_config.order, prime)));
562 // printf("F_19 = %.4e + %.4e i\n", real(F_19(s)), imag(F_19(s)));
563 // printf("F_29 = %.4e + %.4e i\n", real(F_29(s)), imag(F_29(s)));
564 // printf("F_89 = %.4e + %.4e i\n", real(BV::f_89(s)), imag(BV::f_89(s)));
565
566 return (1.+cache.alpha_s_mu_b/PI*sigma_9(s))*C9_eff(s, this->w_config.order, prime)
567 -cache.alpha_s_mu_b/4./PI*(C_0.at(prime ? WCoef::CP1 :WCoef::C1)*F_19(s)+C_0.at(prime ? WCoef::CP2 : WCoef::C2)*F_29(s)+C_0.at(prime ? WCoef::CP8 :WCoef::C8)*BV::f_89(s));
568}
569
570complex_t BXsllDecay::C10_new(double s, bool prime) {
571 s = std::clamp(s, 1e-6, 1. - 1e-6);
572 return (1.+cache.alpha_s_mu_b/PI*sigma_9(s))*(cache.C[prime ? WCoef::CP10 : WCoef::C10]);
573}
574
575double BXsllDecay::W_7(double s) {
576 return pow(abs(C7_new(s, false)), 2) + pow(abs(C7_new(s, true)), 2);
577}
578
579double BXsllDecay::W_9(double s) {
580 return pow(abs(C9_new(s, false)), 2) + pow(abs(C9_new(s, true)), 2);
581}
582
583double BXsllDecay::W_10(double s) {
584 return pow(abs(C10_new(s, false)), 2) + pow(abs(C10_new(s, true)), 2);
585}
586
588 return cache.C[WCoef::C2] * conj(C7_new(s, false)) + cache.C[WCoef::CP2] * conj(C7_new(s, true));
589}
590
592 return cache.C[WCoef::C2] * conj(C9_new(s, false)) + cache.C[WCoef::CP2] * conj(C9_new(s, true));
593}
594
596 return cache.C[WCoef::C2] * conj(C10_new(s, false)) + cache.C[WCoef::CP2] * conj(C10_new(s, true));
597}
598
599double BXsllDecay::W_79(double s) {
600 return real(C7_new(s, false) * conj(C9_new(s, false))) + real(C7_new(s, true) * conj(C9_new(s, true)));
601}
602
603double BXsllDecay::W_710(double s) {
604 return real(C7_new(s, false) * conj(C10_new(s, false))) + real(C7_new(s, true) * conj(C10_new(s, true)));
605}
606
607double BXsllDecay::W_910(double s) {
608 return real(C9_new(s, false) * conj(C10_new(s, false))) + real(C9_new(s, true) * conj(C10_new(s, true)));
609}
610
611double BXsllDecay::dB0_ds(double s, double ml_hat) {
612 double W_Q1 = pow(abs(cache.C[WCoefMapper::cq1_for_lepton_index(static_cast<int>(cfg.gen))]), 2) + pow(abs(cache.C[WCoefMapper::cpq1_for_lepton_index(static_cast<int>(cfg.gen))]), 2);
613 double W_Q2 = pow(abs(cache.C[WCoefMapper::cq2_for_lepton_index(static_cast<int>(cfg.gen))]), 2) + pow(abs(cache.C[WCoefMapper::cpq2_for_lepton_index(static_cast<int>(cfg.gen))]), 2);
614 double W_10Q2 = real(cache.C[WCoefMapper::cq2_for_lepton_index(static_cast<int>(cfg.gen))] * conj(C10_new(s, false))) + real(cache.C[WCoefMapper::cpq2_for_lepton_index(static_cast<int>(cfg.gen))] * conj(C10_new(s, true)));
615
616 double H7 = 4 * (1. + 2. * ml_hat * ml_hat / s) * (1. + 2. / s) * (1. + cache.alpha_s_mu_b / PI * tau_77(s));
617 double H9 = (1. + 2. * ml_hat * ml_hat / s) * (1. + 2. * s) * (1. + cache.alpha_s_mu_b / PI * tau_99(s));
618 double H10 = ((1.+2.*s)+2.*ml_hat*ml_hat/s*(1.-4.*s))*(1.+cache.alpha_s_mu_b/PI*tau_99(s));
619 double H79 = 12 * (1.+2.*ml_hat*ml_hat/s)*(1.+cache.alpha_s_mu_b/PI*tau_79(s));
620
621 // LOG_INFO("W_7 =", W_7(s));
622 // LOG_INFO("W_9 =", W_9(s));
623 // LOG_INFO("W_10 =", W_10(s));
624 // LOG_INFO("W_79 =", W_79(s));
625
626 return pow(1 - s, 2) * sqrt(1 - 4 * pow(ml_hat, 2) / s) * (
627 H7 * W_7(s) +
628 H9 * W_9(s) +
629 H10 * W_10(s) +
630 H79 * W_79(s) +
631 1.5 * (s - 4 * ml_hat * ml_hat) * W_Q1 +
632 1.5 * s * W_Q2 +
633 6 * ml_hat * W_10Q2
634 );
635}
636
637double BXsllDecay::A_FB_0(double s, double ml_hat) {
638 double W_7Q1 = real(C7_new(s, false) * cache.C[WCoefMapper::cq1_for_lepton_index(static_cast<int>(cfg.gen))]) + real(C7_new(s, true) * cache.C[WCoefMapper::cpq1_for_lepton_index(static_cast<int>(cfg.gen))]);
639 double W_9Q1 = real(C9_new(s, false) * cache.C[WCoefMapper::cq1_for_lepton_index(static_cast<int>(cfg.gen))]) + real(C9_new(s, true) * cache.C[WCoefMapper::cpq1_for_lepton_index(static_cast<int>(cfg.gen))]);
640
641 return pow(1 - s, 2) * sqrt(1 - 4 * pow(ml_hat, 2) / s) * (
642 2 * (1 + cache.alpha_s_mu_b * tau_710(s) / PI) * W_710(s) +
643 s * (1 + cache.alpha_s_mu_b * tau_910(s) / PI) * W_910(s) +
644 ml_hat * (2 * W_7Q1 + W_9Q1)
645 );
646}
647
648double BXsllDecay::delta_A_mb2(double s) {
649 return s * (9 + 14 * s - 15 * pow(s, 2)) * W_910(s) + 2 * (7 + 10 * s - 9 * s * s) * W_710(s);
650}
651
652double BXsllDecay::delta_mb2(double s) {
653 return -4 * (6 + 3 * s - 5 * pow(s, 3)) * W_7(s) / s + (1 - 15 * s * s + 10 * pow(s, 3)) * (W_9(s) + W_10(s)) - 4 * (5 + 6 * s - 7 * s * s) * W_79(s);
654}
655
656double BXsllDecay::delta_mb3(double s) {
657 if (s > 0.4)
658 return 0;
659
660 return (5.*pow(s,4.)+19.*pow(s,3.)+9.*s*s-7.*s+22.)/6./(1.-s)*4.*W_7(s)/s+
661 (10.*pow(s,4.)+23.*pow(s,3.)-9.*s*s+13.*s+11.)/6./(1.-s)*(W_9(s) + W_10(s))+
662 4.*(-3.*pow(s,3.)+17.*s*s-s+3.)/2./(1.-s)* W_79(s);
663}
664
665double BXsllDecay::delta_mc2(double s) {
666 complex_t f = F(s / (4. * cache.z));
667 return pow(1 - s, 2) * real((1 + 6 * s - s * s) * f * W_27(s) / s + (2 + s) * f * W_29(s));
668}
669
670double BXsllDecay::delta_A_mc2(double s) {
671 return pow(1 - s, 2) * real((1 + 3 * s) * F(s / (4. * cache.z)) * W_210(s));
672}
673
674double BXsllDecay::delta_bremA(double s) {
675 complex_t C7_0 = cache.C_LO[WCoef::C7];
676 complex_t C8_0 = cache.C_LO[WCoef::C8];
677 complex_t CP7_0 = cache.C_LO[WCoef::CP7];
678 complex_t CP8_0 = cache.C_LO[WCoef::CP8];
679 complex_t C9_0 = C9_eff(s, QCDOrder::LO, false);
680 complex_t CP9_0 = C9_eff(s, QCDOrder::LO, true);
681
682 // printf("C9_0_eff (s = %.4f) = %.5e + %.5e i\n", s, real(C9_0), imag(C9_0));
683 // printf("CP9_0_eff (s = %.4f) = %.5e + %.5e i\n", s, real(CP9_0), imag(CP9_0));
684
685 complex_t c_78 = (*iobs_qcdp).get_constants()->C_F * (C7_0 * conj(C8_0) + CP7_0 * conj(CP8_0));
686 complex_t c_88 = (*iobs_qcdp).get_constants()->C_F * (C8_0 * conj(C8_0) + CP8_0 * conj(CP8_0));
687 complex_t c_89 = (*iobs_qcdp).get_constants()->C_F * (C8_0 * conj(C9_0) + CP8_0 * conj(CP9_0));
688
689 // printf("c_78 (s = %.4f) = %.4e + %.4e i\n", s, real(c_78), imag(c_78));
690 // printf("c_88 (s = %.4f) = %.4e + %.4e i\n", s, real(c_88), imag(c_88));
691 // printf("c_78 (s = %.4f) = %.4e + %.4e i\n", s, real(c_89), imag(c_89));
692
693 return 2. * real(c_78 * tau_78(s) + c_89 * tau_89(s) + 0.5 * c_88 * tau_88(s));
694}
695
697 constexpr double eps = 1e-6;
698 s = std::clamp(s, eps, 1.0 - eps);
699
700 auto C_0 = cache.C_LO;
701 auto CP_0 = cache.C_LO;
702 complex_t C9_0 = C9_eff(s, QCDOrder::LO, false);
703
704 double C_f = (*iobs_qcdp).get_constants()->C_F;
705 double C_tau_1 = C_f / (4. * pow((*iobs_qcdp).get_constants()->Nc, 2));
706 double C_tau_2 = -C_f / (2. * (*iobs_qcdp).get_constants()->Nc);
707
708 complex_t c_11 = C_tau_1 * (C_0[WCoef::C1] * conj(C_0[WCoef::C1]) + CP_0[WCoef::CP1] * conj(CP_0[WCoef::CP1]));
709 complex_t c_12 = 2 * C_tau_2 * real((C_0[WCoef::C1]*conj(C_0[WCoef::C2])) + CP_0[WCoef::CP1] * conj(CP_0[WCoef::CP2]));
710 complex_t c_22 = C_f*(C_0[WCoef::C2]*conj(C_0[WCoef::C2]) + CP_0[WCoef::CP2] * conj(CP_0[WCoef::CP2]));
711 complex_t c_17 = C_tau_2*(C_0[WCoef::C1]*conj(C_0[WCoef::C7]) + CP_0[WCoef::CP1] * conj(CP_0[WCoef::CP7]));
712 complex_t c_27 = C_f*(C_0[WCoef::C2]*conj(C_0[WCoef::C7]) + CP_0[WCoef::CP2] * conj(CP_0[WCoef::CP7]));
713 complex_t c_18 = C_tau_2*(C_0[WCoef::C1]*conj(C_0[WCoef::C8]) + CP_0[WCoef::CP1] * conj(CP_0[WCoef::CP8]));
714 complex_t c_28 = C_f*(C_0[WCoef::C2]*conj(C_0[WCoef::C8]) + CP_0[WCoef::CP2] * conj(CP_0[WCoef::CP8]));
715 complex_t c_19 = C_tau_2*(C_0[WCoef::C1]*conj(C9_0) + CP_0[WCoef::CP1] * conj(CP_0[WCoef::CP9]));
716 complex_t c_29 = C_f*(C_0[WCoef::C2]*conj(C9_0) + CP_0[WCoef::CP2] * conj(CP_0[WCoef::CP9]));
717
718 complex_t w_22 = c_11 + c_12 + c_22;
719 complex_t w_27 = c_17 + c_27;
720 complex_t w_28 = c_18 + c_28;
721 complex_t w_29 = c_19 + c_29;
722
723 // printf("w_22 (s = %.4f) = %.4e + %.4e i\n", s, real(w_22), imag(w_22));
724 // printf("w_27 (s = %.4f) = %.4e + %.4e i\n", s, real(w_27), imag(w_27));
725 // printf("w_28 (s = %.4f) = %.4e + %.4e i\n", s, real(w_28), imag(w_28));
726 // printf("w_29 (s = %.4f) = %.4e + %.4e i\n", s, real(w_29), imag(w_29));
727
728 auto f = [&] (double w) -> double {
729 complex_t D23 = Delta_i_23(s, cache.z, w);
730 complex_t D27 = Delta_i_27(s, cache.z, w);
731 return real(w_22 * tau_22(s, w, D23, D27) + 2.0 * (w_27 * tau_27(s, w, D23, D27) + w_28 * tau_28(s, w, D23, D27) + w_29 * tau_29(s, w, D23, D27)));
732 };
733
734
735 return integrate(f, s, 1.0 - eps, 1e-2);
736 // return integrate(f, s, 1, 1e-2);
737}
738
739double BXsllDecay::delta_bremB(double s) {
740 return lerp(s, cache.delta_brems_lookup, 0, 1);
741}
742
743double BXsllDecay::delta_A_brem(double s) {
744 auto C_0 = cache.C_LO;
745 auto CP_0 = cache.C_LO;
746 complex_t C10 = cache.C[WCoef::C10];
747 complex_t CP10 = cache.C[WCoef::CP10];
748
749 complex_t W_210 = (C_0[WCoef::C2] - C_0[WCoef::C1] / 6.) * C10 + (CP_0[WCoef::CP2] - CP_0[WCoef::CP1] / 6.) * CP10;
750 complex_t W_810 = C_0[WCoef::C8] * C10 + CP_0[WCoef::CP8] * CP10;
751
752 return pow(1 - s, 2) * real(W_810 * tau_810(s) + W_210 * tau_210(s, cache.z));
753}
754
755double BXsllDecay::delta_em(double s, double L_l) {
756 auto C = cache.C;
757 auto Cp = cache.C;
758
759 double C_F = (*iobs_qcdp).get_constants()->C_F;
760 complex_t W_2 = pow(abs(C[WCoef::C2] + C_F * C[WCoef::C1]), 2) + pow(abs(Cp[WCoef::CP2] + C_F * Cp[WCoef::CP1]), 2);
761 complex_t W_7 = pow(abs(C[WCoef::C7]), 2) + pow(abs(Cp[WCoef::CP7]), 2);
762 complex_t W_9 = pow(abs(C[WCoef::C9]), 2) + pow(abs(Cp[WCoef::CP9]), 2);
763 complex_t W_10 = pow(abs(C[WCoef::C10]), 2) + pow(abs(Cp[WCoef::CP10]), 2);
764 complex_t W_27 = (C[WCoef::C2] + C_F * C[WCoef::C1]) * conj(C[WCoef::C7]) + (Cp[WCoef::CP2] + C_F * Cp[WCoef::CP1]) * conj(Cp[WCoef::CP7]);
765 complex_t W_29 = (C[WCoef::C2] + C_F * C[WCoef::C1]) * conj(C[WCoef::C9]) + (Cp[WCoef::CP2] + C_F * Cp[WCoef::CP1]) * conj(Cp[WCoef::CP9]);
766 complex_t W_79 = C[WCoef::C7] * conj(C[WCoef::C9]) + Cp[WCoef::CP7] * conj(Cp[WCoef::CP9]);
767
768 return pow(1 - s, 2) * real(
769 8 * (1 + 2 * s) * (
770 W_9 * omega_99(s, L_l) +
771 W_10 * omega_1010(s, L_l) +
772 real(W_29 * omega_29(s, L_l)) +
773 W_2 * omega_22(s, L_l)
774 ) +
775 96 * real(
776 W_79 * omega_79(s, L_l) +
777 W_27 * omega_27(s, L_l)
778 ) +
779 8 * (4 + 8 / s) * W_7 * omega_77(s, L_l)
780 );
781}
782
783double BXsllDecay::delta_A_em(double s, double L_l) {
784 auto C = cache.C;
785 auto Cp = cache.C;
786
787 double C_F = (*iobs_qcdp).get_constants()->C_F;
788 complex_t W_210 = (C[WCoef::C2] + C_F * C[WCoef::C1]) * conj(C[WCoef::C10]) + (Cp[WCoef::CP2] + C_F * Cp[WCoef::CP1]) * conj(Cp[WCoef::CP10]);
789 complex_t W_710 = C[WCoef::C7] * conj(C[WCoef::C10]) + Cp[WCoef::CP7] * conj(Cp[WCoef::CP10]);
790 complex_t W_910 = C[WCoef::C9] * conj(C[WCoef::C10]) + Cp[WCoef::CP9] * conj(Cp[WCoef::CP10]);
791
792 return pow(1 - s, 2) * real(
793 -48. * W_710 * omega_710(s, L_l) +
794 -24. * s * (
795 W_910 * omega_910(s, L_l) +
796 W_210 * omega_210(s, L_l, cache.z)
797 )
798 );
799}
800
801double BXsllDecay::A_FB(double s, double ml_hat, double L_l) {
802 return -3 * A_FB_0(s, ml_hat) * (cache.pref_A0_0 + cache.pref_A0_1 * s / pow(1 - s, 2)) +
803 cache.pref_delta_mb2 * delta_A_mb2(s) +
804 -3 / 8 * cache.pref_delta_mc2 * delta_A_mc2(s) +
805 8 / 3 * cache.pref_delta_brems * delta_A_brem(s) +
806 cache.pref_delta_em * delta_A_em(s, L_l);
807}
808
809double BXsllDecay::dB_ds(double s, double ml_hat, double L_l) {
810 double dB = cache.pref_dB0_ds * dB0_ds(s, ml_hat) +
811 cache.pref_delta_brems * (delta_bremA(s) + delta_bremB(s)) +
812 cache.pref_delta_em * delta_em(s, L_l);
813
814 if (cfg.gen != BXsllConfig::Lepton::TAU) {
815 dB += cache.pref_delta_mb2 * delta_mb2(s) +
816 cache.pref_delta_mb3 * delta_mb3(s) +
817 cache.pref_delta_mc2 * delta_mc2(s);
818 }
819
820 return dB;
821}
822
823std::vector<ObservableValue> BXsllDecay::BR_B_Xsll(Observables oid) {
824 std::vector<ObservableValue> out;
825
826 auto f = [&] (double s) {
827 return dB_ds(s, cache.m_l_hat, cache.L_l);
828 };
829
830 constexpr double s_eps_low = 1e-8;
831 constexpr double s_eps_high = 1e-3; // ou 1e-2 si encore instable
832
833 for (size_t i = 0; i < this->bins.value().size(); i++) {
834 double rand_err = this->bins.value()[i].second < 8.0 ? cache.rand_err[0] : this->bins.value()[i].first < 12.0 ? cache.rand_err[1] : cache.rand_err[2];
835 const auto [q2_min, q2_max] = this->bins.value()[i];
836
837 double s_min_raw = q2_min / std::pow(cache.m_b_1S, 2);
838 double s_max_raw = q2_max / std::pow(cache.m_b_1S, 2);
839
840 double s_min = std::max(s_min_raw, s_eps_low);
841 double s_max = std::min(s_max_raw, 1.0 - s_eps_high);
842
843 // if (s_max_raw >= 1.0) {
844 // LOG_WARN(
845 // "BXsll bin [", q2_min, ",", q2_max,
846 // "] has s_max = ", s_max_raw,
847 // " >= 1 for m_b_1S = ", cache.m_b_1S,
848 // ". Clipping to ", s_max
849 // );
850 // }
851
852 if (!(s_min < s_max)) {
853 LOG_WARN(
854 "Skipping invalid BXsll bin [", q2_min, ",", q2_max,
855 "] after clipping s from [", s_min_raw, ",", s_max_raw,
856 "] to [", s_min, ",", s_max, "]"
857 );
858
859 out.emplace_back(
861 std::numeric_limits<double>::quiet_NaN(),
862 this->bins.value()[i]
863 );
864 continue;
865 }
866
867 double res = cache.pref_dB_ds * integrate(f, s_min, s_max, 1e-3) * (1 + rand_err);
868
869 out.emplace_back(
871 res,
872 this->bins.value()[i]
873 );
874 }
875
876 return out;
877}
878
879std::vector<ObservableValue> BXsllDecay::A_FB_B_Xsll(Observables oid) {
880 std::vector<ObservableValue> out;
881
882 auto f = [&] (double s) {
883 return A_FB(s, cache.m_l_hat, cache.L_l);
884 };
885
886 constexpr double s_eps_low = 1e-8;
887 constexpr double s_eps_high = 1e-3;
888
889 for (size_t i = 0; i < this->bins.value().size(); i++) {
890 const auto [q2_min, q2_max] = this->bins.value()[i];
891
892 double s_min_raw = q2_min / std::pow(cache.m_b_1S, 2);
893 double s_max_raw = q2_max / std::pow(cache.m_b_1S, 2);
894
895 double s_min = std::max(s_min_raw, s_eps_low);
896 double s_max = std::min(s_max_raw, 1.0 - s_eps_high);
897
898 if (!(s_min < s_max)) {
899 out.emplace_back(
901 std::numeric_limits<double>::quiet_NaN(),
902 this->bins.value()[i]
903 );
904 continue;
905 }
906
907 double res = cache.pref_dB_ds * integrate(f, s_min, s_max, 1e-3);
908
909 out.emplace_back(
911 res,
912 this->bins.value()[i]
913 );
914 }
915
916 return out;
917}
918
919
920std::vector<ObservableValue> BXsllDecay::compute_observable(Observables obs) {
921 switch (obs) {
924 return BR_B_Xsll(obs);
925 break;
928 return BR_B_Xsll(obs);
929 break;
932 return BR_B_Xsll(obs);
933 break;
936 return A_FB_B_Xsll(obs);
937 break;
940 return A_FB_B_Xsll(obs);
941 break;
944 return A_FB_B_Xsll(obs);
945 break;
946 default:
947 LOG_ERROR("IndexError", "Observable", ObservableMapper::str(obs), "doesn't belong to the decay", DecayMapper::str(this->id));
948 }
949}
950
951std::vector<ObservableValue> BXsllDecay::compute_observable(ObservableId obs) {
953}
Observables
Definition GeneralEnum.h:4
QCDOrder
#define LOG_ERROR(type,...)
Macro for logging error messages and terminating the application.
Definition Logger.h:41
#define LOG_WARN(...)
Macro for logging warning messages.
Definition Logger.h:40
void fill_cache(Func &&f, double a, double b, std::array< T, cache_size > &cache, Args &&... args)
Fills a lookup cache for a function on a finite interval [a, b].
Definition Utils.h:241
std::complex< double > complex_t
Convenience alias for std::complex<double>.
Definition Utils.h:35
T lerp(U x, const std::array< T, cache_size > &lookup, double a=0.0, double b=1.0)
Linearly interpolates a cached function on [a, b].
Definition Utils.h:276
complex_t omega_27(double s, double L_l)
complex_t F_19(double s)
double omega_79(double s, double L_l)
double g_lambda(double z)
complex_t W_29(double s)
complex_t G0(double t)
double W_7(double s)
complex_t F_29(double s)
complex_t C7_new(double s, bool prime)
double W_79(double s)
complex_t g_ld(double z, double s)
double delta_bremB(double s)
double delta_mc2(double s)
double omega_910(double s, double L_l)
double f(double z)
double delta_em(double s, double L_l)
double f_9(double s)
void load_cfg_dep_params()
complex_t omega_29(double s, double L_l)
double delta_bremA(double s)
double f_7(double s)
double tau_810(double s)
double dB0_ds(double s, double ml_hat)
double tau_77(double s)
double breit_wigner(double s, double m_V, double br, double gamma_tot, double gamma_had)
double delta_mb2(double s)
double delta_A_em(double s, double L_l)
complex_t omega_210(double s, double L_l, double z)
double kappa(double z)
double W_9(double s)
complex_t C9_eff(double s, QCDOrder order, bool prime)
double sigma_7(double s, double L_mu)
double sigma(double s)
complex_t tau_29(double s, double w, complex_t Delta_23, complex_t Delta_27)
complex_t Sigma_3(double s)
double tau_710(double s)
complex_t F_17(double s)
double PV_R_cc_cont(double s)
complex_t Delta_i_27(double s, double z, double w)
std::vector< ObservableValue > compute_observable(Observables obs) override
Compute an observable given a public observable enum.
void fill_wilson_cache()
complex_t W_27(double s)
complex_t W_210(double s)
double W_710(double s)
double A_FB_0(double s, double ml_hat)
double h(double z)
complex_t F(double r)
void set_cfg_flags(BXsllConfig::Lepton gen)
double delta_A_mb2(double s)
double R_cc_cont(double s)
complex_t tau_210(double s, double z)
double PV_R_cc(double s)
complex_t g(double z, double s)
double dB_ds(double s, double ml_hat, double L_l)
double tau_78(double s)
complex_t Sigma_7(double s, double z)
void load_params() override
Load and cache parameters needed by this decay.
Definition BXsllDecay.cpp:4
double tau_910(double s)
double W_10(double s)
double omega_1010(double s, double L_l)
complex_t Gm1(double t)
double tau_22(double s, double w, complex_t Delta_23, complex_t Delta_27)
complex_t C10_new(double s, bool prime)
double tau_99(double s)
double W_910(double s)
std::vector< ObservableValue > BR_B_Xsll(Observables oid)
double tau_89(double s)
double R_cc(double s)
double omega_99(double s, double L_l)
double tau_88(double s)
double delta_A_mc2(double s)
double delta_mb3(double s)
std::vector< ObservableValue > A_FB_B_Xsll(Observables oid)
complex_t Sigma_1(double s)
complex_t tau_27(double s, double w, complex_t Delta_23, complex_t Delta_27)
double omega_77(double s, double L_l)
double omega_710(double s, double L_l)
double omega_22(double s, double L_l)
double Sigma_2(double s)
double g_rho(double z)
double sigma_9(double s)
complex_t C9_new(double s, bool prime)
double PV_breit_wigner(double s, double m_V, double br, double gamma_tot, double gamma_had)
double tau_79(double s)
complex_t Delta_i_23(double s, double z, double w)
double delta_A_brem(double s)
double A_FB(double s, double ml_hat, double L_l)
double delta_bremB_base(double s)
complex_t tau_28(double s, double w, complex_t Delta_23, complex_t Delta_27)
complex_t F_27(double s)
Definition BWilson.h:149
DecayId id
Unique decay identifier.
WilsonBuildConfig w_config
Wilson build configuration used when enabling this decay (scales, order, groups).
std::optional< std::vector< std::pair< double, double > > > bins
Optional q^2 bins.
std::shared_ptr< IObsParameterProxy< ParamId, DataType, std::string, LhaID > > p
Parameter proxy for SM-like quantities used by the decay (may be SM/BSM depending on wiring).
std::shared_ptr< IObsWilsonProxy > w_proxy
Wilson proxy used at compute-time to query coefficients (matching/run).
Definition DecayParent.h:99
static std::optional< Observables > enum_of(const IdOf< ObservableTag > &id)
Attempts to recover the enum value associated with an identifier.
static IdOf< ObservableTag > to_id(Observables e)
Converts an enum value to an IdOf<Tag>.
static std::string str(const IdOf< ObservableTag > &id)
Returns the string representation of an identifier.
static WCoef cpq1_for_lepton_index(int lepton_index)
static WCoef cq1_for_lepton_index(int lepton_index)
static WCoef cq2_for_lepton_index(int lepton_index)
static std::vector< WCoef > get_group(WGroup g)
Returns the list of Wilson coefficients belonging to a WGroup.
static WCoef cpq2_for_lepton_index(int lepton_index)
constexpr double PI
Definition constants.h:7
constexpr std::complex< double > I
Definition constants.h:20
constexpr double PI2
Definition constants.h:8
constexpr double g
scalar_t c_integrate(ComplexValuedFunction f, double l, double u, double prec)
Performs numerical integration of a complex-valued function of a real variable.
double integrate(RealValuedFunction f, double l, double u, double prec)
Performs numerical integration of a real-valued function of a real variable.
complex_t f_19_1S(double s_hat, double L_b, double z, size_t max_pow=20)
complex_t f_27(double s_hat, double L_b, double z, size_t max_pow=20)
complex_t f_29_1S(double s_hat, double L_b, double z, size_t max_pow=20)
complex_t f_17(double s_hat, double L_b, double z, size_t max_pow=20)
complex_t f_87(double s_hat, double L_b)
complex_t f_89(double s_hat)
scalar_t CLi2(scalar_t x)
Computes the complex dilogarithm function.
Definition polylog.cpp:1312
double Li2(double x)
Computes the dilogarithm function Li2(x).
Definition polylog.cpp:1257
scalar_t pow(const scalar_t &base, const scalar_t &exp)
Definition scalar.cpp:75
double real(const scalar_t &z)
Definition scalar.cpp:83
Configuration for evaluating the strong coupling constant .
Definition Configs.h:206
Lepton gen
Definition BXsllDecay.h:12
std::array< complex_t, 100 > F_29_lookup
Definition BXsllDecay.h:34
double pref_delta_mc2
Definition BXsllDecay.h:24
std::array< double, 100 > delta_brems_lookup
Definition BXsllDecay.h:35
double pref_delta_mb2
Definition BXsllDecay.h:24
std::array< double, 3 > rand_err
Definition BXsllDecay.h:26
std::array< double, 6 > cc_res_width_had
Definition BXsllDecay.h:40
std::array< complex_t, 100 > F_27_lookup
Definition BXsllDecay.h:32
double pref_delta_mb3
Definition BXsllDecay.h:24
std::array< double, 6 > cc_res_width_tot
Definition BXsllDecay.h:39
double pref_dB0_ds
Definition BXsllDecay.h:22
std::array< double, 6 > cc_res_mass
Definition BXsllDecay.h:37
double alpha_s_mu_b
Definition BXsllDecay.h:21
std::array< complex_t, 100 > F_17_lookup
Definition BXsllDecay.h:31
double pref_delta_brems
Definition BXsllDecay.h:24
double pref_dB_ds
Definition BXsllDecay.h:22
std::array< double, 6 > cc_res_br
Definition BXsllDecay.h:38
std::map< WCoef, complex_t > C_LO
Definition BXsllDecay.h:29
std::array< complex_t, 100 > F_19_lookup
Definition BXsllDecay.h:33
double pref_delta_em
Definition BXsllDecay.h:24
std::map< WCoef, complex_t > C
Definition BXsllDecay.h:28
Composite identifier for a single parameter.
Definition ParamID.h:57
QCDOrder order
Perturbative QCD order used for the evolution and matching of Wilson coefficients....
Definition Configs.h:54