Hyperiso 1.0.3
Modular flavour-physics calculations, Wilson coefficients and statistical inference
Loading...
Searching...
No Matches
LbLllDecay.cpp
Go to the documentation of this file.
1#include "LbLllDecay.h"
2
3#include <algorithm>
4#include <array>
5#include <cmath>
6#include <exception>
7#include <limits>
8#include <mutex>
9#include <thread>
10#include <vector>
11
12void LbLllDecay::load_params() {
13 fill_wilson_cache();
14
15 cache.ff_calculator = LbLFFCalculator(p, cfg.ff_src);
16 cache.m_b_mu_b = (*iobs_qcdp)(MassConfig(5, (*p)(ParamId{ParameterType::WILSON, "B_SCALE", 1}, DataType::VALUE), MassType::MSBAR, MassType::POLE));
17 cache.m_Lb = (*p)(ParamId{ParameterType::FLAVOR, "FMASS", 5122}, DataType::VALUE);
18 cache.m_L = (*p)(ParamId{ParameterType::FLAVOR, "FMASS", 3122}, DataType::VALUE);
19 cache.life_L = (*p)(ParamId{ParameterType::FLAVOR, "FLIFE", 5122}, DataType::VALUE) / HBAR;
20 cache.alpha_L = (*p)(ParamId{ParameterType::DECAY, "Lb_L", 8}, DataType::VALUE);
21 cache.N_0 = std::conj((*p)(ParamId{ParameterType::SM, "VCKM", {2, 1}}, DataType::VALUE)) * (*p)(ParamId{ParameterType::SM, "VCKM", {2, 2}}, DataType::VALUE) * (*p)(ParamId{ParameterType::SM, "SMINPUTS", 2}, DataType::VALUE) * (*p)(ParamId{ParameterType::SM, "EW", {1, 2}}, DataType::VALUE) / (std::sqrt(6144. * std::pow(PI, 5) * std::pow(cache.m_Lb, 3)));
22 cache.q2_max = std::pow(cache.m_Lb - cache.m_L, 2);
23
24 load_cfg_dep_params();
25
26 // printf("alpha_em = %.4e\n", (*p)(ParamId{ParameterType::SM, "EW", {1, 2}}, DataType::VALUE).real());
27 // printf("N0 = %.4e + %.4e i\n", cache.N_0.real(), cache.N_0.imag());
28
29 // printf("f_perp = %.4e\n", cache.ff_calculator.get(LbL_FF::F_PERP, 1.0));
30 // printf("h_perp = %.4e\n", cache.ff_calculator.get(LbL_FF::H_PERP, 1.0));
31 // printf("g_perp = %.4e\n", cache.ff_calculator.get(LbL_FF::G_PERP, 1.0));
32 // printf("f_+ = %.4e\n", cache.ff_calculator.get(LbL_FF::F_PLUS, 1.0));
33 // printf("h_+ = %.4e\n", cache.ff_calculator.get(LbL_FF::H_PLUS, 1.0));
34 // printf("g_+ = %.4e\n", cache.ff_calculator.get(LbL_FF::G_PLUS, 1.0));
35 // printf("h_tilde_perp = %.4e\n", cache.ff_calculator.get(LbL_FF::H_TILDE_PERP, 1.0));
36 // printf("h_tilde_+ = %.4e\n", cache.ff_calculator.get(LbL_FF::H_TILDE_PLUS, 1.0));
37}
38
39void LbLllDecay::fill_wilson_cache() {
40 cache.C.clear();
41
42 auto b_wilsons = w_proxy->getAFR(WGroup::B, this->w_config.order);
43 auto bp_wilsons = w_proxy->getAFR(WGroup::BPrime, this->w_config.order);
44
45 cache.C[WCoef::C7] = b_wilsons.at(WCoef::C7);
46 cache.C[WCoef::C9] = b_wilsons.at(WCoef::C9);
47 cache.C[WCoef::C10] = b_wilsons.at(WCoef::C10);
48
49 cache.C[WCoef::CP7] = bp_wilsons.at(WCoef::CP7);
50 cache.C[WCoef::CP9] = bp_wilsons.at(WCoef::CP9);
51 cache.C[WCoef::CP10] = bp_wilsons.at(WCoef::CP10);
52}
53
54void LbLllDecay::set_cfg_flags(LbLllConfig::Lepton gen) {
55 if (cfg.gen != gen) {
56 cfg.gen = gen;
57 load_cfg_dep_params();
58 }
59}
60
61void LbLllDecay::load_cfg_dep_params() {
62 cache.m_l = (*p)(ParamId{ParameterType::SM, "MASS", 11 + 2 * (int)cfg.gen}, DataType::VALUE);
63 cache.q2_min = 4 * std::pow(cache.m_l, 2);
64 compute_binned_K_i();
65}
66
67double LbLllDecay::beta_l(double q2) {
68 return std::sqrt(1 - 4. * cache.m_l * cache.m_l / q2);
69}
70
71double LbLllDecay::lambda(double q2) {
72 double mLb2 = cache.m_Lb * cache.m_Lb;
73 double mL2 = cache.m_L * cache.m_L;
74 return mLb2 * mLb2 + mL2 * mL2 + q2 * q2 - 2. * (mLb2 * mL2 + (mLb2 + mL2) * q2);
75}
76
77double LbLllDecay::s_p(double q2) {
78 return std::pow(cache.m_Lb + cache.m_L, 2) - q2;
79}
80
81double LbLllDecay::s_m(double q2) {
82 return std::pow(cache.m_Lb - cache.m_L, 2) - q2;
83}
84
85complex_t LbLllDecay::N(double q2, bool bar) {
86 complex_t N0 = bar ? std::conj(cache.N_0) : cache.N_0;
87 return N0 * std::sqrt(q2 * beta_l(q2) * std::sqrt(lambda(q2)));
88}
89
90complex_t LbLllDecay::A_perp_1(double q2, double sign, bool bar) {
91 complex_t C_V = (cache.C[WCoef::C9] + cache.C[WCoef::CP9]) + sign * (cache.C[WCoef::C10] + cache.C[WCoef::CP10]);
92 complex_t C_T = cache.C[WCoef::C7] + cache.C[WCoef::CP7];
93 if (bar) {
94 C_V = std::conj(C_V);
95 C_T = std::conj(C_T);
96 }
97 complex_t HVplus = -cache.ff_calculator.get(LbL_FF::F_PERP, q2) * std::sqrt(2 * s_m(q2));
98 complex_t HTplus = cache.ff_calculator.get(LbL_FF::H_PERP, q2) * (cache.m_Lb + cache.m_L) * std::sqrt(2 * s_m(q2));
99 return RT2 * N(q2, bar) * (C_V * HVplus - 2 * cache.m_b_mu_b / q2 * C_T * HTplus);
100}
101
102complex_t LbLllDecay::A_par_1(double q2, double sign, bool bar) {
103 complex_t C_A = (cache.C[WCoef::C9] - cache.C[WCoef::CP9]) + sign * (cache.C[WCoef::C10] - cache.C[WCoef::CP10]);
104 complex_t C_TA = cache.C[WCoef::C7] - cache.C[WCoef::CP7];
105 if (bar) {
106 C_A = std::conj(C_A);
107 C_TA = std::conj(C_TA);
108 }
109 complex_t HAplus = -cache.ff_calculator.get(LbL_FF::G_PERP, q2) * std::sqrt(2 * s_p(q2));
110 complex_t HT5plus = -cache.ff_calculator.get(LbL_FF::H_TILDE_PERP, q2) * (cache.m_Lb - cache.m_L) * std::sqrt(2 * s_p(q2));
111 return -RT2 * N(q2, bar) * (C_A * HAplus + 2 * cache.m_b_mu_b / q2 * C_TA * HT5plus);
112}
113
114complex_t LbLllDecay::A_perp_0(double q2, double sign, bool bar) {
115 complex_t C_V = (cache.C[WCoef::C9] + cache.C[WCoef::CP9]) + sign * (cache.C[WCoef::C10] + cache.C[WCoef::CP10]);
116 complex_t C_T = cache.C[WCoef::C7] + cache.C[WCoef::CP7];
117 if (bar) {
118 C_V = std::conj(C_V);
119 C_T = std::conj(C_T);
120 }
121 complex_t HV0 = cache.ff_calculator.get(LbL_FF::F_PLUS, q2) * (cache.m_Lb + cache.m_L) * std::sqrt(s_m(q2) / q2);
122 complex_t HT0 = -cache.ff_calculator.get(LbL_FF::H_PLUS, q2) * std::sqrt(q2 * s_m(q2));
123 return RT2 * N(q2, bar) * (C_V * HV0 - 2 * cache.m_b_mu_b / q2 * C_T * HT0);
124}
125
126complex_t LbLllDecay::A_par_0(double q2, double sign, bool bar) {
127 complex_t C_A = (cache.C[WCoef::C9] - cache.C[WCoef::CP9]) + sign * (cache.C[WCoef::C10] - cache.C[WCoef::CP10]);
128 complex_t C_TA = cache.C[WCoef::C7] - cache.C[WCoef::CP7];
129 if (bar) {
130 C_A = std::conj(C_A);
131 C_TA = std::conj(C_TA);
132 }
133 complex_t HA0 = cache.ff_calculator.get(LbL_FF::G_PLUS, q2) * (cache.m_Lb - cache.m_L) * std::sqrt(s_p(q2) / q2);
134 complex_t HT50 = cache.ff_calculator.get(LbL_FF::H_TILDE_PLUS, q2) * std::sqrt(q2 * s_p(q2));
135 return -RT2 * N(q2, bar) * (C_A * HA0 + 2 * cache.m_b_mu_b / q2 * C_TA * HT50);
136}
137
138double LbLllDecay::K1ss(double q2, bool bar) {
139 complex_t ALpar1 = A_par_1(q2, -1, bar);
140 complex_t ALperp1 = A_perp_1(q2, -1, bar);
141 complex_t ARpar1 = A_par_1(q2, 1, bar);
142 complex_t ARperp1 = A_perp_1(q2, 1, bar);
143 complex_t ALpar0 = A_par_0(q2, -1, bar);
144 complex_t ALperp0 = A_perp_0(q2, -1, bar);
145 complex_t ARpar0 = A_par_0(q2, 1, bar);
146 complex_t ARperp0 = A_perp_0(q2, 1, bar);
147
148 double b_l = beta_l(q2);
149
150 return std::real(0.25*(ALpar1*conj(ALpar1)+ALperp1*conj(ALperp1)+ARpar1*conj(ARpar1)+ARperp1*conj(ARperp1))
151 + 0.25*(1.+b_l*b_l)*(ALpar0*conj(ALpar0)+ALperp0*conj(ALperp0)+ARpar0*conj(ARpar0)+ARperp0*conj(ARperp0))
152 + 0.5*(1.-b_l*b_l)*(ARpar1*conj(ALpar1)+ARperp1*conj(ALperp1)+ARpar0*conj(ALpar0)+ARperp0*conj(ALperp0)));
153}
154
155double LbLllDecay::K1cc(double q2, bool bar) {
156 complex_t ALpar1 = A_par_1(q2, -1, bar);
157 complex_t ALperp1 = A_perp_1(q2, -1, bar);
158 complex_t ARpar1 = A_par_1(q2, 1, bar);
159 complex_t ARperp1 = A_perp_1(q2, 1, bar);
160 complex_t ALpar0 = A_par_0(q2, -1, bar);
161 complex_t ALperp0 = A_perp_0(q2, -1, bar);
162 complex_t ARpar0 = A_par_0(q2, 1, bar);
163 complex_t ARperp0 = A_perp_0(q2, 1, bar);
164
165 double b_l = beta_l(q2);
166
167 return std::real(0.25*(1.+b_l*b_l)*(ARpar1*conj(ARpar1)+ARperp1*conj(ARperp1)+ALpar1*conj(ALpar1)+ALperp1*conj(ALperp1))
168 +0.25*(1.-b_l*b_l)*(ARpar0*conj(ARpar0)+ARperp0*conj(ARperp0)+ALpar0*conj(ALpar0)+ALperp0*conj(ALperp0))
169 +0.5*(1.-b_l*b_l)*(ARpar1*conj(ALpar1)+ARperp1*conj(ALperp1)+ARpar0*conj(ALpar0)+ARperp0*conj(ALperp0)));
170}
171
172double LbLllDecay::K1c(double q2, bool bar) {
173 complex_t ALpar1 = A_par_1(q2, -1, bar);
174 complex_t ALperp1 = A_perp_1(q2, -1, bar);
175 complex_t ARpar1 = A_par_1(q2, 1, bar);
176 complex_t ARperp1 = A_perp_1(q2, 1, bar);
177 double b_l = beta_l(q2);
178 return -b_l*std::real(ARperp1*conj(ARpar1)-ALperp1*conj(ALpar1));;
179}
180
181double LbLllDecay::K2ss(double q2, bool bar) {
182 complex_t ALpar1 = A_par_1(q2, -1, bar);
183 complex_t ALperp1 = A_perp_1(q2, -1, bar);
184 complex_t ARpar1 = A_par_1(q2, 1, bar);
185 complex_t ARperp1 = A_perp_1(q2, 1, bar);
186 complex_t ALpar0 = A_par_0(q2, -1, bar);
187 complex_t ALperp0 = A_perp_0(q2, -1, bar);
188 complex_t ARpar0 = A_par_0(q2, 1, bar);
189 complex_t ARperp0 = A_perp_0(q2, 1, bar);
190
191 double b_l = beta_l(q2);
192
193 return 0.5 * cache.alpha_L * std::real(
194 ARperp1*conj(ARpar1)+ALperp1*conj(ALpar1)
195 + (1. + b_l*b_l)*(ARperp0*conj(ARpar0)+ALperp0*conj(ALpar0))
196 + (1. - b_l*b_l)*(ARperp1*conj(ALpar1)+ARpar1*conj(ALperp1)+ARperp0*conj(ALpar0)+ARpar0*conj(ALperp0)));
197}
198
199double LbLllDecay::K2cc(double q2, bool bar) {
200 complex_t ALpar1 = A_par_1(q2, -1, bar);
201 complex_t ALperp1 = A_perp_1(q2, -1, bar);
202 complex_t ARpar1 = A_par_1(q2, 1, bar);
203 complex_t ARperp1 = A_perp_1(q2, 1, bar);
204 complex_t ALpar0 = A_par_0(q2, -1, bar);
205 complex_t ALperp0 = A_perp_0(q2, -1, bar);
206 complex_t ARpar0 = A_par_0(q2, 1, bar);
207 complex_t ARperp0 = A_perp_0(q2, 1, bar);
208
209 double b_l = beta_l(q2);
210
211 return 0.5 * cache.alpha_L * std::real(
212 (1.+b_l*b_l)*(ARperp1*conj(ARpar1)+ALperp1*conj(ALpar1))
213 + (1.-b_l*b_l)*(ARpar0*conj(ARperp0)+ALpar0*conj(ALperp0))
214 + (1.-b_l*b_l)*(ARperp1*conj(ALpar1)+ARpar1*conj(ALperp1)+ARperp0*conj(ALpar0)+ARpar0*conj(ALperp0)));
215}
216
217double LbLllDecay::K2c(double q2, bool bar) {
218 complex_t ALpar1 = A_par_1(q2, -1, bar);
219 complex_t ALperp1 = A_perp_1(q2, -1, bar);
220 complex_t ARpar1 = A_par_1(q2, 1, bar);
221 complex_t ARperp1 = A_perp_1(q2, 1, bar);
222 double b_l = beta_l(q2);
223 return -0.5*cache.alpha_L*b_l*std::real(ARpar1*conj(ARpar1)+ARperp1*conj(ARperp1)-ALpar1*conj(ALpar1)-ALperp1*conj(ALperp1));
224}
225
226void LbLllDecay::compute_binned_K_i() {
227 // Same strategy as BKstarllDecay::compute_binned_J_i:
228 // - one quadrature pass per bin and CP-conjugation
229 // - evaluate the 8 transversity amplitudes once per q2 point
230 // - derive all K_i from those amplitudes
231 // - split independent bins across cfg.n_threads workers
232 static constexpr std::array<double, 24> GL24_X {{
233 -0.99518721999702131, -0.97472855597130947, -0.93827455200273280, -0.88641552700440107,
234 -0.82000198597390295, -0.74012419157855436, -0.64809365193697555, -0.54542147138883956,
235 -0.43379350762604513, -0.31504267969616340, -0.19111886747361631, -0.06405689286260563,
236 0.06405689286260563, 0.19111886747361631, 0.31504267969616340, 0.43379350762604513,
237 0.54542147138883956, 0.64809365193697555, 0.74012419157855436, 0.82000198597390295,
238 0.88641552700440107, 0.93827455200273280, 0.97472855597130947, 0.99518721999702131
239 }};
240 static constexpr std::array<double, 24> GL24_W {{
241 0.01234122979998869, 0.02853138862893356, 0.04427743881741941, 0.05929858491543636,
242 0.07334648141108016, 0.08619016153195321, 0.09761865210411393, 0.10744427011596556,
243 0.11550566805372552, 0.12167047292780329, 0.12583745634682825, 0.12793819534675202,
244 0.12793819534675202, 0.12583745634682825, 0.12167047292780329, 0.11550566805372552,
245 0.10744427011596556, 0.09761865210411393, 0.08619016153195321, 0.07334648141108016,
246 0.05929858491543636, 0.04427743881741941, 0.02853138862893356, 0.01234122979998869
247 }};
248
249 struct AmpSet {
250 complex_t ALpar1;
251 complex_t ALperp1;
252 complex_t ARpar1;
253 complex_t ARperp1;
254 complex_t ALpar0;
255 complex_t ALperp0;
256 complex_t ARpar0;
257 complex_t ARperp0;
258 double beta;
259 double beta2;
260 };
261
262 const auto& bins = this->bins.value();
263 const size_t nbins = bins.size();
264
265 auto clear_and_resize = [&] (std::array<std::vector<double>, 6>& dest) {
266 for (auto& v : dest) {
267 v.assign(nbins, std::numeric_limits<double>::quiet_NaN());
268 }
269 };
270
271 clear_and_resize(cache.K_i_binned);
272 clear_and_resize(cache.K_i_bar_binned);
273 cache.bin_widths.assign(nbins, std::numeric_limits<double>::quiet_NaN());
274
275 for (size_t i = 0; i < nbins; ++i) {
276 cache.bin_widths[i] = bins[i].second - bins[i].first;
277 }
278
279 if (nbins == 0) {
280 return;
281 }
282
283 auto eval_amplitudes = [&] (double q2, bool bar) -> AmpSet {
284 const double beta = beta_l(q2);
285 return {
286 A_par_1(q2, -1, bar),
287 A_perp_1(q2, -1, bar),
288 A_par_1(q2, 1, bar),
289 A_perp_1(q2, 1, bar),
290 A_par_0(q2, -1, bar),
291 A_perp_0(q2, -1, bar),
292 A_par_0(q2, 1, bar),
293 A_perp_0(q2, 1, bar),
294 beta,
295 beta * beta
296 };
297 };
298
299 auto eval_integrands = [&] (double q2, bool bar) -> std::array<double, 6> {
300 const AmpSet a = eval_amplitudes(q2, bar);
301
302 const double norm_ALpar1 = std::norm(a.ALpar1);
303 const double norm_ALperp1 = std::norm(a.ALperp1);
304 const double norm_ARpar1 = std::norm(a.ARpar1);
305 const double norm_ARperp1 = std::norm(a.ARperp1);
306 const double norm_ALpar0 = std::norm(a.ALpar0);
307 const double norm_ALperp0 = std::norm(a.ALperp0);
308 const double norm_ARpar0 = std::norm(a.ARpar0);
309 const double norm_ARperp0 = std::norm(a.ARperp0);
310
311 const double one_plus_beta2 = 1.0 + a.beta2;
312 const double one_minus_beta2 = 1.0 - a.beta2;
313 const double alpha_L = cache.alpha_L;
314
315 const complex_t k1_cross =
316 a.ARpar1 * std::conj(a.ALpar1)
317 + a.ARperp1 * std::conj(a.ALperp1)
318 + a.ARpar0 * std::conj(a.ALpar0)
319 + a.ARperp0 * std::conj(a.ALperp0);
320
321 const complex_t k2_cross =
322 a.ARperp1 * std::conj(a.ALpar1)
323 + a.ARpar1 * std::conj(a.ALperp1)
324 + a.ARperp0 * std::conj(a.ALpar0)
325 + a.ARpar0 * std::conj(a.ALperp0);
326
327 const double k1ss = std::real(
328 0.25 * (norm_ALpar1 + norm_ALperp1 + norm_ARpar1 + norm_ARperp1)
329 + 0.25 * one_plus_beta2 * (norm_ALpar0 + norm_ALperp0 + norm_ARpar0 + norm_ARperp0)
330 + 0.5 * one_minus_beta2 * k1_cross
331 );
332
333 const double k1cc = std::real(
334 0.25 * one_plus_beta2 * (norm_ARpar1 + norm_ARperp1 + norm_ALpar1 + norm_ALperp1)
335 + 0.25 * one_minus_beta2 * (norm_ARpar0 + norm_ARperp0 + norm_ALpar0 + norm_ALperp0)
336 + 0.5 * one_minus_beta2 * k1_cross
337 );
338
339 const double k1c = -a.beta * std::real(
340 a.ARperp1 * std::conj(a.ARpar1)
341 - a.ALperp1 * std::conj(a.ALpar1)
342 );
343
344 const double k2ss = 0.5 * alpha_L * std::real(
345 a.ARperp1 * std::conj(a.ARpar1)
346 + a.ALperp1 * std::conj(a.ALpar1)
347 + one_plus_beta2 * (a.ARperp0 * std::conj(a.ARpar0) + a.ALperp0 * std::conj(a.ALpar0))
348 + one_minus_beta2 * k2_cross
349 );
350
351 const double k2cc = 0.5 * alpha_L * std::real(
352 one_plus_beta2 * (a.ARperp1 * std::conj(a.ARpar1) + a.ALperp1 * std::conj(a.ALpar1))
353 + one_minus_beta2 * (a.ARpar0 * std::conj(a.ARperp0) + a.ALpar0 * std::conj(a.ALperp0))
354 + one_minus_beta2 * k2_cross
355 );
356
357 const double k2c = -0.5 * alpha_L * a.beta * (
358 norm_ARpar1 + norm_ARperp1 - norm_ALpar1 - norm_ALperp1
359 );
360
361 return {{k1ss, k1cc, k1c, k2ss, k2cc, k2c}};
362 };
363
364 auto integrate_bin = [&] (double q2_l, double q2_u, bool bar) -> std::array<double, 6> {
365 std::array<double, 6> acc {};
366
367 const double width = q2_u - q2_l;
368 if (!std::isfinite(width) || width <= 0.0) {
369 acc.fill(std::numeric_limits<double>::quiet_NaN());
370 return acc;
371 }
372
373 const double center = 0.5 * (q2_l + q2_u);
374 const double half_width = 0.5 * width;
375
376 for (size_t i = 0; i < GL24_X.size(); ++i) {
377 const double q2 = center + half_width * GL24_X[i];
378 const auto vals = eval_integrands(q2, bar);
379 for (size_t k = 0; k < acc.size(); ++k) {
380 acc[k] += GL24_W[i] * vals[k];
381 }
382 }
383
384 for (double& v : acc) {
385 v *= half_width;
386 }
387
388 return acc;
389 };
390
391 auto compute_and_store_bin = [&] (size_t bin_idx) {
392 const auto& [q2_l, q2_u] = bins[bin_idx];
393
394 const auto integ = integrate_bin(q2_l, q2_u, false);
395 const auto integ_bar = integrate_bin(q2_l, q2_u, true);
396
397 for (size_t k = 0; k < cache.K_i_binned.size(); ++k) {
398 cache.K_i_binned[k][bin_idx] = integ[k];
399 cache.K_i_bar_binned[k][bin_idx] = integ_bar[k];
400 }
401 };
402
403 auto requested_threads = cfg.n_threads;
404 if (requested_threads == 0u) {
405 requested_threads = std::thread::hardware_concurrency();
406 }
407 if (requested_threads == 0u) {
408 requested_threads = 1u;
409 }
410
411 const size_t nworkers = std::min<size_t>(requested_threads, nbins);
412
413 if (nworkers <= 1u) {
414 for (size_t i = 0; i < nbins; ++i) {
415 compute_and_store_bin(i);
416 }
417 return;
418 }
419
420 std::vector<std::thread> workers;
421 workers.reserve(nworkers);
422
423 std::exception_ptr first_exception = nullptr;
424 std::mutex exception_mutex;
425
426 auto worker = [&] (size_t begin, size_t end) {
427 try {
428 for (size_t i = begin; i < end; ++i) {
429 compute_and_store_bin(i);
430 }
431 } catch (...) {
432 std::lock_guard<std::mutex> lock(exception_mutex);
433 if (!first_exception) {
434 first_exception = std::current_exception();
435 }
436 }
437 };
438
439 const size_t chunk = (nbins + nworkers - 1) / nworkers;
440 for (size_t w = 0; w < nworkers; ++w) {
441 const size_t begin = w * chunk;
442 const size_t end = std::min(nbins, begin + chunk);
443 if (begin >= end) {
444 break;
445 }
446 workers.emplace_back(worker, begin, end);
447 }
448
449 for (auto& th : workers) {
450 th.join();
451 }
452
453 if (first_exception) {
454 std::rethrow_exception(first_exception);
455 }
456}
457
458std::vector<ObservableValue> LbLllDecay::dBR_dq2_binned(Observables oid, bool br) {
459 std::vector<ObservableValue> out;
460 double br_factor = br ? cache.life_L : 1.0;
461
462 for (size_t i = 0; i < this->bins.value().size(); i++) {
463 double K1ss = cache.K_i_binned[0][i] + cache.K_i_bar_binned[0][i];
464 double K1cc = cache.K_i_binned[1][i] + cache.K_i_bar_binned[1][i];
465 double integrated_rate = (2 * K1ss + K1cc) / 2;
466 const double requested_width = this->bins.value()[i].second - this->bins.value()[i].first;
467 const double width =
468 (i < cache.bin_widths.size() && std::isfinite(cache.bin_widths[i]) && cache.bin_widths[i] > 0.0)
469 ? cache.bin_widths[i]
470 : requested_width;
471
472 const double res =
473 (std::isfinite(width) && width > 0.0)
474 ? integrated_rate * br_factor / width
475 : std::numeric_limits<double>::quiet_NaN();
476
477 out.emplace_back(ObservableMapper::to_id(oid), res, this->bins.value()[i]);
478 }
479
480 return out;
481}
482
483double LbLllDecay::dG_dq2_avg_bin(size_t bin) {
484 double K1ss = cache.K_i_binned[0][bin] + cache.K_i_bar_binned[0][bin];
485 double K1cc = cache.K_i_binned[1][bin] + cache.K_i_bar_binned[1][bin];
486 return (2 * K1ss + K1cc) / 2;
487}
488
489std::vector<ObservableValue> LbLllDecay::A_FB_l(Observables oid) {
490 std::vector<ObservableValue> out;
491 for (size_t i = 0; i < this->bins.value().size(); i++) {
492 double K1c = cache.K_i_binned[2][i] + cache.K_i_bar_binned[2][i];
493 double res = 0.75 * K1c / dG_dq2_avg_bin(i);
494 out.emplace_back(ObservableMapper::to_id(oid), res, this->bins.value()[i]);
495 }
496 return out;
497}
498
499std::vector<ObservableValue> LbLllDecay::A_FB_h(Observables oid) {
500 std::vector<ObservableValue> out;
501 for (size_t i = 0; i < this->bins.value().size(); i++) {
502 double K2ss = cache.K_i_binned[3][i] + cache.K_i_bar_binned[3][i];
503 double K2cc = cache.K_i_binned[4][i] + cache.K_i_bar_binned[4][i];
504 double res = 0.25 * (2 * K2ss + K2cc) / dG_dq2_avg_bin(i);
505 out.emplace_back(ObservableMapper::to_id(oid), res, this->bins.value()[i]);
506 }
507 return out;
508}
509
510std::vector<ObservableValue> LbLllDecay::A_FB_lh(Observables oid) {
511 std::vector<ObservableValue> out;
512 for (size_t i = 0; i < this->bins.value().size(); i++) {
513 double K2c = cache.K_i_binned[5][i] + cache.K_i_bar_binned[5][i];
514 double res = 0.375 * K2c / dG_dq2_avg_bin(i);
515 out.emplace_back(ObservableMapper::to_id(oid), res, this->bins.value()[i]);
516 }
517 return out;
518}
519
520std::vector<ObservableValue> LbLllDecay::F_L(Observables oid) {
521 std::vector<ObservableValue> out;
522 for (size_t i = 0; i < this->bins.value().size(); i++) {
523 double K1ss = cache.K_i_binned[0][i] + cache.K_i_bar_binned[0][i];
524 double K1cc = cache.K_i_binned[1][i] + cache.K_i_bar_binned[1][i];
525 double res = 0.5 * (2 * K1ss - K1cc) / dG_dq2_avg_bin(i);
526 out.emplace_back(ObservableMapper::to_id(oid), res, this->bins.value()[i]);
527 }
528 return out;
529}
530
531std::vector<ObservableValue> LbLllDecay::F_T(Observables oid) {
532 std::vector<ObservableValue> out;
533 for (size_t i = 0; i < this->bins.value().size(); i++) {
534 double K1cc = cache.K_i_binned[1][i] + cache.K_i_bar_binned[1][i];
535 double res = K1cc / dG_dq2_avg_bin(i);
536 out.emplace_back(ObservableMapper::to_id(oid), res, this->bins.value()[i]);
537 }
538 return out;
539}
540
541std::vector<ObservableValue> LbLllDecay::compute_observable(Observables obs) {
542 switch (obs) {
543 case Observables::DBR_DQ2_LAMBDA_B__LAMBDA_E_E:
544 set_cfg_flags(LbLllConfig::Lepton::E);
545 return dBR_dq2_binned(obs, true);
546 case Observables::DGAMMA_DQ2_LAMBDA_B__LAMBDA_E_E:
547 set_cfg_flags(LbLllConfig::Lepton::E);
548 return dBR_dq2_binned(obs, false);
549 case Observables::A_FB_L_LAMBDA_B__LAMBDA_E_E:
550 set_cfg_flags(LbLllConfig::Lepton::E);
551 return A_FB_l(obs);
552 case Observables::A_FB_H_LAMBDA_B__LAMBDA_E_E:
553 set_cfg_flags(LbLllConfig::Lepton::E);
554 return A_FB_h(obs);
555 case Observables::A_FB_LH_LAMBDA_B__LAMBDA_E_E:
556 set_cfg_flags(LbLllConfig::Lepton::E);
557 return A_FB_lh(obs);
558 case Observables::F_L_LAMBDA_B__LAMBDA_E_E:
559 set_cfg_flags(LbLllConfig::Lepton::E);
560 return F_L(obs);
561 case Observables::F_T_LAMBDA_B__LAMBDA_E_E:
562 set_cfg_flags(LbLllConfig::Lepton::E);
563 return F_T(obs);
564 case Observables::DBR_DQ2_LAMBDA_B__LAMBDA_MU_MU:
565 set_cfg_flags(LbLllConfig::Lepton::MU);
566 return dBR_dq2_binned(obs, true);
567 case Observables::DGAMMA_DQ2_LAMBDA_B__LAMBDA_MU_MU:
568 set_cfg_flags(LbLllConfig::Lepton::MU);
569 return dBR_dq2_binned(obs, false);
570 case Observables::A_FB_L_LAMBDA_B__LAMBDA_MU_MU:
571 set_cfg_flags(LbLllConfig::Lepton::MU);
572 return A_FB_l(obs);
573 case Observables::A_FB_H_LAMBDA_B__LAMBDA_MU_MU:
574 set_cfg_flags(LbLllConfig::Lepton::MU);
575 return A_FB_h(obs);
576 case Observables::A_FB_LH_LAMBDA_B__LAMBDA_MU_MU:
577 set_cfg_flags(LbLllConfig::Lepton::MU);
578 return A_FB_lh(obs);
579 case Observables::F_L_LAMBDA_B__LAMBDA_MU_MU:
580 set_cfg_flags(LbLllConfig::Lepton::MU);
581 return F_L(obs);
582 case Observables::F_T_LAMBDA_B__LAMBDA_MU_MU:
583 set_cfg_flags(LbLllConfig::Lepton::MU);
584 return F_T(obs);
585 case Observables::DBR_DQ2_LAMBDA_B__LAMBDA_TAU_TAU:
586 set_cfg_flags(LbLllConfig::Lepton::TAU);
587 return dBR_dq2_binned(obs, true);
588 case Observables::DGAMMA_DQ2_LAMBDA_B__LAMBDA_TAU_TAU:
589 set_cfg_flags(LbLllConfig::Lepton::TAU);
590 return dBR_dq2_binned(obs, false);
591 case Observables::A_FB_L_LAMBDA_B__LAMBDA_TAU_TAU:
592 set_cfg_flags(LbLllConfig::Lepton::TAU);
593 return A_FB_l(obs);
594 case Observables::A_FB_H_LAMBDA_B__LAMBDA_TAU_TAU:
595 set_cfg_flags(LbLllConfig::Lepton::TAU);
596 return A_FB_h(obs);
597 case Observables::A_FB_LH_LAMBDA_B__LAMBDA_TAU_TAU:
598 set_cfg_flags(LbLllConfig::Lepton::TAU);
599 return A_FB_lh(obs);
600 case Observables::F_L_LAMBDA_B__LAMBDA_TAU_TAU:
601 set_cfg_flags(LbLllConfig::Lepton::TAU);
602 return F_L(obs);
603 case Observables::F_T_LAMBDA_B__LAMBDA_TAU_TAU:
604 set_cfg_flags(LbLllConfig::Lepton::TAU);
605 return F_T(obs);
606 default:
607 LOG_ERROR("IndexError", "Observable", ObservableMapper::str(obs), "doesn't belong to the decay", DecayMapper::str(this->id));
608 }
609}
610
611std::vector<ObservableValue> LbLllDecay::compute_observable(ObservableId obs) {
612 return compute_observable(ObservableMapper::enum_of(obs).value());
613}
614
615
616void LbLllDecay::set_n_threads(size_t n_threads) {
617 unsigned int available_threads = std::thread::hardware_concurrency();
618
619 if (available_threads == 0) {
620 available_threads = 1;
621 }
622
623 if (n_threads == 0) {
624 this->cfg.n_threads = available_threads;
625 return;
626 }
627
628 if (n_threads > available_threads) {
629 LOG_WARN(
630 "Requested", n_threads,
631 "threads, but only", available_threads,
632 "are available. Using", available_threads,
633 "threads instead."
634 );
635
636 this->cfg.n_threads = available_threads;
637 return;
638 }
639
640 this->cfg.n_threads = std::max<size_t>(1, n_threads);
641}
Observables
Observables
Definition GeneralEnum.h:4
Observables::A_FB_L_LAMBDA_B__LAMBDA_MU_MU
@ A_FB_L_LAMBDA_B__LAMBDA_MU_MU
Observables::DGAMMA_DQ2_LAMBDA_B__LAMBDA_TAU_TAU
@ DGAMMA_DQ2_LAMBDA_B__LAMBDA_TAU_TAU
Observables::F_T_LAMBDA_B__LAMBDA_MU_MU
@ F_T_LAMBDA_B__LAMBDA_MU_MU
Observables::A_FB_H_LAMBDA_B__LAMBDA_E_E
@ A_FB_H_LAMBDA_B__LAMBDA_E_E
Observables::A_FB_LH_LAMBDA_B__LAMBDA_E_E
@ A_FB_LH_LAMBDA_B__LAMBDA_E_E
Observables::F_L_LAMBDA_B__LAMBDA_MU_MU
@ F_L_LAMBDA_B__LAMBDA_MU_MU
Observables::A_FB_LH_LAMBDA_B__LAMBDA_MU_MU
@ A_FB_LH_LAMBDA_B__LAMBDA_MU_MU
Observables::F_L_LAMBDA_B__LAMBDA_TAU_TAU
@ F_L_LAMBDA_B__LAMBDA_TAU_TAU
Observables::DGAMMA_DQ2_LAMBDA_B__LAMBDA_MU_MU
@ DGAMMA_DQ2_LAMBDA_B__LAMBDA_MU_MU
Observables::DBR_DQ2_LAMBDA_B__LAMBDA_TAU_TAU
@ DBR_DQ2_LAMBDA_B__LAMBDA_TAU_TAU
Observables::F_L_LAMBDA_B__LAMBDA_E_E
@ F_L_LAMBDA_B__LAMBDA_E_E
Observables::F_T_LAMBDA_B__LAMBDA_TAU_TAU
@ F_T_LAMBDA_B__LAMBDA_TAU_TAU
Observables::A_FB_L_LAMBDA_B__LAMBDA_E_E
@ A_FB_L_LAMBDA_B__LAMBDA_E_E
Observables::A_FB_L_LAMBDA_B__LAMBDA_TAU_TAU
@ A_FB_L_LAMBDA_B__LAMBDA_TAU_TAU
Observables::DGAMMA_DQ2_LAMBDA_B__LAMBDA_E_E
@ DGAMMA_DQ2_LAMBDA_B__LAMBDA_E_E
Observables::DBR_DQ2_LAMBDA_B__LAMBDA_E_E
@ DBR_DQ2_LAMBDA_B__LAMBDA_E_E
Observables::DBR_DQ2_LAMBDA_B__LAMBDA_MU_MU
@ DBR_DQ2_LAMBDA_B__LAMBDA_MU_MU
Observables::F_T_LAMBDA_B__LAMBDA_E_E
@ F_T_LAMBDA_B__LAMBDA_E_E
Observables::A_FB_LH_LAMBDA_B__LAMBDA_TAU_TAU
@ A_FB_LH_LAMBDA_B__LAMBDA_TAU_TAU
Observables::A_FB_H_LAMBDA_B__LAMBDA_TAU_TAU
@ A_FB_H_LAMBDA_B__LAMBDA_TAU_TAU
Observables::A_FB_H_LAMBDA_B__LAMBDA_MU_MU
@ A_FB_H_LAMBDA_B__LAMBDA_MU_MU
LbL_FF::H_TILDE_PERP
@ H_TILDE_PERP
LbL_FF::H_TILDE_PLUS
@ H_TILDE_PLUS
LOG_ERROR
#define LOG_ERROR(type,...)
Macro for logging error messages and terminating the application.
Definition Logger.h:41
LOG_WARN
#define LOG_WARN(...)
Macro for logging warning messages.
Definition Logger.h:40
complex_t
std::complex< double > complex_t
Convenience alias for std::complex<double>.
Definition Utils.h:35
DecayParent::w_config
WilsonBuildConfig w_config
Wilson build configuration used when enabling this decay (scales, order, groups).
Definition DecayParent.h:117
DecayParent::bins
std::optional< std::vector< std::pair< double, double > > > bins
Optional q^2 bins.
Definition DecayParent.h:129
DecayParent::p
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).
Definition DecayParent.h:105
DecayParent::w_proxy
std::shared_ptr< IObsWilsonProxy > w_proxy
Wilson proxy used at compute-time to query coefficients (matching/run).
Definition DecayParent.h:99
GenericMapperWithExt< ObservableTag, Observables, LhaID, std::hash< LhaID >, observable_mapping, observable_flha_mapping >::enum_of
static std::optional< Observables > enum_of(const IdOf< ObservableTag > &id)
Attempts to recover the enum value associated with an identifier.
GenericMapperWithExt< ObservableTag, Observables, LhaID, std::hash< LhaID >, observable_mapping, observable_flha_mapping >::to_id
static IdOf< ObservableTag > to_id(Observables e)
Converts an enum value to an IdOf<Tag>.
GenericMapperWithExt< ObservableTag, Observables, LhaID, std::hash< LhaID >, observable_mapping, observable_flha_mapping >::str
static std::string str(const IdOf< ObservableTag > &id)
Returns the string representation of an identifier.
LbLFFCalculator::get
double get(LbL_FF a, double q2) override
Definition LbLFFCalculator.cpp:20
LbLllDecay::A_perp_0
complex_t A_perp_0(double q2, double sign, bool bar)
Definition LbLllDecay.cpp:114
LbLllDecay::A_FB_l
std::vector< ObservableValue > A_FB_l(Observables oid)
Definition LbLllDecay.cpp:489
LbLllDecay::A_par_0
complex_t A_par_0(double q2, double sign, bool bar)
Definition LbLllDecay.cpp:126
LbLllDecay::s_m
double s_m(double q2)
Definition LbLllDecay.cpp:81
LbLllDecay::A_FB_lh
std::vector< ObservableValue > A_FB_lh(Observables oid)
Definition LbLllDecay.cpp:510
LbLllDecay::A_perp_1
complex_t A_perp_1(double q2, double sign, bool bar)
Definition LbLllDecay.cpp:90
LbLllDecay::K2ss
double K2ss(double q2, bool bar)
Definition LbLllDecay.cpp:181
LbLllDecay::A_par_1
complex_t A_par_1(double q2, double sign, bool bar)
Definition LbLllDecay.cpp:102
LbLllDecay::compute_observable
std::vector< ObservableValue > compute_observable(Observables obs) override
Compute an observable given a public observable enum.
Definition LbLllDecay.cpp:541
LbLllDecay::beta_l
double beta_l(double q2)
Definition LbLllDecay.cpp:67
LbLllDecay::load_params
void load_params() override
Load and cache parameters needed by this decay.
Definition LbLllDecay.cpp:12
LbLllDecay::N
complex_t N(double q2, bool bar)
Definition LbLllDecay.cpp:85
LbLllDecay::dG_dq2_avg_bin
double dG_dq2_avg_bin(size_t bin)
Definition LbLllDecay.cpp:483
LbLllDecay::K2c
double K2c(double q2, bool bar)
Definition LbLllDecay.cpp:217
LbLllDecay::set_cfg_flags
void set_cfg_flags(LbLllConfig::Lepton gen)
Definition LbLllDecay.cpp:54
LbLllDecay::fill_wilson_cache
void fill_wilson_cache()
Definition LbLllDecay.cpp:39
LbLllDecay::set_n_threads
void set_n_threads(size_t n_threads) override
Set the number of worker threads used by decays that support parallel cache filling.
Definition LbLllDecay.cpp:616
LbLllDecay::load_cfg_dep_params
void load_cfg_dep_params()
Definition LbLllDecay.cpp:61
LbLllDecay::s_p
double s_p(double q2)
Definition LbLllDecay.cpp:77
LbLllDecay::F_T
std::vector< ObservableValue > F_T(Observables oid)
Definition LbLllDecay.cpp:531
LbLllDecay::K1c
double K1c(double q2, bool bar)
Definition LbLllDecay.cpp:172
LbLllDecay::A_FB_h
std::vector< ObservableValue > A_FB_h(Observables oid)
Definition LbLllDecay.cpp:499
LbLllDecay::K1cc
double K1cc(double q2, bool bar)
Definition LbLllDecay.cpp:155
LbLllDecay::F_L
std::vector< ObservableValue > F_L(Observables oid)
Definition LbLllDecay.cpp:520
LbLllDecay::compute_binned_K_i
void compute_binned_K_i()
Definition LbLllDecay.cpp:226
LbLllDecay::K2cc
double K2cc(double q2, bool bar)
Definition LbLllDecay.cpp:199
LbLllDecay::K1ss
double K1ss(double q2, bool bar)
Definition LbLllDecay.cpp:138
LbLllDecay::lambda
double lambda(double q2)
Definition LbLllDecay.cpp:71
LbLllDecay::dBR_dq2_binned
std::vector< ObservableValue > dBR_dq2_binned(Observables oid, bool br)
Definition LbLllDecay.cpp:458
PI
constexpr double PI
Definition constants.h:7
HBAR
constexpr double HBAR
Definition constants.h:23
RT2
constexpr double RT2
Definition constants.h:15
LbLllConfig::ff_src
LbL_FF_Src ff_src
Definition LbLllDecay.h:16
LbLllConfig::n_threads
size_t n_threads
Definition LbLllDecay.h:18
LbLllConfig::gen
Lepton gen
Definition LbLllDecay.h:17
LbLllDecayCache::N_0
complex_t N_0
Definition LbLllDecay.h:29
LbLllDecayCache::C
std::map< WCoef, complex_t > C
Definition LbLllDecay.h:23
LbLllDecayCache::ff_calculator
LbLFFCalculator ff_calculator
Definition LbLllDecay.h:24
LbLllDecayCache::bin_widths
std::vector< double > bin_widths
Definition LbLllDecay.h:35
LbLllDecayCache::K_i_binned
std::array< std::vector< double >, 6 > K_i_binned
Definition LbLllDecay.h:33
LbLllDecayCache::K_i_bar_binned
std::array< std::vector< double >, 6 > K_i_bar_binned
Definition LbLllDecay.h:34
MassConfig
Configuration for computing a particle mass at a given scale.
Definition Configs.h:242
ParamId
Composite identifier for a single parameter.
Definition ParamID.h:57
WilsonBuildConfig::order
QCDOrder order
Perturbative QCD order used for the evolution and matching of Wilson coefficients....
Definition Configs.h:54