Hyperiso 1.0.3
Modular flavour-physics calculations, Wilson coefficients and statistical inference
Loading...
Searching...
No Matches
QCDHelper.cpp
Go to the documentation of this file.
1#include "QCDHelper.h"
2
3void QCDHelper::Init() {
4 LOG_DEBUG("Init QCD dependent block");
5 std::unordered_map<ParameterType, std::vector<std::string>> src = {{ParameterType::SM, {"SMINPUTS", "MASS"}}};
6
7 auto func = [] (const BlockSrc& src, std::shared_ptr<DependentBlock> dep_block) {
8 double m_Z = src.get_val("SMINPUTS" , 4);
9 double alpha_s_mZ = src.get_val("SMINPUTS" , 3);
10 double m_b_mb = src.get_val("SMINPUTS" , 5);
11 double m_t_pole = src.get_val("SMINPUTS" , 6);
12 double m_c = src.get_val("MASS" , 4);
13 // double m_s = src.get_val("MASS" , 3);
14 // double m_d = src.get_val("MASS" , 2);
15 // double m_u = src.get_val("MASS" , 1);
16
17 double lambda_5 = match_lambda(alpha_s_mZ, m_Z, 5);
18 double lambda_6_mt_pole = match_lambda(alpha_s_explicit(m_t_pole, lambda_5, 5), m_t_pole, 6);
19 double lambda_4_mb_mb = match_lambda(alpha_s_explicit(m_b_mb, lambda_5, 5), m_b_mb, 4);
20 double lambda_3_mb_mb = match_lambda(alpha_s_explicit(m_c, lambda_4_mb_mb, 4), m_c, 3);
21 double m_t_mt = calc_mt_mt(lambda_6_mt_pole, lambda_5);
22 double lambda_6_mt_mt = match_lambda(alpha_s_explicit(m_t_mt, lambda_5, 5), m_t_mt, 6);
23 double m_b_pole = calc_mb_pole(lambda_5);
24 double m_b_pole_one_loop = calc_mb_pole_one_loop(lambda_5);
25 double lambda_4_mb_pole = match_lambda(alpha_s_explicit(m_b_pole, lambda_5, 5), m_b_pole, 4);
26 double lambda_3_mb_pole = match_lambda(alpha_s_explicit(m_c, lambda_4_mb_pole, 4), m_c, 3);
27 double m_b_1S = calc_mb_1S(lambda_4_mb_pole, m_b_pole);
28 double m_c_pole = calc_mc_pole(lambda_4_mb_mb);
29 double m_c_pole_one_loop = calc_mc_pole_one_loop(lambda_4_mb_mb);
30 double m_b_kin = calc_mb_kinematic(m_b_mb, lambda_3_mb_mb);
31
32 dep_block->store_or_assign(LhaID(1, 3, 1), std::make_shared<Parameter>(ParamId{ParameterType::SM, "QCD", LhaID(1, 3, 1)}, lambda_3_mb_mb, 0., 0.));
33 dep_block->store_or_assign(LhaID(1, 3, 2), std::make_shared<Parameter>(ParamId{ParameterType::SM, "QCD", LhaID(1, 3, 2)}, lambda_3_mb_pole, 0., 0.));
34 dep_block->store_or_assign(LhaID(1, 4, 1), std::make_shared<Parameter>(ParamId{ParameterType::SM, "QCD", LhaID(1, 4, 1)}, lambda_4_mb_mb, 0., 0.));
35 dep_block->store_or_assign(LhaID(1, 4, 2), std::make_shared<Parameter>(ParamId{ParameterType::SM, "QCD", LhaID(1, 4, 2)}, lambda_4_mb_pole, 0., 0.));
36 dep_block->store_or_assign(LhaID(1, 5), std::make_shared<Parameter>(ParamId{ParameterType::SM, "QCD", LhaID(1, 5)}, lambda_5, 0., 0.));
37 dep_block->store_or_assign(LhaID(1, 6, 1), std::make_shared<Parameter>(ParamId{ParameterType::SM, "QCD", LhaID(1, 6, 1)}, lambda_6_mt_mt, 0., 0.));
38 dep_block->store_or_assign(LhaID(1, 6, 2), std::make_shared<Parameter>(ParamId{ParameterType::SM, "QCD", LhaID(1, 6, 2)}, lambda_6_mt_pole, 0., 0.));
39 dep_block->store_or_assign(LhaID(4, 1), std::make_shared<Parameter>(ParamId{ParameterType::SM, "QCD", LhaID(4, 1)}, m_c_pole, 0., 0.));
40 dep_block->store_or_assign(LhaID(4, 2), std::make_shared<Parameter>(ParamId{ParameterType::SM, "QCD", LhaID(4, 2)}, m_c_pole_one_loop, 0., 0.));
41 dep_block->store_or_assign(LhaID(5, 1), std::make_shared<Parameter>(ParamId{ParameterType::SM, "QCD", LhaID(5, 1)}, m_b_mb, 0., 0.));
42 dep_block->store_or_assign(LhaID(5, 2), std::make_shared<Parameter>(ParamId{ParameterType::SM, "QCD", LhaID(5, 2)}, m_b_pole, 0., 0.));
43 dep_block->store_or_assign(LhaID(5, 3), std::make_shared<Parameter>(ParamId{ParameterType::SM, "QCD", LhaID(5, 3)}, m_b_1S, 0., 0.));
44 dep_block->store_or_assign(LhaID(5, 4), std::make_shared<Parameter>(ParamId{ParameterType::SM, "QCD", LhaID(5, 4)}, m_b_kin, 0., 0.));
45 dep_block->store_or_assign(LhaID(5, 5), std::make_shared<Parameter>(ParamId{ParameterType::SM, "QCD", LhaID(5, 5)}, m_b_pole_one_loop, 0., 0.));
46 dep_block->store_or_assign(6, std::make_shared<Parameter>(ParamId{ParameterType::SM, "QCD", 6}, m_t_mt, 0., 0.));
47 };
48
49 DependentBlockManager::addDependentBlock("QCD", src, ParameterType::SM, func);
50}
51
52double QCDHelper::alpha_s(double mu, MassType mass_b_type, MassType mass_t_type) {
53 if (mu < (*Parameters::GetInstance())("MASS", 3)) {
54 std::ostringstream oss;
55 oss << "Renormalisation scale for alpha_s calculation (" << mu
56 << ") is below strange mass (" << (*Parameters::GetInstance())("MASS", 3) << ").";
57 throw std::domain_error(oss.str());
58 // LOG_ERROR("Scale Error", "Renormalisation scale for alpha_s calculation (", mu, ") is below strange mass (" , (*Parameters::GetInstance())("MASS", 3), ").");
59 }
60
61 return alpha_s_explicit(mu, get_lambda(mu, mass_b_type, mass_t_type), get_nf(mu, mass_b_type, mass_t_type));
62}
63
64double QCDHelper::msbar_mass(int pdg_code, double mu, MassType mass_b_type, MassType mass_t_type) {
65 if (pdg_code > 6 || pdg_code < 1) {
66 LOG_ERROR("ValueError", "PDG code", pdg_code, "is not a quark");
67 }
68
69 auto p = Parameters::GetInstance();
70 double quark_mass = pdg_code < 5 ? (*p)("MASS", pdg_code) : pdg_code == 5 ? (*p)("QCD", LhaID(5, 1)) : (*p)("QCD", 6);
71 double Qinit = pdg_code < 4 ? 1 : quark_mass;
72
73 int n_i = get_nf(Qinit, mass_b_type, mass_t_type);
74 int n_f = get_nf(mu, mass_b_type, mass_t_type);
75 auto Q_bounds = getOrderedMasses(mass_b_type, mass_t_type);
76
77 if (fpeq(Qinit, mu))
78 return quark_mass;
79
80 while (n_i > n_f) {
81 quark_mass = runMass(quark_mass, Qinit, Q_bounds.at(n_i - 1), n_i, mass_b_type, mass_t_type);
82 Qinit = Q_bounds.at(n_i - 1);
83 --n_i;
84 }
85
86 while (n_i < n_f) {
87 quark_mass = runMass(quark_mass, Qinit, Q_bounds.at(n_i), n_i, mass_b_type, mass_t_type);
88 Qinit = Q_bounds.at(n_i);
89 ++n_i;
90 }
91
92 return runMass(quark_mass, Qinit, mu, n_f, mass_b_type, mass_t_type);
93}
94
95double QCDHelper::calc_mc_pole(double lambda_4) {
96 double mc = (*Parameters::GetInstance())("MASS", 4);
97 double mu = (*Parameters::GetInstance())("MASS", 2);
98 double md = (*Parameters::GetInstance())("MASS", 1);
99 double ms = (*Parameters::GetInstance())("MASS", 3);
100 double alphas_mc = alpha_s_explicit(mc, lambda_4, 4);
101 return mc * (1 + alphas_mc / PI * (constants->C_F + alphas_mc / PI * ((13.4434 - 1.0414 * 3
102 + 1.0414 * constants->C_F * ((mu + md + ms) / mc)))));
103}
104
105double QCDHelper::calc_mc_pole_one_loop(double lambda_4) {
106 double mc = (*Parameters::GetInstance())("MASS", 4);
107 double alphas_mc = alpha_s_explicit(mc, lambda_4, 4);
108 return mc * (1 + alphas_mc / PI * constants->C_F);
109}
110
111double QCDHelper::calc_mb_pole(double lambda_5) {
112 double mc = (*Parameters::GetInstance())("MASS", 4);
113 double mu = (*Parameters::GetInstance())("MASS", 2);
114 double md = (*Parameters::GetInstance())("MASS", 1);
115 double ms = (*Parameters::GetInstance())("MASS", 3);
116 double mb = (*Parameters::GetInstance())("SMINPUTS", 5);
117 double alphas_mb = alpha_s_explicit(mb, lambda_5, 5);
118 return mb * (1. + alphas_mb / PI * (constants->C_F
119 + alphas_mb / PI * ((13.4434 - 1.0414 * 4. + 1.0414 * constants->C_F * (mu + md + ms + mc) / mb))));
120}
121
122double QCDHelper::calc_mb_pole_one_loop(double lambda_5) {
123 double mb = (*Parameters::GetInstance())("SMINPUTS", 5);
124 double alphas_mb = alpha_s_explicit(mb, lambda_5, 5);
125 return mb * (1 + alphas_mb / PI * constants->C_F);
126}
127
128// From 2005.06487
129double QCDHelper::calc_mb_kinematic(double mb_mb, double lambda_3) {
130 int nl = 3;
131 double mu = 1.0;
132 double mu_s = mb_mb;
133 double lm = std::log(std::pow(mu_s / mb_mb, 2));
134 double lm2 = lm * lm;
135 double alpha_s = alpha_s_explicit(mu_s, lambda_3, nl);
136 double r = mu / mb_mb;
137 double r2 = r * r;
138 double lmu = std::log(2 * mu / mu_s);
139 double l2 = std::log(2.);
140
141 double eps_one_loop = 4. / 3 + lm - 16. / 9 * r - 2. / 3 * r2;
142 double eps_two_loop = 307. / 32 + PI2 / 3 - ZETA3 / 6 + PI2 / 9 * l2 + 509. / 72 * lm + 47 / 24. * lm2
143 - nl * (71. / 144 + PI2 / 18 + 13. / 36 * lm + 1. / 12 * lm2)
144 + r * (-860. / 27 + 8 * PI2 / 9 + 88. / 9 * lmu + nl * (128. / 81 - 16. / 27 * lmu))
145 + r2 * (-83. / 9 + PI2 / 3 + 2. / 3 * lm + 11. / 3 * lmu + nl * (13. / 27 - 2. / 9 * lmu));
146
147 return mb_mb * (1 + alpha_s / PI * eps_one_loop + std::pow(alpha_s / PI, 2) * eps_two_loop);
148}
149
150double QCDHelper::calc_mb_1S(double lambda_4, double mb_pole) {
151 double mu = (*Parameters::GetInstance())("SMINPUTS", 5) / 2.;
152 return mb_pole * (1 - 2. / 9 * pow(alpha_s_explicit(mu, lambda_4, 4), 2.));
153}
154
155// From 0906.5273
156double QCDHelper::calc_mt_mt(double lambda6_mt_pole, double lambda_5) {
157 double mt_pole = (*Parameters::GetInstance())("SMINPUTS", 6);
158 double alpha = alpha_s_explicit(mt_pole, lambda6_mt_pole, 6);
159 double a = 307. / 32. + PI2 / 3. + PI2 / 9. * log(2.) - 1. / 6 * ZETA3 - 71. / 144. * 5.;
160 double mt_mt = mt_pole / (1. + alpha / PI * (4. / 3 + alpha / PI * a));
161 double lambda = match_lambda(alpha_s_explicit(mt_mt, lambda_5, 5), mt_mt, 6);
162 alpha = alpha_s_explicit(mt_mt, lambda, 6);
163 return mt_pole / (1. + alpha / PI * (4. / 3 + alpha / PI * a));
164}
165
166int QCDHelper::get_nf(double Q, MassType mass_b_type, MassType mass_t_type) {
167 auto masses = getOrderedMasses(mass_b_type, mass_t_type);
168 // for (size_t i = 0; i < masses.size(); ++i) {
169 // if (1 - Q / std::abs(masses.at(i)) > 1e-4) {
170 // return i;
171 // }
172 // }
173
174 for (size_t i = 0; i < masses.size(); ++i) {
175 if (Q < masses.at(i))
176 return i;
177 }
178
179 return 6;
180}
181
182double QCDHelper::get_lambda(double mu, MassType mass_b_type, MassType mass_t_type) {
183 int nf = get_nf(mu, mass_b_type, mass_t_type);
184 if (nf == 3)
185 return (*Parameters::GetInstance())("QCD", LhaID(1, 3, mass_b_type == MassType::POLE ? 2 : 1));
186 if (nf == 4)
187 return (*Parameters::GetInstance())("QCD", LhaID(1, 4, mass_b_type == MassType::POLE ? 2 : 1));
188 if (nf == 6)
189 return (*Parameters::GetInstance())("QCD", LhaID(1, 6, mass_t_type == MassType::POLE ? 2 : 1));
190 return (*Parameters::GetInstance())("QCD", LhaID(1, nf));
191}
192
193std::vector<double> QCDHelper::getOrderedMasses(MassType mass_b_type, MassType mass_t_type) {
194 double m_b = (*Parameters::GetInstance())("QCD", LhaID(5, mass_b_type == MassType::POLE ? 2 : 1));
195 double m_t = mass_t_type == MassType::MSBAR ? (*Parameters::GetInstance())("QCD", 6) : (*Parameters::GetInstance())("SMINPUTS", 6);
196 double mc = (*Parameters::GetInstance())("MASS", 4);
197 double mu = (*Parameters::GetInstance())("MASS", 2);
198 double md = (*Parameters::GetInstance())("MASS", 1);
199 double ms = (*Parameters::GetInstance())("MASS", 3);
200 return {md, mu, ms, mc, m_b, m_t};
201}
202
203double QCDHelper::match_lambda(double target_alpha, double Q, int nf) {
204 if (!std::isfinite(target_alpha) || target_alpha <= 0.0 || target_alpha > 1.0) {
205 throw std::domain_error(
206 "match_lambda: unphysical target_alpha=" + std::to_string(target_alpha)
207 );
208 }
209
210 double L_min = 1e-3;
211 double L_max = 1.;
212 double L_moy = L_min;
213 double alphas_min=0;
214 double alphas_moy = alphas_min;
215
216 while ((std::abs(1.- L_min / L_max) > 1e-5)&& (std::abs(1. - alphas_min / target_alpha) > 1e-4) ){
217 alphas_min = alpha_s_explicit(Q, L_min, nf);
218 L_moy = (L_min + L_max) / 2;
219 alphas_moy = alpha_s_explicit(Q, L_moy, nf);
220
221 (target_alpha >= alphas_min && target_alpha <= alphas_moy) ? L_max = L_moy : L_min = L_moy;
222 }
223
224 if (std::abs(1-L_min/L_max) <= 1e-5) {
225 throw std::domain_error(
226 "match_lambda: target alpha outside bracket, alpha=" + std::to_string(target_alpha)
227 );
228 // LOG_ERROR("ValueError", "Unable to find suitable QCD Lambda value to match alpha_s = " + std::to_string(target_alpha)
229 // + " at scale " + std::to_string(Q) + " GeV with " + std::to_string(nf) + " active flavors.");
230 return -1;
231 }
232
233 return L_min;
234}
235
236double QCDHelper::alpha_s_explicit(double mu, double lambda, int nf) {
237 double r = std::pow(mu / lambda, 2);
238 double L = std::log(r);
239 double LL = std::log(L);
240 auto [b0, b1, b2] = constants->beta[nf - 1];
241 double b02 = b0 * b0;
242 double b12 = b1 * b1;
243 return 4 * PI * (1 - 2 * b1 * LL / (b02 * L) + 4 * b12 * (std::pow(LL - .5, 2) + b2 * b0 / 8 / b12 - 1.25) / std::pow(b02 * L, 2)) / (b0 * L);
244}
245
246double QCDHelper::runMass(double mass, double Q_i, double Q_f, int nf, MassType m_b_type, MassType m_t_type) {
247 if (fpeq(Q_i, Q_f)) return mass;
248
249 return mass * R(alpha_s(Q_f, m_b_type, m_t_type), nf)
250 / R(alpha_s(Q_i, m_b_type, m_t_type), nf);
251}
252
253double QCDHelper::R(double alpha, int nf) {
254 auto [b0, b1, b2] = constants->beta[nf - 1];
255 auto [g0, g1, g2] = constants->gamma[nf - 1];
256 double b02 = b0 * b0;
257 double a = std::pow(b0 * alpha / (2 * PI), 2 * g0 / b0);
258 double b = (2 * g1 / b0 - b1 * g0 / b02) * alpha / PI;
259 double c = .5 * (std::pow(2 * g1 / b0 - b1 * g0 / b02, 2)
260 + 2 * g2 / b0 - b1 * g1 / b02 - b2 * g0 / (16 * b02)
261 + b1 * b1 * g0 / (2 * b0 * b02)) * std::pow(alpha / PI, 2);
262 return a * (1 + b + c);
263}
MassType
MassType
Definition GeneralEnum.h:661
LOG_ERROR
#define LOG_ERROR(type,...)
Macro for logging error messages and terminating the application.
Definition Logger.h:41
LOG_DEBUG
#define LOG_DEBUG(...)
Macro for logging debug messages.
Definition Logger.h:45
QCDHelper.h
Core QCD utilities for α_s running and quark mass evolution.
BlockSrc
Lightweight view over a set of source blocks.
Definition SourcesView.h:71
DependentBlockManager::addDependentBlock
static void addDependentBlock(std::string name, std::unordered_map< ParameterType, std::vector< std::string > > source_names, ParameterType dest, DepUpdateFunc recalculateFunc)
Adds a DependentBlock with multiple sources coming from various Parameters instances.
Definition DependentBlockManager.cpp:4
Parameters
Block-based parameter repository for one ParameterType namespace.
Definition Parameters.h:351
Parameters::GetInstance
static std::shared_ptr< Parameters > GetInstance(ParameterType id=ParameterType::SM)
Returns the singleton-like repository for a given parameter type.
Definition Parameters.cpp:11
QCDHelper::alpha_s
static double alpha_s(double mu, MassType mass_b_type=MassType::POLE, MassType mass_t_type=MassType::POLE)
Computes the strong coupling constant α_s at scale μ.
Definition QCDHelper.cpp:52
QCDHelper::msbar_mass
static double msbar_mass(int pdg_code, double mu, MassType mass_b_type=MassType::POLE, MassType mass_t_type=MassType::POLE)
Computes the MS-bar running mass of a quark at scale μ.
Definition QCDHelper.cpp:64
QCDHelper::Init
static void Init()
Initializes the QCD dependent block.
Definition QCDHelper.cpp:3
QCDHelper::constants
static QCDConstants * constants
Pointer to QCD constants (to conform with IQCDProvider interface).
Definition QCDHelper.h:132
QCDHelper::get_nf
static int get_nf(double mu, MassType mass_b_type=MassType::POLE, MassType mass_t_type=MassType::POLE)
Returns the number of active flavors n_f at scale μ.
Definition QCDHelper.cpp:166
PI
constexpr double PI
Definition constants.h:7
PI2
constexpr double PI2
Definition constants.h:8
ZETA3
constexpr double ZETA3
Definition constants.h:14
mty::sm_input::m_b
csl::Expr m_b
Definition sm.h:96
mty::sm_input::m_t
csl::Expr m_t
Definition sm.h:98
test_rng.L
L
Definition test_rng.py:52
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
LhaID
Represents an identifier of a LHA element, possibly containing several sub-ids.
Definition LhaID.h:56
ParamId
Composite identifier for a single parameter.
Definition ParamID.h:57
QCDConstants::C_F
static constexpr double C_F
Fundamental Casimir C_F = (N_c^2 - 1) / (2 N_c).
Definition QCDHelper.h:45
QCDConstants::beta
static constexpr std::array< std::array< double, 3 >, 6 > beta
β-function coefficients β_0, β_1, β_2 for nf = 1..6.
Definition QCDHelper.h:54
QCDConstants::gamma
static constexpr std::array< std::array< double, 3 >, 6 > gamma
γ-function coefficients γ_0, γ_1, γ_2 for nf = 1..6.
Definition QCDHelper.h:66