Hyperiso 1.0.3
Modular flavour-physics calculations, Wilson coefficients and statistical inference
Loading...
Searching...
No Matches
BaseQCDfCalculator.cpp
Go to the documentation of this file.
1#include "BaseQCDfCalculator.h"
2
3double BaseQCDfCalculator::E(double q2) {
4 return (m_B * m_B + m_X * m_X - q2) / (2 * m_B);
5}
6
7BaseQCDfCalculator::BaseQCDfCalculator(int B_id, int X_id, double mu_b, const std::map<WCoef, complex_t> &C, B_FF_Type ff_tp,
8 std::shared_ptr<IObsParameterProxy<ParamId, DataType, std::string, LhaID>> p,
9 std::shared_ptr<IObsQCDProxy> iobs_qcdp) :
10 mu_b(mu_b), C(C), ff_tp(ff_tp), iobs_qcdp(iobs_qcdp)
11{
12 if (!this->allowed_decays.contains({B_id, X_id})) {
13 LOG_ERROR("ValueError", "Wrong meson PDG code in BaseQCDfCalculator constructor:", B_id, ",", X_id);
14 }
15
16 this->src_block = allowed_decays.at({B_id, X_id});
17 this->delta_qu = double (B_id == 521);
18 this->fill_wilson_bar_cache();
19
20 double beta_0 = iobs_qcdp->get_constants()->beta[4][0];
21 auto run = [this, beta_0] (double value_1gev, double eta, double gamma) { return value_1gev * pow(eta, gamma / beta_0); };
22
23 this->Lambda_h = (*p)(ParamId{ParameterType::DECAY, this->src_block, 14}, DataType::VALUE);
24 double mu_f = sqrt(this->mu_b * this->Lambda_h);
25 this->alpha_s_mu_b = (*iobs_qcdp)(AlphasConfig(this->mu_b, MassType::POLE, MassType::POLE));
26 this->alpha_s_mu_f = (*iobs_qcdp)(AlphasConfig(mu_f, MassType::POLE, MassType::POLE));
27 this->loop_f_mu_f = this->alpha_s_mu_f * iobs_qcdp->get_constants()->C_F / (4 * PI);
28 this->loop_f_mu_b = this->alpha_s_mu_b * iobs_qcdp->get_constants()->C_F / (4 * PI);
29 this->m_c_pole = (*p)(ParamId{ParameterType::SM, "QCD", {4, 2}}, DataType::VALUE);
30 this->m_b_pole = (*p)(ParamId{ParameterType::SM, "QCD", {5, 5}}, DataType::VALUE);
31 double eta_f = this->alpha_s_mu_f / (*iobs_qcdp)(AlphasConfig(1.0, MassType::POLE, MassType::POLE));
32 double eta_b = this->alpha_s_mu_b / (*iobs_qcdp)(AlphasConfig(1.0, MassType::POLE, MassType::POLE));
33 double m_b_pole_2loop = (*p)(ParamId{ParameterType::SM, "QCD", {5, 2}}, DataType::VALUE);
34 this->m_b_PS = m_b_pole_2loop - 4 * (*iobs_qcdp)(AlphasConfig(m_b_pole_2loop, MassType::POLE, MassType::POLE)) * mu_f / (3 * PI);
35 this->m_B = (*p)(ParamId{ParameterType::FLAVOR, "FMASS", B_id}, DataType::VALUE);
36 this->m_Bd = (*p)(ParamId{ParameterType::FLAVOR, "FMASS", 511}, DataType::VALUE);
37 this->m_X = (*p)(ParamId{ParameterType::FLAVOR, "FMASS", X_id}, DataType::VALUE);
38 this->f_B = (*p)(ParamId{ParameterType::FLAVOR, "FCONST", {B_id, 1}}, DataType::VALUE);
39 this->f_X_par = (*p)(ParamId{ParameterType::FLAVOR, "FCONST", {X_id, 1}}, DataType::VALUE);
40 this->lambda_B_p = (*p)(ParamId{ParameterType::DECAY, this->src_block, 13}, DataType::VALUE) / (1. - this->alpha_s_mu_f * log(pow(mu_f, 2)) * 1.8 / (3. * PI));
41 this->lambda_hat_u = std::conj((*p)(ParamId{ParameterType::SM, "VCKM", {0, 1}}, DataType::VALUE)) * (*p)(ParamId{ParameterType::SM, "VCKM", {0, 2}}, DataType::VALUE)
42 / (std::conj((*p)(ParamId{ParameterType::SM, "VCKM", {2, 1}}, DataType::VALUE)) * (*p)(ParamId{ParameterType::SM, "VCKM", {2, 2}}, DataType::VALUE));
43 // this->lambda_hat_u = 0.0; // ASK : Why neglected in SI for B > K l l ?
44 this->a_1_par = run((*p)(ParamId{ParameterType::DECAY, this->src_block, {8, 1}}, DataType::VALUE), eta_f, gamma_par(1));
45 this->a_2_par = run((*p)(ParamId{ParameterType::DECAY, this->src_block, {8, 2}}, DataType::VALUE), eta_f, gamma_par(2));
46 this->a_1_par_b = run((*p)(ParamId{ParameterType::DECAY, this->src_block, {8, 1}}, DataType::VALUE), eta_b, gamma_par(1));
47 this->a_2_par_b = run((*p)(ParamId{ParameterType::DECAY, this->src_block, {8, 2}}, DataType::VALUE), eta_b, gamma_par(2));
48
49
50 this->e_q = B_id == 521 ? e_u : e_d;
51 this->z_c = std::pow(this->m_c_pole / this->m_b_PS, 2);
52 this->z_c_IA = std::pow((*iobs_qcdp)(MassConfig(4, mu_b, MassType::MSBAR, MassType::POLE)) / (*iobs_qcdp)(MassConfig(5, mu_b, MassType::MSBAR, MassType::POLE)), 2);
53 this->L_b = std::log(this->mu_b / this->m_b_PS);
54 this->Delta_M = -6. * this->L_b - 4. * (1. - mu_f / this->m_b_PS);
55 int Nc = iobs_qcdp->get_constants()->Nc;
56 this->pref_par = PI2 * this->f_B * this->f_X_par / (Nc * this->m_B);
57
58
59 if (X_id == 333) {
60 this->n_T_par_m_0 = -4. * this->m_B / this->m_b_PS * (this->C_bar[WCoef::C3] + 3. * this->C_bar[WCoef::C4] + 12. * (this->C[WCoef::C3] + 10. * this->C[WCoef::C5]) - this->lambda_hat_u * (4. / 3. * this->C[WCoef::C1] + this->C[WCoef::C2]));
61 this->n_T_par_m_0_bar = -4. * this->m_B / this->m_b_PS * (this->C_bar[WCoef::C3] + 3. * this->C_bar[WCoef::C4] + 12. * (this->C[WCoef::C3] + 10. * this->C[WCoef::C5]) - std::conj(this->lambda_hat_u) * (4. / 3. * this->C[WCoef::C1] + this->C[WCoef::C2]));
62 } else {
63 this->n_T_par_m_0 = 4. * this->m_B / this->m_b_PS * (3. * this->lambda_hat_u * this->delta_qu * this->C[WCoef::C2] - this->C_bar[WCoef::C3] - 3. * this->C_bar[WCoef::C4]);
64 this->n_T_par_m_0_bar = 4. * this->m_B / this->m_b_PS * (3. * std::conj(this->lambda_hat_u) * this->delta_qu * this->C[WCoef::C2] - this->C_bar[WCoef::C3] - 3. * this->C_bar[WCoef::C4]);
65 }
66
67
68 bool isV = (X_id == 313 || X_id == 323 || X_id == 333);
69 if (isV) {
70 this->f_X_perp = run((*p)(ParamId{ParameterType::FLAVOR, "FCONST", {X_id, 2}}, DataType::VALUE), eta_f, iobs_qcdp->get_constants()->C_F);
71 this->f_X_perp_b = run((*p)(ParamId{ParameterType::FLAVOR, "FCONST", {X_id, 2}}, DataType::VALUE), eta_b, iobs_qcdp->get_constants()->C_F);
72 this->a_1_perp = run((*p)(ParamId{ParameterType::DECAY, this->src_block, {7, 1}}, DataType::VALUE), eta_f, gamma_perp(1));
73 this->a_2_perp = run((*p)(ParamId{ParameterType::DECAY, this->src_block, {7, 2}}, DataType::VALUE), eta_f, gamma_perp(2));
74 this->a_1_perp_b = run((*p)(ParamId{ParameterType::DECAY, this->src_block, {7, 1}}, DataType::VALUE), eta_b, gamma_perp(1));
75 this->a_2_perp_b = run((*p)(ParamId{ParameterType::DECAY, this->src_block, {7, 2}}, DataType::VALUE), eta_b, gamma_perp(2));
76
77
78 this->zeta_3_A = (*p)(ParamId{ParameterType::DECAY, this->src_block, 9}, DataType::VALUE);
79 this->zeta_3_V = (*p)(ParamId{ParameterType::DECAY, this->src_block, 10}, DataType::VALUE);
80 this->omega_10_A = (*p)(ParamId{ParameterType::DECAY, this->src_block, 11}, DataType::VALUE);
81 this->delta_t_p = (*p)(ParamId{ParameterType::DECAY, this->src_block, {12, 1}}, DataType::VALUE);
82 this->delta_t_m = (*p)(ParamId{ParameterType::DECAY, this->src_block, {12, 2}}, DataType::VALUE);
83 this->pref_perp = PI2 * this->f_B * this->f_X_perp / (Nc * this->m_B);
84 }
85
86}
87
88double BaseQCDfCalculator::phi_X(double u, double a1, double a2) {
89 double x = 2 * u - 1;
90 double C1 = 3 * x;
91 double C2 = -1.5 + 7.5 * x * x;
92
93 return 6 * u * (1 - u) * (1 + a1 * C1 + a2 * C2);
94}
95
96complex_t BaseQCDfCalculator::Y(double q2) {
97 return BV::h(q2, this->m_c_pole, this->mu_b) * (4./3. * this->C[WCoef::C1] + this->C[WCoef::C2] + 6. * this->C[WCoef::C3] + 60. * this->C[WCoef::C5])
98 - 0.5 * BV::h(q2, this->m_b_pole, this->mu_b) * (7. * this->C[WCoef::C3] + 4./3. * this->C[WCoef::C4] + 76. * this->C[WCoef::C5] + 64./3. * this->C[WCoef::C6])
99 - 0.5 * BV::h(q2, 0., this->mu_b) * (this->C[WCoef::C3] + 4./3. * this->C[WCoef::C4] + 16. * this->C[WCoef::C5] + 64./3. * this->C[WCoef::C6])
100 + 4./3. * this->C[WCoef::C3] + 64./9. * this->C[WCoef::C5] + 64./27. * this->C[WCoef::C6];
101}
102
103complex_t BaseQCDfCalculator::Y_u(double q2) {
104 return (BV::h(q2, this->m_c_pole, this->mu_b) - BV::h(q2, 0., this->mu_b)) * (4. / 3. * this->C[WCoef::C1] + this->C[WCoef::C2]);
105}
106
107complex_t BaseQCDfCalculator::t_perp(double u, double m_q, double q2, double E_Kstar) {
108 double mB2 = this->m_Bd * this->m_Bd;
109 // ASK : Why always m_Bd and not m_B for B+ decay ?
110 if(fpeq(q2, 0.)) {
111 if (fpeq(m_q, 0.)) return 4./(1.-u);
112 double epsilon=1.e-10;
113 complex_t xp=0.5+std::sqrt(0.25-(m_q*m_q-I*epsilon)/((1.-u)*mB2));
114 complex_t xm=0.5-std::sqrt(0.25-(m_q*m_q-I*epsilon)/((1.-u)*mB2));
115 return 4./(1.-u)*(1.+2.*m_q*m_q/(1.-u)/mB2*(BV::L_1(xp)+BV::L_1(xm)));
116 } else {
117 double s_hat = q2 / mB2;
118 double mq_hat = m_q / m_Bd;
119 return 2.*m_Bd/(1.-u)/E_Kstar*BV::I_1(u,s_hat,mq_hat)+q2/(1.-u)/(1.-u)/E_Kstar/E_Kstar*(BV::B_0((1.-u)*mB2+u*q2,m_q)-BV::B_0(q2,m_q));
120 }
121}
122
123complex_t BaseQCDfCalculator::t_par(double u, double m_q, double q2, double E_Kstar) {
124 double mB2 = this->m_Bd * this->m_Bd;
125 double s_hat = q2 / mB2;
126 double mq_hat = m_q / this->m_Bd;
127
128 return 2.*this->m_Bd/(1.-u)/E_Kstar*BV::I_1(u,s_hat,mq_hat)+((1.-u)*mB2+u*q2)/(1.-u)/(1.-u)/E_Kstar/E_Kstar*(BV::B_0((1.-u)*mB2+u*q2,m_q)-BV::B_0(q2,m_q));
129
130}
131
132complex_t BaseQCDfCalculator::T_par_m_0(bool bar) {
133 return bar ? this->n_T_par_m_0_bar : this->n_T_par_m_0;
134}
135
136void BaseQCDfCalculator::fill_wilson_bar_cache() {
137 auto b_ids = WCoefMapper::B_group();
138 for (size_t i = 0; i < 6; i++) {
139 this->C_bar[b_ids[i]] = 0;
140 for (size_t j = 0; j < 6; j++) {
141 this->C_bar[b_ids[i]] += P_bar[i][j] * this->C[b_ids[j]];
142 }
143 }
144}
145
146double BaseQCDfCalculator::gv_dga_4(double u) {
147 double a1 = -60. * this->zeta_3_A * (this->omega_10_A + 4.) + 1680. * this->zeta_3_V;
148 double a2 = 30. * this->zeta_3_A * (15. * this->omega_10_A + 32.) - 12600. * this->zeta_3_V + 36. * this->a_1_par - 72. * this->a_2_par - 12.;
149 double a3 = -100. * this->zeta_3_A * (9. * this->omega_10_A + 8.) + 25200. * this->zeta_3_V - 48. * this->a_1_par + 240. * this->a_2_par;
150 double a4 = 525. * this->zeta_3_A * this->omega_10_A - 14700. * this->zeta_3_V - 180. * this->a_2_par;
151 return -u * (a1 + u * (a2 + u * (a3 + u * a4))) / 4. + this->delta_t_p * (9. * u - 1.5) + this->delta_t_m * 6. * u + 3. * (this->delta_t_p + this->delta_t_m) * log(1 - u);
152}
153
154complex_t BaseQCDfCalculator::F_V(double v, bool bar) {
155 complex_t l_u = bar ? std::conj(this->lambda_hat_u) : this->lambda_hat_u;
156 return .75 * (
157 BV::h(v, this->m_c_pole, this->mu_b) * (this->C_bar[WCoef::C2] + this->C_bar[WCoef::C4] + this->C_bar[WCoef::C6] + l_u * (this->C[WCoef::C2] - this->C[WCoef::C1] / 6.))
158 + BV::h(v, this->m_b_pole, this->mu_b) * (this->C_bar[WCoef::C3] + this->C_bar[WCoef::C4] + this->C_bar[WCoef::C6])
159 + BV::h(v, 0., this->mu_b) * (this->C_bar[WCoef::C3] + 3. * this->C_bar[WCoef::C4] + 3. * this->C_bar[WCoef::C6] - l_u * (this->C[WCoef::C2] - this->C[WCoef::C1] / 6.))
160 - 8. / 27. * (this->C_bar[WCoef::C3] - this->C_bar[WCoef::C5] - 15. * this->C_bar[WCoef::C6])
161 );
162}
163
164double BaseQCDfCalculator::L(double q2) {
165 if (fpeq(q2, 0.0)) return 1.0;
166
167 double mb2 = std::pow(this->m_b_PS, 2);
168 return (q2 - mb2) * std::log(1. - q2 / mb2) / q2;
169}
170
171complex_t BaseQCDfCalculator::C_perp_0(double q2, double sign, bool bar) {
172 complex_t l_u = bar ? std::conj(this->lambda_hat_u) : this->lambda_hat_u;
173 complex_t C7 = this->C[WCoef::C7] + sign * this->C[WCoef::CP7];
174 if (bar) C7 = std::conj(C7);
175
176
177 if (ff_tp == B_FF_Type::FULL) {
178 if (fpeq(q2, 0.0)) return 0.0;
179 return q2 * l_u * Y_u(q2) / (2. * this->m_b_PS * this->m_B);
180 } else {
181 if (fpeq(q2, 0.0)) return C7;
182 return C7 + q2 * (Y(q2) + l_u * Y_u(q2)) / (2. * this->m_b_PS * this->m_B);
183 }
184}
185
186complex_t BaseQCDfCalculator::C_par_0(double q2, double sign, bool bar) {
187 complex_t l_u = bar ? std::conj(this->lambda_hat_u) : this->lambda_hat_u;
188 complex_t C7 = this->C[WCoef::C7] + sign * this->C[WCoef::CP7];
189 if (bar) C7 = std::conj(C7);
190
191 if (ff_tp == B_FF_Type::FULL) {
192 return -this->m_B * l_u * Y_u(q2) / (2. * this->m_b_PS);
193 } else {
194 return -C7 - this->m_B * (Y(q2) + l_u * Y_u(q2)) / (2. * this->m_b_PS);
195 }
196}
197
198complex_t BaseQCDfCalculator::C_perp_f(double q2, double sign, bool bar) {
199 if (ff_tp == B_FF_Type::FULL) return 0.0;
200
201 complex_t C7 = this->C[WCoef::C7] + sign * this->C[WCoef::CP7];
202 if (bar) C7 = std::conj(C7);
203
204 return C7 * (-2. * this->L_b - L(q2) + this->Delta_M);
205}
206
207complex_t BaseQCDfCalculator::C_par_f(double q2, double sign, bool bar) {
208 if (ff_tp == B_FF_Type::FULL) return 0.0;
209
210 complex_t C7 = this->C[WCoef::C7] + sign * this->C[WCoef::CP7];
211 if (bar) C7 = std::conj(C7);
212 return -C7 * (2. * std::log(this->m_b_PS / this->mu_b) + 2. * L(q2) + this->Delta_M);
213}
214
215complex_t BaseQCDfCalculator::C_perp_nf(double q2, bool bar) {
216 double s_hat = q2 / (this->m_b_PS * this->m_b_PS);
217 complex_t l_u = bar ? std::conj(this->lambda_hat_u) : this->lambda_hat_u;
218 complex_t F_27 = BV::f_27(s_hat, this->L_b, this->z_c) * (1. + l_u) + BV::f_27_u(s_hat, this->L_b) * l_u;
219
220 if (fpeq(q2, 0.0))
221 return -(this->C_bar[WCoef::C2] * F_27 + this->C[WCoef::C8] * BV::f_87(s_hat, this->L_b)) / iobs_qcdp->get_constants()->C_F;
222
223 complex_t F_19 = BV::f_19_PS(s_hat, this->L_b, this->z_c) * (1. + l_u) + BV::f_19_u(s_hat, this->L_b) * l_u;
224 complex_t F_29 = BV::f_29_PS(s_hat, this->L_b, this->z_c) * (1. + l_u) + BV::f_29_u(s_hat, this->L_b) * l_u;
225
226 return -(
227 this->C_bar[WCoef::C2] * F_27
228 + this->C[WCoef::C8] * BV::f_87(s_hat, this->L_b)
229 + q2 / (2. * this->m_b_PS * this->m_B) * (
230 (this->C_bar[WCoef::C2] + this->C_bar[WCoef::C1] / 3.) * F_29
231 + 2. * this->C_bar[WCoef::C1] * F_19
232 + this->C[WCoef::C8] * BV::f_89(s_hat)
233 )
234 ) / iobs_qcdp->get_constants()->C_F;
235}
236
237complex_t BaseQCDfCalculator::C_par_nf(double q2, bool bar) {
238 double s_hat = q2 / (this->m_b_PS * this->m_b_PS);
239 complex_t l_u = bar ? std::conj(this->lambda_hat_u) : this->lambda_hat_u;
240 complex_t F_27 = BV::f_27(s_hat, this->L_b, this->z_c) * (1. + l_u) + BV::f_27_u(s_hat, this->L_b) * l_u;
241 complex_t F_19 = BV::f_19_PS(s_hat, this->L_b, this->z_c) * (1. + l_u) + BV::f_19_u(s_hat, this->L_b) * l_u;
242 complex_t F_29 = BV::f_29_PS(s_hat, this->L_b, this->z_c) * (1. + l_u) + BV::f_29_u(s_hat, this->L_b) * l_u;
243
244 return (
245 this->C_bar[WCoef::C2] * F_27
246 + this->C[WCoef::C8] * BV::f_87(s_hat, this->L_b)
247 + this->m_B / (2 * this->m_b_PS) * (
248 (this->C_bar[WCoef::C2] + this->C_bar[WCoef::C1] / 3.) * F_29
249 + 2. * this->C_bar[WCoef::C1] * F_19
250 + this->C[WCoef::C8] * BV::f_89(s_hat))
251 ) / iobs_qcdp->get_constants()->C_F;
252}
253
254complex_t BaseQCDfCalculator::T_par_p_p_f(double u, double q2, bool bar) {
255 return 2. * T_perp_p_p_f(u, q2, bar);
256}
257
258complex_t BaseQCDfCalculator::T_par_p_m_f(double u, double q2, bool bar) {
259 return 2. * T_perp_p_m_f(u, q2, bar);
260}
261
262complex_t BaseQCDfCalculator::T_perp_p_p_f(double u, double q2, bool bar) {
263 complex_t C7 = this->C[WCoef::C7] + this->C[WCoef::CP7];
264 if (bar) C7 = std::conj(C7);
265 return 2. * this->m_B / (1. - u) / this->E(q2) * C7;
266}
267
268complex_t BaseQCDfCalculator::T_perp_p_m_f(double u, double q2, bool bar) {
269 complex_t C7 = this->C[WCoef::C7] - this->C[WCoef::CP7];
270 if (bar) C7 = std::conj(C7);
271 return 2. * this->m_B / (1. - u) / this->E(q2) * C7;
272}
273
274complex_t BaseQCDfCalculator::T_perp_p_nf(double u, double q2, bool bar) {
275 double E = this->E(q2);
276 complex_t l_u = bar ? std::conj(this->lambda_hat_u) : this->lambda_hat_u;
277 complex_t t_perp_mc = t_perp(u, this->m_c_pole, q2, E);
278 complex_t t_perp_mb = t_perp(u, this->m_b_PS, q2, E);
279 complex_t t_perp_0 = t_perp(u, 0.0, q2, E);
280
281 complex_t c8_term = fpeq(q2, 0.0) ? 0.0 : -4.0 * e_d * this->C[WCoef::C8] / (u + (1 - u) * q2 / (this->m_B * this->m_B));
282 return c8_term + this->m_B / (2 * this->m_b_PS) * (
283 e_u * (
284 t_perp_mc * (this->C_bar[WCoef::C2] + this->C_bar[WCoef::C4] - this->C_bar[WCoef::C6] + l_u * (this->C[WCoef::C2] - this->C[WCoef::C1] / 6.))
285 - t_perp_0 * l_u * (this->C[WCoef::C2] - this->C[WCoef::C1] / 6.)
286 )
287 + e_d * (
288 t_perp_mb * (this->C_bar[WCoef::C3] + this->C_bar[WCoef::C4] - this->C_bar[WCoef::C6] - 4. * this->m_b_PS / this->m_B * this->C_bar[WCoef::C5])
289 + t_perp_0 * this->C_bar[WCoef::C3]
290 )
291 );
292}
293
294complex_t BaseQCDfCalculator::T_par_p_nf(double u, double q2, bool bar) {
295 double E = this->E(q2);
296 complex_t l_u = bar ? std::conj(this->lambda_hat_u) : this->lambda_hat_u;
297 complex_t t_par_mc = t_par(u, this->m_c_pole, q2, E);
298 complex_t t_par_mb = t_par(u, this->m_b_PS, q2, E);
299 complex_t t_par_0 = t_par(u, 0., q2, E);
300
301 return this->m_B / this->m_b_PS * (
302 e_u * (
303 t_par_mc * (this->C_bar[WCoef::C2] + this->C_bar[WCoef::C4] - this->C_bar[WCoef::C6] + l_u * (this->C[WCoef::C2] - this->C[WCoef::C1] / 6.))
304 - t_par_0 * l_u * (this->C[WCoef::C2] - this->C[WCoef::C1] / 6.)
305 )
306 + e_d * (
307 t_par_mb * (this->C_bar[WCoef::C3] + this->C_bar[WCoef::C4] - this->C_bar[WCoef::C6])
308 + t_par_0 * this->C_bar[WCoef::C3]
309 )
310 );
311}
312
313complex_t BaseQCDfCalculator::T_par_m_nf(double u, double q2, bool bar) {
314 double v = this->m_B * this->m_B * (1 - u) + q2 * u;
315
316 return 8. * this->m_B * this->m_B * this->C[WCoef::C8] / v + 8. * this->m_B / this->m_b_PS * F_V(v, bar);
317}
318
319complex_t BaseQCDfCalculator::inv_lambda_B_m(double q2) {
320 double omega_0 = 2. * (this->m_B - this->m_b_PS) / 3.;
321 double x = q2 / (this->m_B * omega_0);
322 return std::exp(-x) / omega_0 * (I * PI - Ei(x));
323}
324
325complex_t BaseQCDfCalculator::I_perp_p(double q2, bool bar) {
326 double pref = this->loop_f_mu_f / this->lambda_B_p;
327
328 if (this->ff_tp == B_FF_Type::SOFT) {
329 auto f_soft = [q2, bar, this] (double u) {
330 return phi_X(u, this->a_1_perp, this->a_2_perp) * (T_perp_p_p_f(u, q2, bar) + T_perp_p_nf(u, q2, bar));
331 };
332 return pref * c_integrate(f_soft, 0, 1, 1e-2);
333 } else {
334 auto f_full = [q2, bar, this] (double u) {
335 return phi_X(u, this->a_1_perp, this->a_2_perp) * T_perp_p_nf(u, q2, bar);
336 };
337 return pref * c_integrate(f_full, 0, 1, 1e-2);
338 }
339}
340
341complex_t BaseQCDfCalculator::I_perp_m(double q2, bool bar) {
342 double pref = this->loop_f_mu_f / this->lambda_B_p;
343
344 if (this->ff_tp == B_FF_Type::SOFT) {
345 auto f_soft = [q2, bar, this] (double u) {
346 return phi_X(u, this->a_1_perp, this->a_2_perp) * (T_perp_p_m_f(u, q2, bar) + T_perp_p_nf(u, q2, bar));
347 };
348 return pref * c_integrate(f_soft, 0, 1, 1e-2);
349 } else {
350 auto f_full = [q2, bar, this] (double u) {
351 return phi_X(u, this->a_1_perp, this->a_2_perp) * T_perp_p_nf(u, q2, bar);
352 };
353 return pref * c_integrate(f_full, 0, 1, 1e-2);
354 }
355}
356
357complex_t BaseQCDfCalculator::I_par_p(double q2, bool bar) {
358 if (this->ff_tp == B_FF_Type::SOFT) {
359 auto f_soft = [q2, bar, this] (double u) {
360 double phi = phi_X(u, this->a_1_par, this->a_2_par);
361 complex_t i1 = phi * (T_par_p_p_f(u, q2, bar) + T_par_p_nf(u, q2, bar));
362 complex_t i2 = phi * (T_par_m_0(bar) + this->loop_f_mu_f * T_par_m_nf(u, q2, bar));
363 return this->loop_f_mu_f / this->lambda_B_p * i1 + this->e_q * inv_lambda_B_m(q2) * i2;
364 };
365 return c_integrate(f_soft, 0, 1, 1e-2);
366 } else {
367 auto f_full = [q2, bar, this] (double u) {
368 double phi = phi_X(u, this->a_1_par, this->a_2_par);
369 complex_t i1 = phi * T_par_p_nf(u, q2, bar);
370 complex_t i2 = phi * (T_par_m_0(bar) + this->loop_f_mu_f * T_par_m_nf(u, q2, bar));
371 return this->loop_f_mu_f / this->lambda_B_p * i1 + this->e_q * inv_lambda_B_m(q2) * i2;
372 };
373 return c_integrate(f_full, 0, 1, 1e-2);
374 }
375}
376
377complex_t BaseQCDfCalculator::I_par_m(double q2, bool bar) {
378 if (this->ff_tp == B_FF_Type::SOFT) {
379 auto f_soft = [q2, bar, this] (double u) {
380 double phi = phi_X(u, this->a_1_par, this->a_2_par);
381 complex_t i1 = phi * (T_par_p_m_f(u, q2, bar) + T_par_p_nf(u, q2, bar));
382 complex_t i2 = phi * (T_par_m_0(bar) + this->loop_f_mu_f * T_par_m_nf(u, q2, bar));
383 return this->loop_f_mu_f / this->lambda_B_p * i1 + this->e_q * inv_lambda_B_m(q2) * i2;
384 };
385 return c_integrate(f_soft, 0, 1, 1e-2);
386 } else {
387 auto f_full = [q2, bar, this] (double u) {
388 double phi = phi_X(u, this->a_1_par, this->a_2_par);
389 complex_t i1 = phi * T_par_p_nf(u, q2, bar);
390 complex_t i2 = phi * (T_par_m_0(bar) + this->loop_f_mu_f * T_par_m_nf(u, q2, bar));
391 return this->loop_f_mu_f / this->lambda_B_p * i1 + this->e_q * inv_lambda_B_m(q2) * i2;
392 };
393 return c_integrate(f_full, 0, 1, 1e-2);
394 }
395}
LOG_ERROR
#define LOG_ERROR(type,...)
Macro for logging error messages and terminating the application.
Definition Logger.h:41
complex_t
std::complex< double > complex_t
Convenience alias for std::complex<double>.
Definition Utils.h:35
BaseQCDfCalculator::t_perp
complex_t t_perp(double u, double m_q, double q2, double E_Kstar)
Definition BaseQCDfCalculator.cpp:107
BaseQCDfCalculator::inv_lambda_B_m
complex_t inv_lambda_B_m(double q2)
Definition BaseQCDfCalculator.cpp:319
BaseQCDfCalculator::I_perp_m
complex_t I_perp_m(double q2, bool bar)
Definition BaseQCDfCalculator.cpp:341
BaseQCDfCalculator::BaseQCDfCalculator
BaseQCDfCalculator()=default
BaseQCDfCalculator::C_par_0
complex_t C_par_0(double q2, double sign, bool bar)
Definition BaseQCDfCalculator.cpp:186
BaseQCDfCalculator::phi_X
double phi_X(double u, double a1, double a2)
Definition BaseQCDfCalculator.cpp:88
BaseQCDfCalculator::C_par_nf
complex_t C_par_nf(double q2, bool bar)
Definition BaseQCDfCalculator.cpp:237
BaseQCDfCalculator::T_par_m_nf
complex_t T_par_m_nf(double u, double q2, bool bar)
Definition BaseQCDfCalculator.cpp:313
BaseQCDfCalculator::e_u
static constexpr double e_u
Definition BaseQCDfCalculator.h:37
BaseQCDfCalculator::T_par_p_p_f
complex_t T_par_p_p_f(double u, double q2, bool bar)
Definition BaseQCDfCalculator.cpp:254
BaseQCDfCalculator::iobs_qcdp
std::shared_ptr< IObsQCDProxy > iobs_qcdp
Definition BaseQCDfCalculator.h:35
BaseQCDfCalculator::E
double E(double q2)
Definition BaseQCDfCalculator.cpp:3
BaseQCDfCalculator::e_d
static constexpr double e_d
Definition BaseQCDfCalculator.h:38
BaseQCDfCalculator::T_perp_p_p_f
complex_t T_perp_p_p_f(double u, double q2, bool bar)
Definition BaseQCDfCalculator.cpp:262
BaseQCDfCalculator::Y
complex_t Y(double q2)
Definition BaseQCDfCalculator.cpp:96
BaseQCDfCalculator::C_bar
std::map< WCoef, complex_t > C_bar
Definition BaseQCDfCalculator.h:29
BaseQCDfCalculator::gv_dga_4
double gv_dga_4(double u)
Definition BaseQCDfCalculator.cpp:146
BaseQCDfCalculator::P_bar
static constexpr std::array< std::array< double, 6 >, 6 > P_bar
Definition BaseQCDfCalculator.h:39
BaseQCDfCalculator::C_perp_0
complex_t C_perp_0(double q2, double sign, bool bar)
Definition BaseQCDfCalculator.cpp:171
BaseQCDfCalculator::T_par_m_0
complex_t T_par_m_0(bool bar)
Definition BaseQCDfCalculator.cpp:132
BaseQCDfCalculator::I_par_p
complex_t I_par_p(double q2, bool bar)
Definition BaseQCDfCalculator.cpp:357
BaseQCDfCalculator::F_V
complex_t F_V(double v, bool bar)
Definition BaseQCDfCalculator.cpp:154
BaseQCDfCalculator::C_perp_nf
complex_t C_perp_nf(double q2, bool bar)
Definition BaseQCDfCalculator.cpp:215
BaseQCDfCalculator::T_perp_p_nf
complex_t T_perp_p_nf(double u, double q2, bool bar)
Definition BaseQCDfCalculator.cpp:274
BaseQCDfCalculator::t_par
complex_t t_par(double u, double m_q, double q2, double E_Kstar)
Definition BaseQCDfCalculator.cpp:123
BaseQCDfCalculator::T_par_p_nf
complex_t T_par_p_nf(double u, double q2, bool bar)
Definition BaseQCDfCalculator.cpp:294
BaseQCDfCalculator::C_par_f
complex_t C_par_f(double q2, double sign, bool bar)
Definition BaseQCDfCalculator.cpp:207
BaseQCDfCalculator::T_perp_p_m_f
complex_t T_perp_p_m_f(double u, double q2, bool bar)
Definition BaseQCDfCalculator.cpp:268
BaseQCDfCalculator::I_par_m
complex_t I_par_m(double q2, bool bar)
Definition BaseQCDfCalculator.cpp:377
BaseQCDfCalculator::C
std::map< WCoef, complex_t > C
Definition BaseQCDfCalculator.h:28
BaseQCDfCalculator::Y_u
complex_t Y_u(double q2)
Definition BaseQCDfCalculator.cpp:103
BaseQCDfCalculator::I_perp_p
complex_t I_perp_p(double q2, bool bar)
Definition BaseQCDfCalculator.cpp:325
BaseQCDfCalculator::T_par_p_m_f
complex_t T_par_p_m_f(double u, double q2, bool bar)
Definition BaseQCDfCalculator.cpp:258
BaseQCDfCalculator::C_perp_f
complex_t C_perp_f(double q2, double sign, bool bar)
Definition BaseQCDfCalculator.cpp:198
BaseQCDfCalculator::allowed_decays
static const std::map< LhaID, std::string > allowed_decays
Definition BaseQCDfCalculator.h:48
C1
Definition BWilson.h:6
C2
Definition BWilson.h:20
C7
Definition BWilson.h:99
WCoefMapper::B_group
static const std::vector< WCoef > & B_group()
B-group (C1–C10) Wilson coefficients.
Definition wcoef_ids.hpp:239
PI
constexpr double PI
Definition constants.h:7
I
constexpr std::complex< double > I
Definition constants.h:20
PI2
constexpr double PI2
Definition constants.h:8
E
constexpr double E
Definition constants.h:13
Ei
double Ei(double x)
Computes the exponential integral function Ei(x).
Definition exponential_integral.cpp:106
epsilon
constexpr double epsilon
Definition exponential_integral.cpp:3
c_integrate
scalar_t c_integrate(ComplexValuedFunction f, double l, double u, double prec)
Performs numerical integration of a complex-valued function of a real variable.
Definition integration.cpp:63
BV::h
complex_t h(double s, double m_q, double mu_b)
Definition ObsMath.cpp:11530
BV::I_1
complex_t I_1(double u, double s_hat, double m_q_hat)
Definition ObsMath.cpp:11552
BV::f_19_u
complex_t f_19_u(double s_hat, double L_b)
Definition ObsMath.cpp:11456
BV::f_27_u
complex_t f_27_u(double s_hat, double L_b)
Definition ObsMath.cpp:11429
BV::f_19_PS
complex_t f_19_PS(double s_hat, double L_b, double z, size_t max_pow=20)
Definition ObsMath.cpp:11443
BV::f_27
complex_t f_27(double s_hat, double L_b, double z, size_t max_pow=20)
Definition ObsMath.cpp:11406
BV::B_0
complex_t B_0(double s, double m_q)
Definition ObsMath.cpp:11543
BV::f_29_u
complex_t f_29_u(double s_hat, double L_b)
Definition ObsMath.cpp:11485
BV::f_87
complex_t f_87(double s_hat, double L_b)
Definition ObsMath.cpp:11508
BV::L_1
complex_t L_1(complex_t x)
Definition ObsMath.cpp:11548
BV::f_29_PS
complex_t f_29_PS(double s_hat, double L_b, double z, size_t max_pow=20)
Definition ObsMath.cpp:11447
BV::f_89
complex_t f_89(double s_hat)
Definition ObsMath.cpp:11521
pow
scalar_t pow(const scalar_t &base, const scalar_t &exp)
Definition scalar.cpp:75
fpeq
std::enable_if_t< not std::numeric_limits< T >::is_integer, bool > fpeq(T, T, std::size_t n=10)
Compares two floating point numbers with a given precision.
Definition special_generic.h:107
AlphasConfig
Configuration for evaluating the strong coupling constant .
Definition Configs.h:206
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