From a0eebeea25ba7c3da89faa6b91d06ad9d15ccaab Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Mon, 17 Jan 2022 14:21:09 +0100 Subject: [PATCH 01/24] Adding N3LO msbar mass running --- src/eko/msbar_masses.py | 51 +++++++++++++++++++++++++++++++---------- 1 file changed, 39 insertions(+), 12 deletions(-) diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py index 6cfef54d9..d13bc5fb1 100644 --- a/src/eko/msbar_masses.py +++ b/src/eko/msbar_masses.py @@ -30,8 +30,8 @@ def msbar_ker_exact(a0, a1, order, nf): ker: float Exact |MSbar| kernel: - ..math: - k_{exact} = e^{\int_{a_s(\mu_{h,0}^2)}^{a_s(\mu^2)} \gamma(a_s) / \beta(a_s) da_s} + .. math:: + k_{exact} = e^{- \int_{a_s(\mu_{h,0}^2)}^{a_s(\mu^2)} \gamma_m(a_s) / \beta(a_s) da_s} """ b_vec = [beta(0, nf)] g_vec = [gamma(0, nf)] @@ -41,6 +41,9 @@ def msbar_ker_exact(a0, a1, order, nf): if order >= 2: b_vec.append(beta(2, nf)) g_vec.append(gamma(2, nf)) + if order >= 3: + b_vec.append(beta(3, nf)) + g_vec.append(gamma(3, nf)) # quad ker def integrand(a, b_vec, g_vec): @@ -83,12 +86,15 @@ def msbar_ker_expanded(a0, a1, order, nf): ker: float Expanded |MSbar| kernel: - ..math: + .. math:: k_{expanded} &= \left (\frac{a_s(\mu^2)}{a_s(\mu_{h,0}^2)} \right )^{c_0} \frac{j_{exp}(a_s(\mu^2))}{j_{exp}(a_s(\mu_{h,0}^2))} \\ - j_{exp}(a_s) &= 1 + a_s \left [ c_1 - b_1 c_0 \right ] - + \frac{a_s^2}{2} - \left [c_2 - c_1 b_1 - b_2 c_0 + b_1^2 c_0 + (c_1 - b_1 c_0)^2 right] + j_{exp}(a_s) &= 1 + a_s \left [ c_1 - b_1 c_0 \right ] \\ + & + \frac{a_s^2}{2} + \left [c_2 - c_1 b_1 - b_2 c_0 + b_1^2 c_0 + (c_1 - b_1 c_0)^2 \right ] \\ + & + \frac{a_s^3}{6} [ -2 b_3 c_0 - b_1^3 c_0 (1 + c_0) (2 + c_0) - 2 b_2 c_1 \\ + & - 3 b_2 c_0 c_1 + b_1^2 (2 + 3 c_0 (2 + c_0)) c_1 + c_1^3 + 3 c_1 c_2 \\ + & + b_1 (b_2 c_0 (4 + 3 c_0) - 3 (1 + c_0) c_1^2 - (2 + 3 c_0) c_2) + 2 c_3 ] """ b0 = beta(0, nf) c0 = gamma(0, nf) / b0 @@ -98,15 +104,36 @@ def msbar_ker_expanded(a0, a1, order, nf): if order >= 1: b1 = b(1, nf) c1 = gamma(1, nf) / b0 - r = c1 - b1 * c0 - num += a1 * r - den += a0 * r + u = c1 - b1 * c0 + num += a1 * u + den += a0 * u if order >= 2: b2 = b(2, nf) c2 = gamma(2, nf) / b0 - r = (c2 - c1 * b1 - b2 * c0 + b1 ** 2 * c0 + (c1 - b1 * c0) ** 2) / 2.0 - num += a1 ** 2 * r - den += a0 ** 2 * r + u = (c2 - c1 * b1 - b2 * c0 + b1 ** 2 * c0 + (c1 - b1 * c0) ** 2) / 2.0 + num += a1 ** 2 * u + den += a0 ** 2 * u + if order >= 3: + b3 = b(3, nf) + c3 = gamma(3, nf) / b0 + u = ( + 1 + / 6 + * ( + -2 * b3 * c0 + - b1 ** 3 * c0 * (1 + c0) * (2 + c0) + - 2 * b2 * c1 + - 3 * b2 * c0 * c1 + + b1 ** 2 * (2 + 3 * c0 * (2 + c0)) * c1 + + c1 ** 3 + + 3 * c1 * c2 + + b1 + * (b2 * c0 * (4 + 3 * c0) - 3 * (1 + c0) * c1 ** 2 - (2 + 3 * c0) * c2) + + 2 * c3 + ) + ) + num += a1 ** 3 * u + den += a0 ** 3 * u return ev_mass * num / den From ae73822286708e8083e1d24f454d1652afb89071 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Mon, 17 Jan 2022 14:21:57 +0100 Subject: [PATCH 02/24] more msbar tests --- tests/benchmark_msbar_evolution.py | 168 +++++++++++++++++------------ tests/test_msbar_masses.py | 4 +- 2 files changed, 100 insertions(+), 72 deletions(-) diff --git a/tests/benchmark_msbar_evolution.py b/tests/benchmark_msbar_evolution.py index 3a93df3b2..df218748a 100644 --- a/tests/benchmark_msbar_evolution.py +++ b/tests/benchmark_msbar_evolution.py @@ -24,77 +24,103 @@ def benchmark_APFEL_msbar_evolution(self): Q2m = np.power([2.0, 4.5, 175], 2) m2 = np.power((1.4, 4.5, 175), 2) apfel_vals_dict = { - 0: np.array( - [ - [1.6046128320144937, 5.82462147889985, 319.1659462836651], - [0.9184216270407894, 3.3338000474743867, 182.67890037277968], - [0.8892735505271812, 3.2279947658871553, 176.88119438599423], - ] - ), - 1: np.array( - [ - [1.7606365126844892, 6.707425163030458, 392.9890591164956], - [0.8227578662769597, 3.134427109507681, 183.6465604399247], - [0.7934449726444709, 3.0227549733595973, 177.10367302092334], - ] - ), - 2: np.array( - [ - [1.8315619347807859, 7.058348563796064, 417.06081590596904], - [0.8078493428925942, 3.113242546931765, 183.76436225206226], - [0.7788276351007536, 3.0014003869088755, 177.16269119698978], - ] - ), + "exact": { + 0: np.array( + [ + [1.6046128320144937, 5.82462147889985, 319.1659462836651], + [0.9184216270407894, 3.3338000474743867, 182.67890037277968], + [0.8892735505271812, 3.2279947658871553, 176.88119438599423], + ] + ), + 1: np.array( + [ + [1.7606365126844892, 6.707425163030458, 392.9890591164956], + [0.8227578662769597, 3.134427109507681, 183.6465604399247], + [0.7934449726444709, 3.0227549733595973, 177.10367302092334], + ] + ), + 2: np.array( + [ + [1.8315619347807859, 7.058348563796064, 416.630860855048], + [0.8078493428925942, 3.113242546931765, 183.76436225206226], + [0.7788276351007536, 3.0014003869088755, 177.16269119698978], + ] + ), + }, + "expanded": { + 0: np.array( + [ + [1.6046128320144941, 5.824621478899853, 319.1659462836651], + [0.9184216270407891, 3.3338000474743863, 182.67890037277968], + [0.8892735505271808, 3.227994765887154, 176.88119438599423], + ] + ), + 1: np.array( + [ + [1.7533055503995305, 6.672122790503439, 390.2255025944903], + [0.8251533585949282, 3.1400827586994855, 183.65075271254497], + [0.7957427000040153, 3.028161864236207, 177.104951823859], + ] + ), + 2: np.array( + [ + [1.8268480938423455, 7.037891257825139, 415.2638988980831], + [0.8084237419306857, 3.1144420987483485, 183.76461377981423], + [0.7793806824526442, 3.002554084550691, 177.16277079680097], + ] + ), + }, } # collect my values - for order in [0, 1, 2]: - as_VFNS = StrongCoupling( - alphas_ref, - scale_ref, - m2, - thresholds_ratios, - order=order, - method="exact", - hqm_scheme="MSBAR", - ) - my_vals = [] - for Q2 in Q2s: - my_masses = [] - for n in [3, 4, 5]: - my_masses.append( - evolve_msbar_mass( - m2[n - 3], - Q2m[n - 3], - strong_coupling=as_VFNS, - config=dict(fact_to_ren=1), - q2_to=Q2, - ) - ) - my_vals.append(my_masses) - # get APFEL numbers - if available else use cache - apfel_vals = apfel_vals_dict[order] - if use_APFEL: - # run apfel - apfel.CleanUp() - apfel.SetTheory("QCD") - apfel.SetPerturbativeOrder(order) - apfel.SetAlphaEvolution("exact") - apfel.SetAlphaQCDRef(alphas_ref, np.sqrt(scale_ref)) - apfel.SetVFNS() - apfel.SetMSbarMasses(*np.sqrt(m2)) - apfel.SetMassScaleReference(*np.sqrt(Q2m)) - apfel.SetRenFacRatio(1) - apfel.InitializeAPFEL() - # collect apfel masses - apfel_vals_cur = [] + for method in ["exact", "expanded"]: + for order in [0, 1, 2]: + as_VFNS = StrongCoupling( + alphas_ref, + scale_ref, + m2, + thresholds_ratios, + order=order, + method=method, + hqm_scheme="MSBAR", + ) + my_vals = [] for Q2 in Q2s: - masses = [] - for n in [4, 5, 6]: - masses.append(apfel.HeavyQuarkMass(n, np.sqrt(Q2))) - apfel_vals_cur.append(masses) - print(apfel_vals_cur) - np.testing.assert_allclose(apfel_vals, np.array(apfel_vals_cur)) - # check myself to APFEL - np.testing.assert_allclose( - apfel_vals, np.sqrt(np.array(my_vals)), rtol=2e-3 - ) + my_masses = [] + for n in [3, 4, 5]: + my_masses.append( + evolve_msbar_mass( + m2[n - 3], + Q2m[n - 3], + strong_coupling=as_VFNS, + config=dict(fact_to_ren=1), + q2_to=Q2, + ) + ) + my_vals.append(my_masses) + # get APFEL numbers - if available else use cache + apfel_vals = apfel_vals_dict[method][order] + if use_APFEL: + # run apfel + apfel.CleanUp() + apfel.SetTheory("QCD") + apfel.SetPerturbativeOrder(order) + apfel.SetAlphaEvolution(method) + apfel.SetAlphaQCDRef(alphas_ref, np.sqrt(scale_ref)) + apfel.SetVFNS() + apfel.SetMSbarMasses(*np.sqrt(m2)) + apfel.SetMassScaleReference(*np.sqrt(Q2m)) + apfel.SetRenFacRatio(1) + apfel.InitializeAPFEL() + # collect apfel masses + apfel_vals_cur = [] + for Q2 in Q2s: + masses = [] + for n in [4, 5, 6]: + masses.append(apfel.HeavyQuarkMass(n, np.sqrt(Q2))) + apfel_vals_cur.append(masses) + print(apfel_vals_cur) + np.testing.assert_allclose(apfel_vals, np.array(apfel_vals_cur)) + # check myself to APFEL + np.testing.assert_allclose( + apfel_vals, np.sqrt(np.array(my_vals)), rtol=2e-3 + ) diff --git a/tests/test_msbar_masses.py b/tests/test_msbar_masses.py index 4721ee02d..70dab67e2 100644 --- a/tests/test_msbar_masses.py +++ b/tests/test_msbar_masses.py @@ -38,7 +38,9 @@ class TestMsbarMasses: def test_compute_msbar_mass(self): # Test solution of msbar(m) = m for method in ["EXA", "EXP"]: - for order in [1, 2]: + # TODO: this will fail until the N3LO alphas is not allowed, + # see N3LO matching PR + for order in [1, 2, 3]: theory_dict.update({"ModEv": method, "PTO": order}) strong_coupling = StrongCoupling.from_dict(theory_dict) From 9d211e36226d4027a382c34ef5f8d8d1f43ab2d7 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Mon, 17 Jan 2022 14:22:53 +0100 Subject: [PATCH 03/24] Fix msbar docs #84 --- doc/source/theory/pQCD.rst | 18 +++++++++++------- src/eko/gamma.py | 2 +- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/doc/source/theory/pQCD.rst b/doc/source/theory/pQCD.rst index 3a9aba25e..2fee66a4a 100644 --- a/doc/source/theory/pQCD.rst +++ b/doc/source/theory/pQCD.rst @@ -97,7 +97,7 @@ For each heavy quark :math:`h` we solve for :math:`m_h`: where the evolved |MSbar| mass is given by: .. math :: - m_{\overline{MS},h}(\mu^2) = m_{h,0} \int_{a_s(\mu_{h,0}^2)}^{a_s(\mu^2)} \frac{\gamma(a_s)}{\beta(a_s)} d a_s + m_{\overline{MS},h}(\mu^2) = m_{h,0} \exp \left[ - \int_{a_s(\mu_{h,0}^2)}^{a_s(\mu^2)} \frac{\gamma_m(a_s)}{\beta(a_s)} d a_s \right ] and :math:`m_{h,0}` is the given initial condition at the scale :math:`\mu_{h,0}`. Here there is a subtle complication since the solution @@ -127,24 +127,28 @@ In doing so EKO takes advantage of the monotony of the |RGE| solution Now, being able to evaluate :math:`a_s(\mu_{h,0}^2)`, there are two ways of solving the previous integral and finally compute the evolved -:math:`m_{\overline{MS},h}`. In fact, the function :math:`\gamma(a_s)` is the +:math:`m_{\overline{MS},h}`. In fact, the function :math:`\gamma_m(a_s)` is the anomalous QCD mass dimension and, as the :math:`\beta` function, it can be evaluated -perturbatively in :math:`a_s` up to :math:`\mathcal{O}(a_s^3)`: +perturbatively in :math:`a_s` up to :math:`\mathcal{O}(a_s^4)`: .. math :: - \gamma(a_s) &= - \sum\limits_{n=0} \gamma_n a_s^{n+1} \\ + \gamma_m(a_s) &= \sum\limits_{n=0} \gamma_{m_n} a_s^{n} \\ -Even here it is useful to define :math:`c_k = \gamma_k/\beta_0, k>0`. +Even here it is useful to define :math:`c_k = \gamma_{m_k}/\beta_0, k \ge 0`. Therefore the two solution strategies are: - ``method = "exact"``: the integral is solved exactly using the expression of - :math:`\beta,\gamma` up to the specified perturbative order + :math:`\beta,\gamma_m` up to the specified perturbative order - ``method = "expanded"``: the integral is approximate by the following expansion: .. math :: m_{\overline{MS},h}(\mu^2) & = m_{h,0} \left ( \frac{a_s(\mu^2)}{a_s(\mu_{h,0}^2)} \right )^{c_0} \frac{j_{exp}(a_s(\mu^2))}{j_{exp}(a_s(\mu_{h,0}^2))} \\ - j_{exp}(a_s) &= 1 + a_s \left [ c_1 - b_1 c_0 \right ] + \frac{a_s^2}{2} \left [c_2 - c_1 b_1 - b_2 c_0 + b_1^2 c_0 + (c_1 - b_1 c_0)^2 \right] + j_{exp}(a_s) &= 1 + a_s \left [ c_1 - b_1 c_0 \right ] \\ + & + \frac{a_s^2}{2} \left [c_2 - c_1 b_1 - b_2 c_0 + b_1^2 c_0 + (c_1 - b_1 c_0)^2 \right] \\ + & + \frac{a_s^3}{6} [ -2 b_3 c_0 - b_1^3 c_0 (1 + c_0) (2 + c_0) - 2 b_2 c_1 \\ + & - 3 b_2 c_0 c_1 + b_1^2 (2 + 3 c_0 (2 + c_0)) c_1 + c_1^3 + 3 c_1 c_2 \\ + & + b_1 (b_2 c_0 (4 + 3 c_0) - 3 (1 + c_0) c_1^2 - (2 + 3 c_0) c_2) + 2 c_3 ] The procedure is iterated on all the heavy quarks, updating the temporary instance diff --git a/src/eko/gamma.py b/src/eko/gamma.py index 468c54693..cafada3c3 100644 --- a/src/eko/gamma.py +++ b/src/eko/gamma.py @@ -66,5 +66,5 @@ def gamma(order, nf): elif order == 3: _gamma = gamma_3(nf) else: - raise ValueError("Gamma coefficients beyond N3LO are not implemented!") + raise ValueError("QCD Gamma coefficients beyond N3LO are not implemented!") return _gamma From 1b7aa2cc0ab02705f72ed5a4a1f9902329eef000 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Tue, 18 Jan 2022 09:30:51 +0100 Subject: [PATCH 04/24] Adding beta3 and allow for N3LO alpha_s --- doc/source/theory/pQCD.rst | 4 ++-- src/eko/beta.py | 31 ++++++++++++++++++++++++++++++- src/eko/strong_coupling.py | 18 ++++++++++++++---- tests/test_beta.py | 13 ++++++++++++- tests/test_strong_coupling.py | 4 ++-- 5 files changed, 60 insertions(+), 10 deletions(-) diff --git a/doc/source/theory/pQCD.rst b/doc/source/theory/pQCD.rst index 2fee66a4a..78fa027dd 100644 --- a/doc/source/theory/pQCD.rst +++ b/doc/source/theory/pQCD.rst @@ -132,9 +132,9 @@ anomalous QCD mass dimension and, as the :math:`\beta` function, it can be evalu perturbatively in :math:`a_s` up to :math:`\mathcal{O}(a_s^4)`: .. math :: - \gamma_m(a_s) &= \sum\limits_{n=0} \gamma_{m_n} a_s^{n} \\ + \gamma_m(a_s) &= \sum\limits_{n=0} \gamma_{m,n} a_s^{n+1} \\ -Even here it is useful to define :math:`c_k = \gamma_{m_k}/\beta_0, k \ge 0`. +Even here it is useful to define :math:`c_k = \gamma_{m,k}/\beta_0, k \ge 0`. Therefore the two solution strategies are: diff --git a/src/eko/beta.py b/src/eko/beta.py index 8dba687be..b7ed021ad 100644 --- a/src/eko/beta.py +++ b/src/eko/beta.py @@ -8,6 +8,7 @@ import numba as nb from eko import constants +from eko.anomalous_dimensions.harmonics import zeta3 @nb.njit("f8(u1)", cache=True) @@ -85,6 +86,32 @@ def beta_2(nf: int): return beta_2 +@nb.njit("f8(u1)", cache=True) +def beta_3(nf: int): + """ + Computes the fourth coefficient of the QCD beta function + + Implements Eq. (3.6) of :cite:`Herzog:2017ohr`. + + Parameters + ---------- + nf : int + number of active flavors + + Returns + ------- + beta_3 : float + fourth coefficient of the QCD beta function :math:`\\beta_3^{n_f}` + """ + beta_3 = ( + 149753.0 / 6.0 + + 3564.0 * zeta3 + nf * (-1078361.0 / 162.0 - 6508.0 / 27.0 * zeta3) + + nf ** 2 * (50065.0 / 162.0 + 6472.0 / 81.0 * zeta3) + + 1093.0 / 729.0 * nf ** 3 + ) + return beta_3 + + @nb.njit("f8(u1,u1)", cache=True) def beta(k, nf): """ @@ -109,8 +136,10 @@ def beta(k, nf): beta_ = beta_1(nf) elif k == 2: beta_ = beta_2(nf) + elif k == 3: + beta_ = beta_3(nf) else: - raise ValueError("Beta coefficients beyond NNLO are not implemented!") + raise ValueError("Beta coefficients beyond N3LO are not implemented!") return beta_ diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index 1109a5273..93ffa2d77 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -72,7 +72,7 @@ def as_expanded(order, as_ref, nf, scale_from, scale_to): as_NLO = as_LO * (1 - b1 * as_LO * np.log(den)) res = as_NLO # NNLO - if order == 2: + if order >= 2: beta2 = beta(2, nf) b2 = beta2 / beta0 res = as_LO * ( @@ -80,6 +80,10 @@ def as_expanded(order, as_ref, nf, scale_from, scale_to): + as_LO * (as_LO - as_ref) * (b2 - b1 ** 2) + as_NLO * b1 * np.log(as_NLO / as_ref) ) + #N3LO + if order >= 3: + pass + # TODO: add N3LO expanded return res @@ -146,8 +150,8 @@ def __init__( raise ValueError(f"alpha_s_ref has to be positive - got {alpha_s_ref}") if scale_ref <= 0: raise ValueError(f"scale_ref has to be positive - got {scale_ref}") - if order not in [0, 1, 2]: - raise NotImplementedError("a_s beyond NNLO is not implemented") + if order not in [0, 1, 2, 3]: + raise NotImplementedError("a_s beyond N3LO is not implemented") self.order = order if method not in ["expanded", "exact"]: raise ValueError(f"Unknown method {method}") @@ -268,6 +272,11 @@ def compute_exact(self, as_ref, nf, scale_from, scale_to): beta2 = beta(2, nf) b2 = beta2 / beta0 b_vec.append(b2) + # N3LO + if self.order >= 3: + beta3 = beta(3, nf) + b3 = beta3 / beta0 + b_vec.append(b3) # integration kernel def rge(_t, a, b_vec): return -(a ** 2) * np.sum([a ** k * b for k, b in enumerate(b_vec)]) @@ -393,7 +402,8 @@ def compute_matching_coeffs_up(mass_scheme): matching_coeffs_down: forward matching coefficient matrix """ - matching_coeffs_up = np.zeros((3, 3)) + matching_coeffs_up = np.zeros((4, 4)) + # TODO: add N3LO parameters if mass_scheme == "MSBAR": matching_coeffs_up[2, 0] = -22.0 / 3.0 matching_coeffs_up[2, 1] = 22.0 / 3.0 diff --git a/tests/test_beta.py b/tests/test_beta.py index a8cf1894f..3ad36e31f 100644 --- a/tests/test_beta.py +++ b/tests/test_beta.py @@ -7,6 +7,7 @@ import pytest from eko import beta +from eko.anomalous_dimensions.harmonics import zeta3 def _flav_test(function): @@ -38,14 +39,24 @@ def test_beta_2(): np.testing.assert_approx_equal(beta.beta_2(5), 4 ** 3 * 9769 / 3456) +def test_beta_3(): + """Test fourth beta function coefficient""" + _flav_test(beta.beta_3) + # from hep-ph/9706430 + np.testing.assert_allclose( + beta.beta_3(5), 4 ** 4 * (11027.0 / 648.0 * zeta3 - 598391.0 / 373248.0) + ) + + def test_beta(): """beta-wrapper""" nf = 3 np.testing.assert_allclose(beta.beta(0, nf), beta.beta_0(nf)) np.testing.assert_allclose(beta.beta(1, nf), beta.beta_1(nf)) np.testing.assert_allclose(beta.beta(2, nf), beta.beta_2(nf)) + np.testing.assert_allclose(beta.beta(3, nf), beta.beta_3(nf)) with pytest.raises(ValueError): - beta.beta(3, 3) + beta.beta(4, 3) def test_b(): diff --git a/tests/test_strong_coupling.py b/tests/test_strong_coupling.py index 20ef41af8..2d689de87 100644 --- a/tests/test_strong_coupling.py +++ b/tests/test_strong_coupling.py @@ -85,7 +85,7 @@ def test_init(self): scale_ref, threshold_holder.area_walls[1:-1], (1.0, 1.0, 1.0), - 3, + 4, ) with pytest.raises(ValueError): StrongCoupling( @@ -127,7 +127,7 @@ def test_ref(self): alphas_ref = 0.118 scale_ref = 91.0 ** 2 for thresh_setup in thresh_setups: - for order in [0, 1, 2]: + for order in [0, 1, 2, 3]: for method in ["exact", "expanded"]: # create sc = StrongCoupling( From 1dbd5ce768f1e8722e8eb4c49b4dccec04cb736e Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Wed, 19 Jan 2022 08:56:04 +0100 Subject: [PATCH 05/24] Add alpha_s N3LO matching coeffs --- doc/source/refs.bib | 13 ++++++++ doc/source/theory/pQCD.rst | 2 +- src/eko/msbar_masses.py | 2 +- src/eko/strong_coupling.py | 50 ++++++++++++++++++++++-------- tests/benchmark_msbar_evolution.py | 2 +- tests/benchmark_strong_coupling.py | 38 ++++++++++++++++------- tests/test_msbar_masses.py | 2 -- 7 files changed, 80 insertions(+), 29 deletions(-) diff --git a/doc/source/refs.bib b/doc/source/refs.bib index a760f287e..f035ee5f9 100644 --- a/doc/source/refs.bib +++ b/doc/source/refs.bib @@ -405,3 +405,16 @@ @article{Gribov:1972ri reportnumber = {IPTI-381-71}, slaccitation = {%%CITATION = SJNCA,15,438;%%} } + +@article{Chetyrkin:1997sg, + author = "Chetyrkin, K. G. and Kniehl, Bernd A. and Steinhauser, M.", + title = "{Strong coupling constant with flavor thresholds at four loops in the MS scheme}", + eprint = "hep-ph/9706430", + archivePrefix = "arXiv", + reportNumber = "MPI-PHT-97-025", + doi = "10.1103/PhysRevLett.79.2184", + journal = "Phys. Rev. Lett.", + volume = "79", + pages = "2184--2187", + year = "1997" +} \ No newline at end of file diff --git a/doc/source/theory/pQCD.rst b/doc/source/theory/pQCD.rst index 78fa027dd..b5af9bfb0 100644 --- a/doc/source/theory/pQCD.rst +++ b/doc/source/theory/pQCD.rst @@ -29,7 +29,7 @@ We implement two different strategies to solve the |RGE|: & \hspace{60pt} \left. + a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2) b_1 \ln\left(a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2)/a_s(\mu_0^2)\right)\right] When the renormalization scale crosses a flavor threshold matching conditions -have to be applied :cite:`Schroder:2005hy,Chetyrkin:2005ia`. +have to be applied :cite:`Chetyrkin:1997sg,Schroder:2005hy,Chetyrkin:2005ia`. Splitting Functions ------------------- diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py index d13bc5fb1..ceedf2a28 100644 --- a/src/eko/msbar_masses.py +++ b/src/eko/msbar_masses.py @@ -32,7 +32,7 @@ def msbar_ker_exact(a0, a1, order, nf): .. math:: k_{exact} = e^{- \int_{a_s(\mu_{h,0}^2)}^{a_s(\mu^2)} \gamma_m(a_s) / \beta(a_s) da_s} - """ + """ # pylint:disable=line-too-long b_vec = [beta(0, nf)] g_vec = [gamma(0, nf)] if order >= 1: diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index 93ffa2d77..e324d7190 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -80,7 +80,7 @@ def as_expanded(order, as_ref, nf, scale_from, scale_to): + as_LO * (as_LO - as_ref) * (b2 - b1 ** 2) + as_NLO * b1 * np.log(as_NLO / as_ref) ) - #N3LO + # N3LO if order >= 3: pass # TODO: add N3LO expanded @@ -368,9 +368,9 @@ def a_s( self.thresholds.thresholds_ratios[seg.nf - shift] ) m_coeffs = ( - compute_matching_coeffs_down(self.hqm_scheme) + compute_matching_coeffs_down(self.hqm_scheme, seg.nf-1) if is_downward_path - else compute_matching_coeffs_up(self.hqm_scheme) + else compute_matching_coeffs_up(self.hqm_scheme, seg.nf) ) fact = 1.0 # shift @@ -383,10 +383,15 @@ def a_s( return final_as -def compute_matching_coeffs_up(mass_scheme): +def compute_matching_coeffs_up(mass_scheme, nf): r""" - Matching coefficients :cite:`Schroder:2005hy,Chetyrkin:2005ia,Vogt:2004ns` at threshold - when moving to a regime with *more* flavors. + Matching coefficients :cite:`Chetyrkin:1997sg,Schroder:2005hy,Chetyrkin:2005ia,Vogt:2004ns` + at threshold when moving to a regime with *more* flavors. + + We follow notation of :cite:`Vogt:2004ns` (eq 2.43) + + The *inverse* |MSbar| values instead are given in :cite:`Chetyrkin:1997sg` (eq 6) + multiplied by a factor of 4 (and 4^2 ...) .. math:: a_s^{(n_l+1)} = a_s^{(n_l)} + \sum\limits_{n=1} (a_s^{(n_l)})^n @@ -396,6 +401,8 @@ def compute_matching_coeffs_up(mass_scheme): ---------- mass_scheme: Heavy quark mass scheme: "POLE" or "MSBAR" + nf: + number of active flavors in the lower patch Returns ------- @@ -403,23 +410,34 @@ def compute_matching_coeffs_up(mass_scheme): forward matching coefficient matrix """ matching_coeffs_up = np.zeros((4, 4)) - # TODO: add N3LO parameters if mass_scheme == "MSBAR": - matching_coeffs_up[2, 0] = -22.0 / 3.0 - matching_coeffs_up[2, 1] = 22.0 / 3.0 + matching_coeffs_up[2, 0] = -22.0 / 9.0 + # TODO: here apfel and :cite:`Chetyrkin:1997sg` (eq 6) do not agree... + # Apfel has: matching_coeffs_up[2, 1] = 22.0 / 3.0 + + # c30 = -d30 + matching_coeffs_up[3, 0] = -62.2116 + 5.4177 * nf + # c31 = -d31 + 5 c11 * c20 + matching_coeffs_up[3, 1] = 243.444 - 10.4074 * nf + elif mass_scheme == "POLE": matching_coeffs_up[2, 0] = 14.0 / 3.0 - matching_coeffs_up[2, 1] = 38.0 / 3.0 + matching_coeffs_up[3, 0] = 340.729 - 16.7981 * nf + matching_coeffs_up[3, 1] = 8941.0 / 27.0 - 409.0 / 27.0 * nf matching_coeffs_up[1, 1] = 4.0 / 3.0 * constants.TR + matching_coeffs_up[2, 1] = 38.0 / 3.0 matching_coeffs_up[2, 2] = 4.0 / 9.0 + matching_coeffs_up[3, 2] = 511.0 / 9.0 + matching_coeffs_up[3, 3] = 8.0 / 27.0 + return matching_coeffs_up # inversion of the matching coefficients -def compute_matching_coeffs_down(mass_scheme): +def compute_matching_coeffs_down(mass_scheme, nf): """ - Matching coefficients :cite:`Schroder:2005hy` :cite:`Chetyrkin:2005ia` at threshold + Matching coefficients :cite:`Chetyrkin:1997sg,Schroder:2005hy,Chetyrkin:2005ia` at threshold when moving to a regime with *less* flavors. This is the perturbative inverse of :data:`matching_coeffs_up` and has been obtained via @@ -442,16 +460,22 @@ def compute_matching_coeffs_down(mass_scheme): ---------- mass_scheme: Heavy quark mass scheme: "POLE" or "MSBAR" + nf: + number of active flavors in the lower patch Returns ------- matching_coeffs_down: downward matching coefficient matrix """ - _c = compute_matching_coeffs_up(mass_scheme) + _c = compute_matching_coeffs_up(mass_scheme, nf) matching_coeffs_down = np.zeros_like(_c) matching_coeffs_down[1, 1] = -_c[1, 1] matching_coeffs_down[2, 0] = -_c[2, 0] matching_coeffs_down[2, 1] = -_c[2, 1] matching_coeffs_down[2, 2] = 2.0 * _c[1, 1] ** 2 - _c[2, 2] + matching_coeffs_down[3, 0] = -_c[3, 0] + matching_coeffs_down[3, 1] = 5 * _c[1, 1] * _c[2, 0] - _c[3, 1] + matching_coeffs_down[3, 2] = 5 * _c[1, 1] * _c[2, 1] - _c[3, 2] + matching_coeffs_down[3, 3] = -5 * _c[1, 1] ** 3 + 5 * _c[1, 1] * _c[2, 2] - _c[3, 3] return matching_coeffs_down diff --git a/tests/benchmark_msbar_evolution.py b/tests/benchmark_msbar_evolution.py index df218748a..1bf939177 100644 --- a/tests/benchmark_msbar_evolution.py +++ b/tests/benchmark_msbar_evolution.py @@ -122,5 +122,5 @@ def benchmark_APFEL_msbar_evolution(self): np.testing.assert_allclose(apfel_vals, np.array(apfel_vals_cur)) # check myself to APFEL np.testing.assert_allclose( - apfel_vals, np.sqrt(np.array(my_vals)), rtol=2e-3 + apfel_vals, np.sqrt(np.array(my_vals)), rtol=5e-3 ) diff --git a/tests/benchmark_strong_coupling.py b/tests/benchmark_strong_coupling.py index 01813a5c7..e42a73b64 100644 --- a/tests/benchmark_strong_coupling.py +++ b/tests/benchmark_strong_coupling.py @@ -165,10 +165,18 @@ def benchmark_pegasus_ffns(self): 0.009245273177012448, ] ), + 3: np.array( + [ + 0.02224042559454458, + 0.014988806602725375, + 0.01140950716164355, + 0.009245168439298405 + ] + ), } # collect my values threshold_holder = thresholds.ThresholdsAtlas.ffns(nf) - for order in [0, 1, 2]: + for order in [0, 1, 2, 3]: as_FFNS = StrongCoupling( alphas_ref, scale_ref, @@ -200,11 +208,11 @@ def benchmark_pegasus_ffns(self): ) # print(pegasus_vals_cur) np.testing.assert_allclose( - pegasus_vals, np.array(pegasus_vals_cur), err_msg="order=%d" % order + pegasus_vals, np.array(pegasus_vals_cur), err_msg=f"order={order}" ) - # check myself to APFEL + # check myself to PEGASUS np.testing.assert_allclose( - pegasus_vals, np.array(my_vals), rtol=1.5e-7, err_msg="order=%d" % order + pegasus_vals, np.array(my_vals), rtol=1.5e-7, err_msg=f"order={order}" ) def benchmark_APFEL_vfns(self): @@ -474,18 +482,18 @@ def benchmark_APFEL_vfns_threshold(self): np.testing.assert_allclose(apfel_vals, np.array(my_vals)) def benchmark_APFEL_vfns_msbar(self): - Q2s = np.power([30, 96, 150], 2) + Q2s = np.power([3, 96, 150], 2) alphas_ref = 0.118 scale_ref = 91.0 ** 2 thresholds_ratios = np.power((1.0, 1.0, 1.0), 2) - Q2m = np.power([2.0, 2.0, 175], 2) - m2 = np.power((1.4, 2.0, 175), 2) + Q2m = np.power([2.0, 4.0, 175], 2) + m2 = np.power((1.4, 4.0, 175), 2) apfel_vals_dict = { 1: np.array( - [0.011285356789817751, 0.009315017128518715, 0.00873290495922747] + [0.01985110421500437, 0.009315017128518715, 0.00873290495922747] ), 2: np.array( - [0.011292233300595585, 0.009314871933701218, 0.008731890252153528] + [0.02005213805142418, 0.009314871933701218, 0.008731890252153528] ), } # collect my values @@ -524,7 +532,7 @@ def benchmark_APFEL_vfns_msbar(self): print(apfel_vals_cur) np.testing.assert_allclose(apfel_vals, np.array(apfel_vals_cur)) # check myself to APFEL - np.testing.assert_allclose(apfel_vals, np.array(my_vals), rtol=2e-4) + np.testing.assert_allclose(apfel_vals, np.array(my_vals)) def _get_Lambda2_LO(self, as_ref, scale_ref, nf): """Transformation to Lambda_QCD""" @@ -675,8 +683,16 @@ def benchmark_lhapdf_exact(self): 0.009214183512196827, ] ), + 3: np.array( + [ + 0.025955428065113105, + 0.015862752802768762, + 0.011614793662449371, + 0.009214009540864635 + ] + ), } - for order in range(2 + 1): + for order in range(3 + 1): sc = StrongCoupling( alphas_ref, scale_ref, diff --git a/tests/test_msbar_masses.py b/tests/test_msbar_masses.py index 70dab67e2..5bfee6748 100644 --- a/tests/test_msbar_masses.py +++ b/tests/test_msbar_masses.py @@ -38,8 +38,6 @@ class TestMsbarMasses: def test_compute_msbar_mass(self): # Test solution of msbar(m) = m for method in ["EXA", "EXP"]: - # TODO: this will fail until the N3LO alphas is not allowed, - # see N3LO matching PR for order in [1, 2, 3]: theory_dict.update({"ModEv": method, "PTO": order}) strong_coupling = StrongCoupling.from_dict(theory_dict) From a2c6b28599fe8deeaec5d62e777cb467e29b3cd3 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Wed, 19 Jan 2022 17:39:06 +0100 Subject: [PATCH 06/24] Fix alphas matching conditions and remove old reference --- doc/source/refs.bib | 13 ------------- doc/source/theory/pQCD.rst | 2 +- src/eko/strong_coupling.py | 20 ++++++++++---------- 3 files changed, 11 insertions(+), 24 deletions(-) diff --git a/doc/source/refs.bib b/doc/source/refs.bib index f035ee5f9..a760f287e 100644 --- a/doc/source/refs.bib +++ b/doc/source/refs.bib @@ -405,16 +405,3 @@ @article{Gribov:1972ri reportnumber = {IPTI-381-71}, slaccitation = {%%CITATION = SJNCA,15,438;%%} } - -@article{Chetyrkin:1997sg, - author = "Chetyrkin, K. G. and Kniehl, Bernd A. and Steinhauser, M.", - title = "{Strong coupling constant with flavor thresholds at four loops in the MS scheme}", - eprint = "hep-ph/9706430", - archivePrefix = "arXiv", - reportNumber = "MPI-PHT-97-025", - doi = "10.1103/PhysRevLett.79.2184", - journal = "Phys. Rev. Lett.", - volume = "79", - pages = "2184--2187", - year = "1997" -} \ No newline at end of file diff --git a/doc/source/theory/pQCD.rst b/doc/source/theory/pQCD.rst index b5af9bfb0..78fa027dd 100644 --- a/doc/source/theory/pQCD.rst +++ b/doc/source/theory/pQCD.rst @@ -29,7 +29,7 @@ We implement two different strategies to solve the |RGE|: & \hspace{60pt} \left. + a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2) b_1 \ln\left(a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2)/a_s(\mu_0^2)\right)\right] When the renormalization scale crosses a flavor threshold matching conditions -have to be applied :cite:`Chetyrkin:1997sg,Schroder:2005hy,Chetyrkin:2005ia`. +have to be applied :cite:`Schroder:2005hy,Chetyrkin:2005ia`. Splitting Functions ------------------- diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index e324d7190..793f790bd 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -368,7 +368,7 @@ def a_s( self.thresholds.thresholds_ratios[seg.nf - shift] ) m_coeffs = ( - compute_matching_coeffs_down(self.hqm_scheme, seg.nf-1) + compute_matching_coeffs_down(self.hqm_scheme, seg.nf - 1) if is_downward_path else compute_matching_coeffs_up(self.hqm_scheme, seg.nf) ) @@ -385,12 +385,12 @@ def a_s( def compute_matching_coeffs_up(mass_scheme, nf): r""" - Matching coefficients :cite:`Chetyrkin:1997sg,Schroder:2005hy,Chetyrkin:2005ia,Vogt:2004ns` + Matching coefficients :cite:`Schroder:2005hy,Chetyrkin:2005ia,Vogt:2004ns` at threshold when moving to a regime with *more* flavors. - We follow notation of :cite:`Vogt:2004ns` (eq 2.43) + We follow notation of :cite:`Vogt:2004ns` (eq 2.43) for POLE masses - The *inverse* |MSbar| values instead are given in :cite:`Chetyrkin:1997sg` (eq 6) + The *inverse* |MSbar| values instead are given in :cite:`Schroder:2005hy` (eq 3.1) multiplied by a factor of 4 (and 4^2 ...) .. math:: @@ -412,23 +412,23 @@ def compute_matching_coeffs_up(mass_scheme, nf): matching_coeffs_up = np.zeros((4, 4)) if mass_scheme == "MSBAR": matching_coeffs_up[2, 0] = -22.0 / 9.0 - # TODO: here apfel and :cite:`Chetyrkin:1997sg` (eq 6) do not agree... - # Apfel has: matching_coeffs_up[2, 1] = 22.0 / 3.0 + matching_coeffs_up[2, 1] = 22.0 / 3.0 # c30 = -d30 matching_coeffs_up[3, 0] = -62.2116 + 5.4177 * nf # c31 = -d31 + 5 c11 * c20 - matching_coeffs_up[3, 1] = 243.444 - 10.4074 * nf + matching_coeffs_up[3, 1] = 365.0 / 3.0 - 67.0 / 9.0 * nf + matching_coeffs_up[3, 2] = 109.0 / 3.0 + 16.0 / 9.0 * nf elif mass_scheme == "POLE": matching_coeffs_up[2, 0] = 14.0 / 3.0 + matching_coeffs_up[2, 1] = 38.0 / 3.0 matching_coeffs_up[3, 0] = 340.729 - 16.7981 * nf matching_coeffs_up[3, 1] = 8941.0 / 27.0 - 409.0 / 27.0 * nf + matching_coeffs_up[3, 2] = 511.0 / 9.0 matching_coeffs_up[1, 1] = 4.0 / 3.0 * constants.TR - matching_coeffs_up[2, 1] = 38.0 / 3.0 matching_coeffs_up[2, 2] = 4.0 / 9.0 - matching_coeffs_up[3, 2] = 511.0 / 9.0 matching_coeffs_up[3, 3] = 8.0 / 27.0 return matching_coeffs_up @@ -437,7 +437,7 @@ def compute_matching_coeffs_up(mass_scheme, nf): # inversion of the matching coefficients def compute_matching_coeffs_down(mass_scheme, nf): """ - Matching coefficients :cite:`Chetyrkin:1997sg,Schroder:2005hy,Chetyrkin:2005ia` at threshold + Matching coefficients :cite:`Schroder:2005hy,Chetyrkin:2005ia` at threshold when moving to a regime with *less* flavors. This is the perturbative inverse of :data:`matching_coeffs_up` and has been obtained via From b198af02be581814871e318cd99d190b5715671e Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Wed, 19 Jan 2022 21:51:05 +0100 Subject: [PATCH 07/24] Adding N3LO expansion --- doc/source/theory/pQCD.rst | 23 ++++++++++++++---- src/eko/strong_coupling.py | 45 ++++++++++++++++++++++++++++++----- tests/test_strong_coupling.py | 33 +++++++++++++++++++++++++ 3 files changed, 91 insertions(+), 10 deletions(-) diff --git a/doc/source/theory/pQCD.rst b/doc/source/theory/pQCD.rst index 78fa027dd..47bba9bc8 100644 --- a/doc/source/theory/pQCD.rst +++ b/doc/source/theory/pQCD.rst @@ -23,14 +23,29 @@ We implement two different strategies to solve the |RGE|: - ``method="expanded"``: using approximate solutions: .. math :: - a^{\text{LO}}_s(\mu_R^2) &= \frac{a_s(\mu_0^2)}{1 + a_s(\mu_0^2) \beta_0 \ln(\mu_R^2/\mu_0^2)} \\ - a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2) &= a^{\text{LO}}_s(\mu_R^2)-b_1 \left[a^{\text{LO}}_s(\mu_R^2)\right]^2 \ln\left(1+a_s(\mu_0^2) \beta_0 \ln(\mu_R^2/\mu_0^2)\right) \\ + a^{\text{LO}}_s(\mu_R^2) &= \frac{a_s(\mu_0^2)}{1 + a_s(\mu_0^2) \beta_0 L_{\mu}} \\ + a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2) &= a^{\text{LO}}_s(\mu_R^2)-b_1 \left[a^{\text{LO}}_s(\mu_R^2)\right]^2 \ln\left(1+a_s(\mu_0^2) \beta_0 L_{\mu}\right) \\ a^{\text{NNLO}}_{s,\text{exp}}(\mu_R^2) &= a^{\text{LO}}_s(\mu_R^2)\left[1 + a^{\text{LO}}_s(\mu_R^2)\left(a^{\text{LO}}_s(\mu_R^2) - a_s(\mu_0^2)\right)(b_2 - b_1^2) \right.\\ - & \hspace{60pt} \left. + a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2) b_1 \ln\left(a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2)/a_s(\mu_0^2)\right)\right] + & \hspace{60pt} \left. + a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2) b_1 \ln\left(a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2)/a_s(\mu_0^2)\right)\right] \\ + a^{\text{N3LO}}_{s,\text{exp}}(\mu_R^2) &= a^{\text{NNLO}}_s(\mu_R^2) + \frac{a^{\text{LO}}_s(\mu_R^2)^4}{2 b_0^3 \left(1 + b_0 L_{\mu} a^{\text{LO}}_s(\mu_R^2)\right)^4} \left\{ \right. \\ + & -2 b_1^3 L_{\text{LO}}^3 + 5 b_1^3 L_{\alpha_s}^2 + 2 b_1^3 L_{\alpha_s}^3 + b_1^3 L_{\text{LO}}^2 \left(5 + 6 L_{\alpha_s} \right) \\ + & + 2 b_0 b_1 L_{\alpha_s} \left[ b_2 + 2 \left(b_1^2 - b_0 b_2 \right) L_{\mu} a^{\text{LO}}_s(\mu_R^2) \right] \\ + & - b_0^2 L_{\mu} a^{\text{LO}}_s(\mu_R^2) \left[ -2 b_1 b_2 + 2 b_0 b_3 + \left( b_1^3 - 2 b_0 b_1 b_2 + b_0^2 b_3 \right) L_{\mu} a^{\text{LO}}_s(\mu_R^2) \right] \\ + & - 2 b_1 L_{\text{LO}} \left[ 5 b_1^2 L_{\alpha_s} + 3 b_1^2 L_{\alpha_s}^2 + b_0 \left[b_2 + 2 \left(b_1^2 - b_0 b_2\right) L_{\mu} a^{\text{LO}}_s(\mu_R^2)\right] \right ] \\ + & \left. \right\} + +being: + +.. math :: + L_{\mu} &= \ln(\mu_R^2/\mu_0^2) \\ + L_{\text{LO}} &= \ln(a^{\text{LO}}_s(\mu_R^2)) \\ + L_{\alpha_s} &= \ln \left( \frac{a^{\text{LO}}_s(\mu_R^2)}{(1 + b_0 L_{\mu} a^{\text{LO}}_s(\mu_R^2)} \right) \\ When the renormalization scale crosses a flavor threshold matching conditions have to be applied :cite:`Schroder:2005hy,Chetyrkin:2005ia`. - +In particular the matching when passing form :math:`n_f` to :math:`n_f-1` schemes +is presented in equation 3.1 of :cite:`Schroder:2005hy` for |MSbar| masses, while the +same expression for POLE masses is reported in Appendix A. Splitting Functions ------------------- diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index 793f790bd..e6caf4ef7 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -80,11 +80,44 @@ def as_expanded(order, as_ref, nf, scale_from, scale_to): + as_LO * (as_LO - as_ref) * (b2 - b1 ** 2) + as_NLO * b1 * np.log(as_NLO / as_ref) ) - # N3LO + # N3LO expansion is taken from Luca Rottoli and it's simpler + # than the APFEL expansion which is more accurate if order >= 3: - pass - # TODO: add N3LO expanded - + b3 = beta(3, nf) / beta0 + log_fact = np.log(as_LO / (1 + beta0 * lmu * as_LO)) + res += ( + as_LO ** 4 + / (2 * beta0 ** 3 * (1 + beta0 * lmu * as_LO) ** 4) + * ( + -2 * b1 ** 3 * np.log(as_LO) ** 3 + + 5 * b1 ** 3 * log_fact ** 2 + + 2 * b1 ** 3 * log_fact ** 3 + + b1 ** 3 * np.log(as_LO) ** 2 * (5 + 6 * log_fact) + + 2 + * beta0 + * b1 + * log_fact + * (b2 + 2 * (b1 ** 2 - beta0 * b2) * lmu * as_LO) + - beta0 ** 2 + * lmu + * as_LO + * ( + -2 * b1 * b2 + + 2 * beta0 * b3 + + (b1 ** 3 - 2 * beta0 * b1 * b2 + beta0 ** 2 * b3) + * lmu + * as_LO + ) + - 2 + * b1 + * np.log(as_LO) + * ( + 5 * b1 ** 2 * log_fact + + 3 * b1 ** 2 * log_fact ** 2 + + beta0 * (b2 + 2 * (b1 ** 2 - beta0 * b2) * lmu * as_LO) + ) + ) + ) return res @@ -445,8 +478,8 @@ def compute_matching_coeffs_down(mass_scheme, nf): .. code-block:: Mathematica Module[{f, g, l, sol}, - f[a_] := a + Sum[d[n, k]*L^k*a^(1 + n), {n, 3}, {k, 0, n}]; - g[a_] := a + Sum[c[n, k]*L^k*a^(1 + n), {n, 3}, {k, 0, n}] /. {c[1, 0] -> 0}; + f[a_] := a + Sum[d[n, k)*L^k*a^(1 + n), {n, 3}, {k, 0, n}]; + g[a_] := a + Sum[c[n, k)*L^k*a^(1 + n), {n, 3}, {k, 0, n}] /. {c[1, 0] -> 0}; l = CoefficientList[Normal@Series[f[g[a]], {a, 0, 5}], {a, L}]; sol = First@ Solve[{l[[3]] == 0, l[[4]] == 0, l[[5]] == 0}, diff --git a/tests/test_strong_coupling.py b/tests/test_strong_coupling.py index 2d689de87..6fb231a7e 100644 --- a/tests/test_strong_coupling.py +++ b/tests/test_strong_coupling.py @@ -163,3 +163,36 @@ def test_exact_LO(self): np.testing.assert_allclose( sc_expanded.a_s(q2), sc_exact.a_s(q2), rtol=5e-4 ) + + def benchmark_expanded_n3lo(self): + """test N3LO - NNLO expansion with some references value from Mathematica""" + Q2 = 100 ** 2 + # use a big alpha_s to enlarge th difference + alphas_ref = 0.9 + scale_ref = 90 ** 2 + m2c = 2 + m2b = 25 + m2t = 30625 + threshold_list = [m2c, m2b, m2t] + # Reference values from Mathematica cache + mathematica_val = -0.000101654 + # collect my values + as_NNLO = StrongCoupling( + alphas_ref, + scale_ref, + threshold_list, + (1.0, 1.0, 1.0), + order=2, + method="expanded", + ) + as_N3LO = StrongCoupling( + alphas_ref, + scale_ref, + threshold_list, + (1.0, 1.0, 1.0), + order=3, + method="expanded", + ) + np.testing.assert_allclose( + mathematica_val, as_N3LO.a_s(Q2) - as_NNLO.a_s(Q2), rtol=3e-6 + ) From 49b4841faa86b83e498b3552c7bef1bdd05408f6 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Thu, 20 Jan 2022 11:13:05 +0100 Subject: [PATCH 08/24] Fix some typos --- doc/source/theory/pQCD.rst | 3 ++- src/eko/strong_coupling.py | 4 ++-- tests/test_strong_coupling.py | 5 ++--- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/source/theory/pQCD.rst b/doc/source/theory/pQCD.rst index 47bba9bc8..e856f37a5 100644 --- a/doc/source/theory/pQCD.rst +++ b/doc/source/theory/pQCD.rst @@ -45,7 +45,8 @@ When the renormalization scale crosses a flavor threshold matching conditions have to be applied :cite:`Schroder:2005hy,Chetyrkin:2005ia`. In particular the matching when passing form :math:`n_f` to :math:`n_f-1` schemes is presented in equation 3.1 of :cite:`Schroder:2005hy` for |MSbar| masses, while the -same expression for POLE masses is reported in Appendix A. +same expression for POLE masses is reported in Appendix A. + Splitting Functions ------------------- diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index e6caf4ef7..859540b41 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -478,8 +478,8 @@ def compute_matching_coeffs_down(mass_scheme, nf): .. code-block:: Mathematica Module[{f, g, l, sol}, - f[a_] := a + Sum[d[n, k)*L^k*a^(1 + n), {n, 3}, {k, 0, n}]; - g[a_] := a + Sum[c[n, k)*L^k*a^(1 + n), {n, 3}, {k, 0, n}] /. {c[1, 0] -> 0}; + f[a_] := a + Sum[d[n, k]*L^k*a^(1 + n), {n, 3}, {k, 0, n}]; + g[a_] := a + Sum[c[n, k]*L^k*a^(1 + n), {n, 3}, {k, 0, n}] /. {c[1, 0] -> 0}; l = CoefficientList[Normal@Series[f[g[a]], {a, 0, 5}], {a, L}]; sol = First@ Solve[{l[[3]] == 0, l[[4]] == 0, l[[5]] == 0}, diff --git a/tests/test_strong_coupling.py b/tests/test_strong_coupling.py index 6fb231a7e..98408bd11 100644 --- a/tests/test_strong_coupling.py +++ b/tests/test_strong_coupling.py @@ -165,16 +165,15 @@ def test_exact_LO(self): ) def benchmark_expanded_n3lo(self): - """test N3LO - NNLO expansion with some references value from Mathematica""" + """test N3LO - NNLO expansion with some reference value from Mathematica""" Q2 = 100 ** 2 - # use a big alpha_s to enlarge th difference + # use a big alpha_s to enlarge the difference alphas_ref = 0.9 scale_ref = 90 ** 2 m2c = 2 m2b = 25 m2t = 30625 threshold_list = [m2c, m2b, m2t] - # Reference values from Mathematica cache mathematica_val = -0.000101654 # collect my values as_NNLO = StrongCoupling( From 95fd2c69bbd6b2453bbc4c04e4bbe8a4eb825867 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Tue, 25 Jan 2022 14:43:30 +0100 Subject: [PATCH 09/24] Fix msbar config --- src/eko/msbar_masses.py | 8 ++++---- src/eko/runner.py | 6 +++--- tests/benchmark_msbar_evolution.py | 2 +- tests/test_msbar_masses.py | 6 +++--- 4 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py index ceedf2a28..ba2463eb1 100644 --- a/src/eko/msbar_masses.py +++ b/src/eko/msbar_masses.py @@ -172,8 +172,8 @@ def msbar_ker_dispatcher(q2_to, q2m_ref, strong_coupling, fact_to_ren, nf): def evolve_msbar_mass( m2_ref, q2m_ref, + config, nf_ref=None, - config=None, strong_coupling=None, q2_to=None, ): @@ -189,11 +189,11 @@ def evolve_msbar_mass( squared initial mass reference q2m_ref: float squared initial scale - nf_ref: int, optional (not used when q2_to is given) - number of active flavours at the scale q2m_ref, where the solution is searched config: dict |MSbar| configuration dictionary - strong_coupling: eko.strong_coupling.StrongCoupling + nf_ref: int, optional (not used when q2_to is given) + number of active flavours at the scale q2m_ref, where the solution is searched + strong_coupling: eko.strong_coupling.StrongCoupling, optional Instance of :class:`~eko.strong_coupling.StrongCoupling` able to generate a_s for any q q2_to: float, optional diff --git a/src/eko/runner.py b/src/eko/runner.py index c56efb5c5..afaf475e8 100644 --- a/src/eko/runner.py +++ b/src/eko/runner.py @@ -239,11 +239,11 @@ def compute_msbar_mass(theory_card): # len(masses[q2m_ref > masses]) + 3 is the nf at the given reference scale if nf_target != len(masses[q2m_ref > masses]) + 3: q2_to = masses[q_idx + shift] - m2_ref = evolve_msbar_mass(m2_ref, q2m_ref, config=config, q2_to=q2_to) + m2_ref = evolve_msbar_mass(m2_ref, q2m_ref, config, q2_to=q2_to) q2m_ref = q2_to # now solve the RGE - masses[q_idx] = evolve_msbar_mass(m2_ref, q2m_ref, nf_target, config=config) + masses[q_idx] = evolve_msbar_mass(m2_ref, q2m_ref, config, nf_target) # Check the msbar ordering for m2_msbar, hq in zip(masses[:-1], quark_names[4:]): @@ -251,7 +251,7 @@ def compute_msbar_mass(theory_card): m2_ref = np.power(theory_card[f"m{hq}"], 2) config["thr_masses"] = masses # check that m_msbar_hq < msbar_hq+1 (m_msbar_hq) - m2_test = evolve_msbar_mass(m2_ref, q2m_ref, config=config, q2_to=m2_msbar) + m2_test = evolve_msbar_mass(m2_ref, q2m_ref, config, q2_to=m2_msbar) if m2_msbar > m2_test: raise ValueError( "The MSBAR masses do not preserve the correct ordering,\ diff --git a/tests/benchmark_msbar_evolution.py b/tests/benchmark_msbar_evolution.py index 1bf939177..cd57fc850 100644 --- a/tests/benchmark_msbar_evolution.py +++ b/tests/benchmark_msbar_evolution.py @@ -91,8 +91,8 @@ def benchmark_APFEL_msbar_evolution(self): evolve_msbar_mass( m2[n - 3], Q2m[n - 3], + dict(fact_to_ren=1), strong_coupling=as_VFNS, - config=dict(fact_to_ren=1), q2_to=Q2, ) ) diff --git a/tests/test_msbar_masses.py b/tests/test_msbar_masses.py index 5bfee6748..2afe52913 100644 --- a/tests/test_msbar_masses.py +++ b/tests/test_msbar_masses.py @@ -53,10 +53,10 @@ def test_compute_msbar_mass(self): evolve_msbar_mass( m2_ref, Q2m_ref, - strong_coupling=strong_coupling, - config=dict( + dict( fact_to_ren=theory_dict["fact_to_ren_scale_ratio"] ), + strong_coupling=strong_coupling, q2_to=m2_computed[nf - 3], ) ) @@ -87,8 +87,8 @@ def test_compute_msbar_mass_VFNS(self): evolve_msbar_mass( m2_ref, Q2m_ref, + dict(fact_to_ren=theory_dict["fact_to_ren_scale_ratio"]), strong_coupling=strong_coupling, - config=dict(fact_to_ren=theory_dict["fact_to_ren_scale_ratio"]), q2_to=m2_computed[nf - 3], ) ) From 9e0c6a34db5a765186b15af03327b5a2c351ceef Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Tue, 25 Jan 2022 14:45:22 +0100 Subject: [PATCH 10/24] Some pylint fixes --- src/ekomark/benchmark/external/apfel_utils.py | 4 ++-- src/ekomark/benchmark/runner.py | 2 +- src/ekomark/data/operators.py | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ekomark/benchmark/external/apfel_utils.py b/src/ekomark/benchmark/external/apfel_utils.py index 39d58cf9a..dbc94a9de 100644 --- a/src/ekomark/benchmark/external/apfel_utils.py +++ b/src/ekomark/benchmark/external/apfel_utils.py @@ -72,14 +72,14 @@ def compute_apfel_data( # init evolution apfel.InitializeAPFEL() - print("Loading APFEL took %f s" % (time.perf_counter() - apf_start)) + print(f"Loading APFEL took {(time.perf_counter() - apf_start)} s") # Run apf_tabs = {} for q2 in operators["Q2grid"]: apfel.EvolveAPFEL(theory["Q0"], np.sqrt(q2)) - print("Executing APFEL took %f s" % (time.perf_counter() - apf_start)) + print(f"Executing APFEL took {(time.perf_counter() - apf_start)} s") tab = {} for pid in br.flavor_basis_pids: diff --git a/src/ekomark/benchmark/runner.py b/src/ekomark/benchmark/runner.py index 1b23db6f9..b479e3b0d 100644 --- a/src/ekomark/benchmark/runner.py +++ b/src/ekomark/benchmark/runner.py @@ -87,7 +87,7 @@ def run_me(self, theory, ocard, _pdf): else: # load print(f"Using cached eko data: {os.path.relpath(path,os.getcwd())}") - with open(path) as o: + with open(path, encoding="utf-8") as o: out = eko.output.Output.load_yaml(o) if self.plot_operator: diff --git a/src/ekomark/data/operators.py b/src/ekomark/data/operators.py index b25b6f919..12bd14309 100644 --- a/src/ekomark/data/operators.py +++ b/src/ekomark/data/operators.py @@ -66,7 +66,7 @@ def build(update=None): if update is None: update = {} for c in cartesian_product(update): - card = dict() + card = {} card.update(c) cards.append(card) return cards From 297510a00192375cdf796844747fe5e197d2d568 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Wed, 26 Jan 2022 15:35:36 +0100 Subject: [PATCH 11/24] Adding docstring and numba to gamma --- doc/source/refs.bib | 13 ++++ src/eko/anomalous_dimensions/harmonics.py | 1 + src/eko/gamma.py | 74 ++++++++++++++++++++--- src/eko/msbar_masses.py | 3 +- src/eko/strong_coupling.py | 22 +++---- 5 files changed, 89 insertions(+), 24 deletions(-) diff --git a/doc/source/refs.bib b/doc/source/refs.bib index 347b62de3..e9fdcec6e 100644 --- a/doc/source/refs.bib +++ b/doc/source/refs.bib @@ -419,3 +419,16 @@ @article{Bertone:2013vaa pages = {1647--1668}, year = {2014} } + +@article{Vermaseren:1997fq, + author = "Vermaseren, J. A. M. and Larin, S. A. and van Ritbergen, T.", + title = "{The four loop quark mass anomalous dimension and the invariant quark mass}", + eprint = "hep-ph/9703284", + archivePrefix = "arXiv", + reportNumber = "UM-TH-97-03, NIKHEF-97-012", + doi = "10.1016/S0370-2693(97)00660-6", + journal = "Phys. Lett. B", + volume = "405", + pages = "327--333", + year = "1997" +} \ No newline at end of file diff --git a/src/eko/anomalous_dimensions/harmonics.py b/src/eko/anomalous_dimensions/harmonics.py index 4aae320f1..6270f0f0b 100644 --- a/src/eko/anomalous_dimensions/harmonics.py +++ b/src/eko/anomalous_dimensions/harmonics.py @@ -13,6 +13,7 @@ zeta2 = scipy.special.zeta(2) zeta3 = scipy.special.zeta(3) zeta4 = scipy.special.zeta(4) +zeta5 = scipy.special.zeta(5) @nb.njit("c16(c16,u1)", cache=True) diff --git a/src/eko/gamma.py b/src/eko/gamma.py index 4c974e123..aa7e1346d 100644 --- a/src/eko/gamma.py +++ b/src/eko/gamma.py @@ -4,26 +4,83 @@ See :doc:`pQCD ingredients `. """ +import numba as nb -import scipy.special - -from .anomalous_dimensions.harmonics import zeta3, zeta4 +from .anomalous_dimensions.harmonics import zeta3, zeta4, zeta5 +@nb.njit("f8()", cache=True) def gamma_0(): + """ + Computes the first coefficient of the QCD gamma function. + + Implements Eq. (15) of :cite:`Vermaseren:1997fq`. + + Returns + ------- + gamma_0 : float + first coefficient of the QCD gamma function :math:`\\gamma_{m,0}^{n_f}` + """ return 4.0 -def gamma_1(nf): +@nb.njit("f8(u1)", cache=True) +def gamma_1(nf: int): + """ + Computes the second coefficient of the QCD gamma function. + + Implements Eq. (15) of :cite:`Vermaseren:1997fq`. + + Parameters + ---------- + nf : int + number of active flavors + + Returns + ------- + gamma_1 : float + second coefficient of the QCD gamma function :math:`\\gamma_{m,1}^{n_f}` + """ return 202.0 / 3.0 - 20.0 / 9.0 * nf -def gamma_2(nf): +@nb.njit("f8(u1)", cache=True) +def gamma_2(nf: int): + """ + Computes the third coefficient of the QCD gamma function. + + Implements Eq. (15) of :cite:`Vermaseren:1997fq`. + + Parameters + ---------- + nf : int + number of active flavors + + Returns + ------- + gamma_2 : float + third coefficient of the QCD gamma function :math:`\\gamma_{m,2}^{n_f}` + """ return 1249.0 - (2216.0 / 27.0 + 160.0 / 3.0 * zeta3) * nf - 140.0 / 81.0 * nf ** 2 -def gamma_3(nf): - zeta5 = scipy.special.zeta(5) +@nb.njit("f8(u1)", cache=True) +def gamma_3(nf: int): + """ + Computes the fourth coefficient of the QCD gamma function. + + Implements Eq. (15) of :cite:`Vermaseren:1997fq`. + + Parameters + ---------- + nf : int + number of active flavors + + Returns + ------- + gamma_3 : float + fourth coefficient of the QCD gamma function :math:`\\gamma_{m,3}^{n_f}` + """ return ( 4603055.0 / 162.0 + 135680.0 * zeta3 / 28.0 @@ -40,7 +97,8 @@ def gamma_3(nf): ) -def gamma(order, nf): +@nb.njit("f8(u1,u1)", cache=True) +def gamma(order, nf: int): """ Compute the value of a gamma coefficient diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py index d86001323..8d1f908db 100644 --- a/src/eko/msbar_masses.py +++ b/src/eko/msbar_masses.py @@ -2,6 +2,7 @@ r""" This module contains the |RGE| for the |MSbar| masses """ +import numba as nb import numpy as np from scipy import integrate, optimize @@ -66,7 +67,7 @@ def integrand(a, b_vec, g_vec): val, _ = res[:2] return np.exp(val) - +@nb.njit("f8(f8,f8,u1,u1)", cache=True) def msbar_ker_expanded(a0, a1, order, nf): r""" Expanded |MSbar| |RGE| kernel diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index 859540b41..a3e0f7459 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -13,7 +13,7 @@ import scipy from . import constants, thresholds -from .beta import beta +from .beta import beta, b logger = logging.getLogger(__name__) @@ -67,14 +67,12 @@ def as_expanded(order, as_ref, nf, scale_from, scale_to): res = as_LO # NLO if order >= 1: - beta1 = beta(1, nf) - b1 = beta1 / beta0 + b1 = b(1, nf) as_NLO = as_LO * (1 - b1 * as_LO * np.log(den)) res = as_NLO # NNLO if order >= 2: - beta2 = beta(2, nf) - b2 = beta2 / beta0 + b2 = b(2, nf) res = as_LO * ( 1.0 + as_LO * (as_LO - as_ref) * (b2 - b1 ** 2) @@ -83,7 +81,7 @@ def as_expanded(order, as_ref, nf, scale_from, scale_to): # N3LO expansion is taken from Luca Rottoli and it's simpler # than the APFEL expansion which is more accurate if order >= 3: - b3 = beta(3, nf) / beta0 + b3 = b(3, nf) log_fact = np.log(as_LO / (1 + beta0 * lmu * as_LO)) res += ( as_LO ** 4 @@ -297,19 +295,13 @@ def compute_exact(self, as_ref, nf, scale_from, scale_to): b_vec = [1] # NLO if self.order >= 1: - beta1 = beta(1, nf) - b1 = beta1 / beta0 - b_vec.append(b1) + b_vec.append(b(1, nf)) # NNLO if self.order >= 2: - beta2 = beta(2, nf) - b2 = beta2 / beta0 - b_vec.append(b2) + b_vec.append(b(2, nf)) # N3LO if self.order >= 3: - beta3 = beta(3, nf) - b3 = beta3 / beta0 - b_vec.append(b3) + b_vec.append(b(3, nf)) # integration kernel def rge(_t, a, b_vec): return -(a ** 2) * np.sum([a ** k * b for k, b in enumerate(b_vec)]) From 1a7f5036785e9183c6e4c28389cd4d9152254868 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Tue, 8 Feb 2022 09:28:30 +0100 Subject: [PATCH 12/24] Split evolve msbar and solve msbar. Fix Sc todo --- src/eko/msbar_masses.py | 165 +++++++++++++++-------------- tests/benchmark_msbar_evolution.py | 2 +- tests/test_msbar_masses.py | 8 +- 3 files changed, 93 insertions(+), 82 deletions(-) diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py index 8d1f908db..5bccb02bb 100644 --- a/src/eko/msbar_masses.py +++ b/src/eko/msbar_masses.py @@ -9,7 +9,7 @@ from .beta import b, beta from .evolution_operator.flavors import quark_names from .gamma import gamma -from .strong_coupling import StrongCoupling, strong_coupling_mod_ev +from .strong_coupling import StrongCoupling def msbar_ker_exact(a0, a1, order, nf): @@ -34,7 +34,7 @@ def msbar_ker_exact(a0, a1, order, nf): .. math:: k_{exact} = e^{- \int_{a_s(\mu_{h,0}^2)}^{a_s(\mu^2)} \gamma_m(a_s) / \beta(a_s) da_s} - """ # pylint:disable=line-too-long + """ # pylint:disable=line-too-long b_vec = [beta(0, nf)] g_vec = [gamma(0, nf)] if order >= 1: @@ -67,6 +67,7 @@ def integrand(a, b_vec, g_vec): val, _ = res[:2] return np.exp(val) + @nb.njit("f8(f8,f8,u1,u1)", cache=True) def msbar_ker_expanded(a0, a1, order, nf): r""" @@ -171,19 +172,16 @@ def msbar_ker_dispatcher(q2_to, q2m_ref, strong_coupling, fact_to_ren, nf): return msbar_ker_exact(a0, a1, order, nf) -def evolve_msbar_mass( +def solve_msbar_mass( m2_ref, q2m_ref, strong_coupling, - nf_ref=None, - fact_to_ren=1.0, - q2_to=None, + nf_ref, + fact_to_ren, ): r""" - Compute the |MSbar| mass. - If the final scale is not given it solves the equation :math:`m_{\overline{MS}}(m) = m` + Compute the |MSbar| mass solving the equation :math:`m_{\overline{MS}}(m) = m` for a fixed number of nf - Else perform the mass evolution. Parameters ---------- @@ -191,13 +189,56 @@ def evolve_msbar_mass( squared initial mass reference q2m_ref: float squared initial scale - nf_ref: int, optional (not used when q2_to is given) + strong_coupling: eko.strong_coupling.StrongCoupling + Instance of :class:`~eko.strong_coupling.StrongCoupling` able to generate a_s for + any q + nf_ref: int number of active flavours at the scale q2m_ref, where the solution is searched fact_to_ren: float :math:`\mu_F^2/\muR^2` + + Returns + ------- + m2 : float + :math:`m_{\overline{MS}}(\mu_2)^2` + """ + + def rge(m2, q2m_ref, strong_coupling, fact_to_ren, nf_ref): + return ( + m2_ref + * msbar_ker_dispatcher(m2, q2m_ref, strong_coupling, fact_to_ren, nf_ref) + ** 2 + - m2 + ) + + msbar_mass = optimize.fsolve( + rge, q2m_ref, args=(q2m_ref, strong_coupling, fact_to_ren, nf_ref) + ) + return float(msbar_mass) + + +def evolve_msbar_mass( + m2_ref, + q2m_ref, + strong_coupling, + fact_to_ren, + q2_to, +): + r""" + Perform the |MSbar| mass evolution up to given scale. + It allows for different number of active flavors. + + Parameters + ---------- + m2_ref: float + squared initial mass reference + q2m_ref: float + squared initial scale strong_coupling: eko.strong_coupling.StrongCoupling Instance of :class:`~eko.strong_coupling.StrongCoupling` able to generate a_s for any q + fact_to_ren: float + :math:`\mu_F^2/\muR^2` q2_to: float, optional scale at which the mass is computed @@ -206,58 +247,38 @@ def evolve_msbar_mass( m2 : float :math:`m_{\overline{MS}}(\mu_2)^2` """ - - # TODO: split function: 1 function, 1 job - # TODO: this should resolve the argument priority - if q2_to is None: - - def rge(m2, q2m_ref, strong_coupling, fact_to_ren, nf_ref): - return ( - m2_ref - * msbar_ker_dispatcher( - m2, q2m_ref, strong_coupling, fact_to_ren, nf_ref - ) - ** 2 - - m2 - ) - - msbar_mass = optimize.fsolve( - rge, q2m_ref, args=(q2m_ref, strong_coupling, fact_to_ren, nf_ref) + # evolution might involve different number of active flavors. + # Find out the evolution path (always sorted) + q2_low, q2_high = sorted([q2m_ref, q2_to]) + area_walls = np.array(strong_coupling.thresholds.area_walls) + path = np.concatenate( + ( + [q2_low], + area_walls[(q2_low < area_walls) & (area_walls < q2_high)], + [q2_high], ) - return float(msbar_mass) - else: - - # evolution might involve different number of active flavors. - # Find out the evolution path (always sorted) - q2_low, q2_high = sorted([q2m_ref, q2_to]) - area_walls = np.array(strong_coupling.thresholds.area_walls) - path = np.concatenate( - ( - [q2_low], - area_walls[(q2_low < area_walls) & (area_walls < q2_high)], - [q2_high], - ) + ) + nf_init = len(area_walls[q2_low >= area_walls]) + 2 + nf_final = len(area_walls[q2_high > area_walls]) + 2 + + ev_mass = 1.0 + for i, n in enumerate(np.arange(nf_init, nf_final + 1)): + q2_init = path[i] + q2_final = path[i + 1] + # if you are going backward + # need to reverse the evolution in each path segment + if q2m_ref > q2_to: + q2_init, q2_final = q2_final, q2_init + ev_mass *= msbar_ker_dispatcher( + q2_final, q2_init, strong_coupling, fact_to_ren, n ) - nf_init = len(area_walls[q2_low >= area_walls]) + 2 - nf_final = len(area_walls[q2_high > area_walls]) + 2 - - ev_mass = 1.0 - for i, n in enumerate(np.arange(nf_init, nf_final + 1)): - q2_init = path[i] - q2_final = path[i + 1] - # if you are going backward - # need to reverse the evolution in each path segment - if q2m_ref > q2_to: - q2_init, q2_final = q2_final, q2_init - ev_mass *= msbar_ker_dispatcher( - q2_final, q2_init, strong_coupling, fact_to_ren, n - ) - return m2_ref * ev_mass ** 2 + return m2_ref * ev_mass ** 2 def compute_msbar_mass(theory_card): r""" Compute the |MSbar| masses solving the equation :math:`m_{\bar{MS}}(\mu) = \mu` + for each heavy quark and consistent boundary contitions. Parameters ---------- @@ -276,18 +297,8 @@ def compute_msbar_mass(theory_card): masses = np.concatenate((np.zeros(nfa_ref - 3), np.full(6 - nfa_ref, np.inf))) fact_to_ren = theory_card["fact_to_ren_scale_ratio"] ** 2 - # TODO: why is this different than calling - # `StrongCoupling.from_dict(theory_card, masses=thr_masses)` def sc(thr_masses): - return StrongCoupling( - theory_card["alphas"], - q2_ref, - thr_masses, - thresholds_ratios=[1, 1, 1], - order=theory_card["PTO"], - method=strong_coupling_mod_ev(theory_card["ModEv"]), - nf_ref=nfa_ref, - ) + return StrongCoupling.from_dict(theory_card, masses=thr_masses) # First you need to look for the thr around the given as_ref heavy_quarks = quark_names[3:] @@ -350,19 +361,19 @@ def sc(thr_masses): m2_ref = evolve_msbar_mass( m2_ref, q2m_ref, - strong_coupling=sc(masses), - fact_to_ren=fact_to_ren, - q2_to=q2_to, + sc(masses), + fact_to_ren, + q2_to, ) q2m_ref = q2_to # now solve the RGE - masses[q_idx] = evolve_msbar_mass( + masses[q_idx] = solve_msbar_mass( m2_ref, q2m_ref, - strong_coupling=sc(masses), - nf_ref=nf_target, - fact_to_ren=fact_to_ren, + sc(masses), + nf_target, + fact_to_ren, ) # Check the msbar ordering @@ -373,11 +384,11 @@ def sc(thr_masses): m2_test = evolve_msbar_mass( m2_ref, q2m_ref, - strong_coupling=sc(masses), - fact_to_ren=fact_to_ren, - q2_to=m2_msbar, + sc(masses), + fact_to_ren, + m2_msbar, ) - if m2_msbar > m2_test: + if m2_msbar >= m2_test: raise ValueError( "The MSBAR masses do not preserve the correct ordering,\ check the initial reference values" diff --git a/tests/benchmark_msbar_evolution.py b/tests/benchmark_msbar_evolution.py index 80c95cb60..544f27cec 100644 --- a/tests/benchmark_msbar_evolution.py +++ b/tests/benchmark_msbar_evolution.py @@ -103,7 +103,7 @@ def benchmark_APFEL_msbar_evolution(self): apfel.CleanUp() apfel.SetTheory("QCD") apfel.SetPerturbativeOrder(order) - apfel.SetAlphaEvolution("exact") + apfel.SetAlphaEvolution(method) apfel.SetAlphaQCDRef(alphas_ref, np.sqrt(scale_ref)) apfel.SetVFNS() apfel.SetMSbarMasses(*np.sqrt(m2)) diff --git a/tests/test_msbar_masses.py b/tests/test_msbar_masses.py index 002a04647..31731f593 100644 --- a/tests/test_msbar_masses.py +++ b/tests/test_msbar_masses.py @@ -96,10 +96,10 @@ def test_errors(self): with pytest.raises(ValueError, match="do not preserve the correct ordering"): theory_dict.update( dict( - mc=1.0, - mb=1.00001, - Qmc=1.0, - Qmb=1.00001, + mc=1.009, + mb=1.01, + Qmc=1.012, + Qmb=1.01, ) ) compute_msbar_mass(theory_dict) From 82e2c2b9fddc0389f5399da89e3dbe39b140dcd0 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Tue, 8 Feb 2022 09:29:19 +0100 Subject: [PATCH 13/24] Fix pylint import order --- src/eko/output.py | 10 +++++----- src/ekomark/benchmark/external/LHA_utils.py | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/eko/output.py b/src/eko/output.py index 45c2f387c..f6b7aa963 100644 --- a/src/eko/output.py +++ b/src/eko/output.py @@ -9,9 +9,9 @@ import tempfile import warnings -import lz4.frame import numpy as np import yaml +from lz4 import frame from . import basis_rotation as br from . import interpolation, version @@ -308,7 +308,7 @@ def get_raw(self, binarize=True, skip_q2_grid=False): out["Q2grid"][q2][k] = float(v) continue if binarize: - out["Q2grid"][q2][k] = lz4.frame.compress(v.tobytes()) + out["Q2grid"][q2][k] = frame.compress(v.tobytes()) else: out["Q2grid"][q2][k] = v.tolist() else: @@ -393,7 +393,7 @@ def dump_tar(self, tarname): stream = io.BytesIO() np.save(stream, operator) stream.seek(0) - with lz4.frame.open( + with frame.open( (tmpdir / kind).with_suffix(".npy.lz4"), "wb" ) as fo: fo.write(stream.read()) @@ -436,7 +436,7 @@ def load_yaml(cls, stream, skip_q2_grid=False): elif isinstance(v, list): v = np.array(v) elif isinstance(v, bytes): - v = np.frombuffer(lz4.frame.decompress(v)) + v = np.frombuffer(frame.decompress(v)) v = v.reshape(len_tpids, len_tgrid, len_ipids, len_igrid) op[k] = v return cls(obj) @@ -500,7 +500,7 @@ def load_tar(cls, tarname): grids = {} for fp in innerdir.glob("*.npy.lz4"): - with lz4.frame.open(fp, "rb") as fd: + with frame.open(fp, "rb") as fd: stream = io.BytesIO(fd.read()) stream.seek(0) grids[pathlib.Path(fp.stem).stem] = np.load(stream) diff --git a/src/ekomark/benchmark/external/LHA_utils.py b/src/ekomark/benchmark/external/LHA_utils.py index 89589b831..d91624180 100644 --- a/src/ekomark/benchmark/external/LHA_utils.py +++ b/src/ekomark/benchmark/external/LHA_utils.py @@ -4,9 +4,9 @@ """ import pathlib -import matplotlib.pyplot as plt import numpy as np import yaml +from matplotlib import pyplot as plt from matplotlib.backends.backend_pdf import PdfPages from eko import basis_rotation as br From c148c1a54df67d51eb2a544b90b12508211fd799 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Wed, 9 Feb 2022 19:29:53 +0100 Subject: [PATCH 14/24] Adding MSBar metching relations with kthr!=1 --- benchmarks/apfel_bench.py | 3 +- doc/source/refs.bib | 14 ++++ src/eko/msbar_masses.py | 115 +++++++++++++++++++++------ src/eko/strong_coupling.py | 45 +++++++---- tests/benchmark_msbar_evolution.py | 122 ++++++++++++++++++++++++++++- tests/test_msbar_masses.py | 22 +++--- 6 files changed, 271 insertions(+), 50 deletions(-) diff --git a/benchmarks/apfel_bench.py b/benchmarks/apfel_bench.py index 928b9c708..db69ac93c 100644 --- a/benchmarks/apfel_bench.py +++ b/benchmarks/apfel_bench.py @@ -95,7 +95,8 @@ def benchmark_msbar(self, pto): when passing kthr != 1 both apfel and eko use ``kThr * msbar``, as thr scale, where ``msbar`` is the usual ms_bar solution. However apfel and eko mange the alpha_s thr differently - (apfel uses the given mass parameters as thr), so the + (apfel uses the given mass parameters as thr multipied by the kthr), + while in eko only the already computed thr matters, so the benchmark is not a proper comparison with this option allowed. """ th = self.vfns_theory.copy() diff --git a/doc/source/refs.bib b/doc/source/refs.bib index e9fdcec6e..a46b8394f 100644 --- a/doc/source/refs.bib +++ b/doc/source/refs.bib @@ -431,4 +431,18 @@ @article{Vermaseren:1997fq volume = "405", pages = "327--333", year = "1997" +} + +@article{Liu:2015fxa, + author = "Liu, Tao and Steinhauser, Matthias", + title = "{Decoupling of heavy quarks at four loops and effective Higgs-fermion coupling}", + eprint = "1502.04719", + archivePrefix = "arXiv", + primaryClass = "hep-ph", + reportNumber = "TTP15-05", + doi = "10.1016/j.physletb.2015.05.023", + journal = "Phys. Lett. B", + volume = "746", + pages = "330--334", + year = "2015" } \ No newline at end of file diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py index 5bccb02bb..a57d6e0e3 100644 --- a/src/eko/msbar_masses.py +++ b/src/eko/msbar_masses.py @@ -9,7 +9,7 @@ from .beta import b, beta from .evolution_operator.flavors import quark_names from .gamma import gamma -from .strong_coupling import StrongCoupling +from .strong_coupling import StrongCoupling, invert_matching_coeffs def msbar_ker_exact(a0, a1, order, nf): @@ -250,6 +250,7 @@ def evolve_msbar_mass( # evolution might involve different number of active flavors. # Find out the evolution path (always sorted) q2_low, q2_high = sorted([q2m_ref, q2_to]) + area_walls = np.array(strong_coupling.thresholds.area_walls) path = np.concatenate( ( @@ -262,17 +263,83 @@ def evolve_msbar_mass( nf_final = len(area_walls[q2_high > area_walls]) + 2 ev_mass = 1.0 - for i, n in enumerate(np.arange(nf_init, nf_final + 1)): + is_downward_path = bool(q2m_ref > q2_to) + for i, nf in enumerate(np.arange(nf_init, nf_final + 1)): q2_init = path[i] q2_final = path[i + 1] # if you are going backward # need to reverse the evolution in each path segment - if q2m_ref > q2_to: + if is_downward_path: + m_coeffs = compute_matching_coeffs_down(nf - 1) q2_init, q2_final = q2_final, q2_init - ev_mass *= msbar_ker_dispatcher( - q2_final, q2_init, strong_coupling, fact_to_ren, n + shift = 4 + else: + m_coeffs = compute_matching_coeffs_up(nf) + shift = 3 + fact = 1.0 + # shift + for pto in range(1, strong_coupling.order + 1): + for l in range(pto + 1): + as_thr = strong_coupling.a_s(q2_final / fact_to_ren, q2_final) + # TODO: do we need to add np.log(fac_to_ren) here ??? + L = np.log(strong_coupling.thresholds.thresholds_ratios[pto - shift]) + fact += as_thr ** pto * L ** l * m_coeffs[pto, l] + ev_mass *= ( + fact + * msbar_ker_dispatcher(q2_final, q2_init, strong_coupling, fact_to_ren, nf) + ** 2 ) - return m2_ref * ev_mass ** 2 + return m2_ref * ev_mass + + +def compute_matching_coeffs_up(nf): + """ + |MSbar| matching coefficients :cite:`Liu:2015fxa` at threshold + when moving to a regime with *more* flavors. + + Note :cite:`Liu:2015fxa` (eq 5.) reports the backward relation, + multiplied by a factor of 4 (and 4^2 ...) + + Parameters + ---------- + nf: + number of active flavors in the lower patch + + Returns + ------- + matching_coeffs_up: + downward matching coefficient matrix + """ + matching_coeffs_up = np.zeros((4, 4)) + matching_coeffs_up[2, 0] = -89.0 / 27.0 + matching_coeffs_up[2, 1] = 20.0 / 9.0 + matching_coeffs_up[2, 2] = -4.0 / 3.0 + + matching_coeffs_up[3, 0] = -118.248 - 1.58257 * nf + matching_coeffs_up[3, 1] = 71.7887 + 7.85185 * nf + matching_coeffs_up[3, 2] = -700.0 / 27.0 + matching_coeffs_up[3, 3] = -232.0 / 27.0 + 16.0 / 27.0 * nf + + return matching_coeffs_up + + +def compute_matching_coeffs_down(nf): + """ + |MSbar| matching coefficients :cite:`Liu:2015fxa` (eq 5) at threshold + when moving to a regime with *less* flavors. + + Parameters + ---------- + nf: + number of active flavors in the lower patch + + Returns + ------- + matching_coeffs_down: + downward matching coefficient matrix + """ + _c = compute_matching_coeffs_up(nf) + return invert_matching_coeffs(_c) def compute_msbar_mass(theory_card): @@ -376,21 +443,25 @@ def sc(thr_masses): fact_to_ren, ) + # TODO: this test seems to be quite hard to contradict + # what if we check just the mass ordering? # Check the msbar ordering - for m2_msbar, hq in zip(masses[:-1], quark_names[4:]): - q2m_ref = np.power(theory_card[f"Qm{hq}"], 2) - m2_ref = np.power(theory_card[f"m{hq}"], 2) - # check that m_msbar_hq < msbar_hq+1 (m_msbar_hq) - m2_test = evolve_msbar_mass( - m2_ref, - q2m_ref, - sc(masses), - fact_to_ren, - m2_msbar, - ) - if m2_msbar >= m2_test: - raise ValueError( - "The MSBAR masses do not preserve the correct ordering,\ - check the initial reference values" - ) + if not (masses == np.sort(masses)).all(): + raise ValueError("masses need to be sorted") + # for m2_msbar, hq in zip(masses[:-1], quark_names[4:]): + # q2m_ref = np.power(theory_card[f"Qm{hq}"], 2) + # m2_ref = np.power(theory_card[f"m{hq}"], 2) + # # check that m_msbar_hq < msbar_hq+1 (m_msbar_hq) + # m2_test = evolve_msbar_mass( + # m2_ref, + # q2m_ref, + # sc(masses), + # fact_to_ren, + # m2_msbar, + # ) + # if m2_msbar >= m2_test: + # raise ValueError( + # "The MSBAR masses do not preserve the correct ordering,\ + # check the initial reference values" + # ) return masses diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index a3e0f7459..1b6120a0c 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -13,7 +13,7 @@ import scipy from . import constants, thresholds -from .beta import beta, b +from .beta import b, beta logger = logging.getLogger(__name__) @@ -465,6 +465,24 @@ def compute_matching_coeffs_down(mass_scheme, nf): Matching coefficients :cite:`Schroder:2005hy,Chetyrkin:2005ia` at threshold when moving to a regime with *less* flavors. + Parameters + ---------- + mass_scheme: + Heavy quark mass scheme: "POLE" or "MSBAR" + nf: + number of active flavors in the lower patch + + Returns + ------- + matching_coeffs_down: + downward matching coefficient matrix + """ + _c = compute_matching_coeffs_up(mass_scheme, nf) + return invert_matching_coeffs(_c) + + +def invert_matching_coeffs(c_up): + """ This is the perturbative inverse of :data:`matching_coeffs_up` and has been obtained via .. code-block:: Mathematica @@ -483,24 +501,21 @@ def compute_matching_coeffs_down(mass_scheme, nf): Parameters ---------- - mass_scheme: - Heavy quark mass scheme: "POLE" or "MSBAR" - nf: - number of active flavors in the lower patch + c_up : numpy.ndarray + forward matching coefficient matrix Returns ------- matching_coeffs_down: downward matching coefficient matrix """ - _c = compute_matching_coeffs_up(mass_scheme, nf) - matching_coeffs_down = np.zeros_like(_c) - matching_coeffs_down[1, 1] = -_c[1, 1] - matching_coeffs_down[2, 0] = -_c[2, 0] - matching_coeffs_down[2, 1] = -_c[2, 1] - matching_coeffs_down[2, 2] = 2.0 * _c[1, 1] ** 2 - _c[2, 2] - matching_coeffs_down[3, 0] = -_c[3, 0] - matching_coeffs_down[3, 1] = 5 * _c[1, 1] * _c[2, 0] - _c[3, 1] - matching_coeffs_down[3, 2] = 5 * _c[1, 1] * _c[2, 1] - _c[3, 2] - matching_coeffs_down[3, 3] = -5 * _c[1, 1] ** 3 + 5 * _c[1, 1] * _c[2, 2] - _c[3, 3] + matching_coeffs_down = np.zeros_like(c_up) + matching_coeffs_down[1, 1] = - c_up[1, 1] + matching_coeffs_down[2, 0] = - c_up[2, 0] + matching_coeffs_down[2, 1] = - c_up[2, 1] + matching_coeffs_down[2, 2] = 2.0 * c_up[1, 1] ** 2 - c_up[2, 2] + matching_coeffs_down[3, 0] = - c_up[3, 0] + matching_coeffs_down[3, 1] = 5 *c_up[1, 1] * c_up[2, 0] - c_up[3, 1] + matching_coeffs_down[3, 2] = 5 * c_up[1, 1] * c_up[2, 1] - c_up[3, 2] + matching_coeffs_down[3, 3] = -5 * c_up[1, 1] ** 3 + 5 * c_up[1, 1] * c_up[2, 2] - c_up[3, 3] return matching_coeffs_down diff --git a/tests/benchmark_msbar_evolution.py b/tests/benchmark_msbar_evolution.py index 544f27cec..c837b5711 100644 --- a/tests/benchmark_msbar_evolution.py +++ b/tests/benchmark_msbar_evolution.py @@ -2,7 +2,7 @@ """This module benchmarks MSbar mass evolution against APFEL.""" import numpy as np -from eko.msbar_masses import evolve_msbar_mass +from eko.msbar_masses import compute_msbar_mass, evolve_msbar_mass, compute_msbar_mass from eko.strong_coupling import StrongCoupling # try to load APFEL - if not available, we'll use the cached values @@ -13,6 +13,27 @@ except ImportError: use_APFEL = False +theory_dict = { + "alphas": 0.1180, + "Qref": 91, + "nfref": 5, + "MaxNfPdf": 6, + "MaxNfAs": 6, + "Q0": 1, + "fact_to_ren_scale_ratio": 1.0, + "mc": 1.5, + "mb": 4.1, + "mt": 175.0, + "kcThr": 1.0, + "kbThr": 1.0, + "ktThr": 1.0, + "HQ": "MSBAR", + "Qmc": 18, + "Qmb": 20, + "Qmt": 175.0, + "PTO": 2, + "ModEv": "EXA", +} class BenchmarkMSbar: def benchmark_APFEL_msbar_evolution(self): @@ -125,3 +146,102 @@ def benchmark_APFEL_msbar_evolution(self): np.testing.assert_allclose( apfel_vals, np.sqrt(np.array(my_vals)), rtol=5e-3 ) + + def benchmark_APFEL_msbar_solution(self): + apfel_vals_dict = { + "EXA": { + 0: np.array( + [1.9855, 4.8062, 175.0000], + ), + 1: np.array( + [2.1308, 4.9656, 175.0000], + ), + 2: np.array([2.1566, 4.9841, 175.0000]), + }, + } + # collect my values + for order in [0, 1, 2]: + theory_dict["PTO"] = order + my_masses = compute_msbar_mass(theory_dict) + # get APFEL numbers - if available else use cache + apfel_vals = apfel_vals_dict[theory_dict["ModEv"]][order] + if use_APFEL: + # run apfel + apfel.CleanUp() + apfel.SetTheory("QCD") + apfel.SetPerturbativeOrder(order) + apfel.SetAlphaEvolution("exact") + apfel.SetAlphaQCDRef(theory_dict["alphas"], theory_dict["Qref"]) + apfel.SetVFNS() + apfel.SetMSbarMasses( + theory_dict["mc"], theory_dict["mb"], theory_dict["mt"] + ) + apfel.SetMassScaleReference( + theory_dict["Qmc"], theory_dict["Qmb"], theory_dict["Qmt"] + ) + apfel.SetRenFacRatio(theory_dict["fact_to_ren_scale_ratio"]) + apfel.InitializeAPFEL() + apfel.EnableWelcomeMessage(1) + # check myself to APFEL + np.testing.assert_allclose( + apfel_vals, np.sqrt(np.array(my_masses)), rtol=5e-5 + ) + + # With this test you can see that in Apfel + # the solution value of mb is affected by "kbThr", while + # in EKO this doesn't happend, since mb is searched + # with an Nf=5 larger range. + + # def benchmark_APFEL_msbar_solution_kthr(self): + # theory_dict.update({ + # "mc": 1.5, + # "mb": 4.1, + # "mt": 175.0, + # "kcThr": 1.2, + # "kbThr": 1.8, + # "ktThr": 1.0, + # "Qmc": 18, + # "Qmb": 20, + # "Qmt": 175.0, + # }) + # apfel_vals_dict = { + # "EXA": { + # 0: np.array( + # [1.9891,4.5102,175.0000], + # ), + # 1: np.array( + # [2.1316, 4.6047, 175.0000], + # ), + # 2: np.array([2.1574, 4.6153, 175.0000]), + # }, + # } + # # collect my values + # for order in [0, 1, 2]: + # theory_dict["PTO"] = order + # my_masses = compute_msbar_mass(theory_dict) + # # get APFEL numbers - if available else use cache + # apfel_vals = apfel_vals_dict[theory_dict["ModEv"]][order] + # if use_APFEL: + # # run apfel + # apfel.CleanUp() + # apfel.SetTheory("QCD") + # apfel.SetPerturbativeOrder(order) + # apfel.SetAlphaEvolution("exact") + # apfel.SetAlphaQCDRef(theory_dict["alphas"], theory_dict["Qref"]) + # apfel.SetVFNS() + # apfel.SetMSbarMasses( + # theory_dict["mc"], theory_dict["mb"], theory_dict["mt"] + # ) + # apfel.SetMassScaleReference( + # theory_dict["Qmc"], theory_dict["Qmb"], theory_dict["Qmt"] + # ) + # apfel.SetMassMatchingScales( + # theory_dict["kcThr"], theory_dict["kbThr"], theory_dict["ktThr"] + # ) + # apfel.SetRenFacRatio(theory_dict["fact_to_ren_scale_ratio"]) + # apfel.InitializeAPFEL() + # apfel.EnableWelcomeMessage(1) + # #check myself to APFEL + # np.testing.assert_allclose( + # apfel_vals, np.sqrt(np.array(my_masses)), rtol=5e-5 + # ) diff --git a/tests/test_msbar_masses.py b/tests/test_msbar_masses.py index 31731f593..a0093a884 100644 --- a/tests/test_msbar_masses.py +++ b/tests/test_msbar_masses.py @@ -56,7 +56,7 @@ def test_compute_msbar_mass(self): q2_to=m2_computed[nf - 3], ) ) - np.testing.assert_allclose(m2_computed, m2_test, rtol=5e-3) + np.testing.assert_allclose(m2_computed, m2_test, rtol=9e-3 if order==3 else 5e-3) def test_compute_msbar_mass_VFNS(self): # test the solution now with some initial contition @@ -93,16 +93,16 @@ def test_compute_msbar_mass_VFNS(self): def test_errors(self): # test mass ordering - with pytest.raises(ValueError, match="do not preserve the correct ordering"): - theory_dict.update( - dict( - mc=1.009, - mb=1.01, - Qmc=1.012, - Qmb=1.01, - ) - ) - compute_msbar_mass(theory_dict) + # with pytest.raises(ValueError, match="do not preserve the correct ordering"): + # theory_dict.update( + # dict( + # mc=1.009, + # mb=1.01, + # Qmc=1.012, + # Qmb=1.01, + # ) + # ) + # compute_msbar_mass(theory_dict) with pytest.raises(ValueError, match="masses need to be sorted"): theory_dict.update( dict( From 924078c23dcd68bdf1a7447b7cef2f3c838b6287 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Thu, 10 Feb 2022 09:34:37 +0100 Subject: [PATCH 15/24] Improve N3LO alpha_s expanded --- doc/source/theory/pQCD.rst | 12 ++++++------ src/eko/strong_coupling.py | 32 +++++++++++++++++--------------- tests/test_strong_coupling.py | 2 +- 3 files changed, 24 insertions(+), 22 deletions(-) diff --git a/doc/source/theory/pQCD.rst b/doc/source/theory/pQCD.rst index e856f37a5..967f7018e 100644 --- a/doc/source/theory/pQCD.rst +++ b/doc/source/theory/pQCD.rst @@ -27,19 +27,19 @@ We implement two different strategies to solve the |RGE|: a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2) &= a^{\text{LO}}_s(\mu_R^2)-b_1 \left[a^{\text{LO}}_s(\mu_R^2)\right]^2 \ln\left(1+a_s(\mu_0^2) \beta_0 L_{\mu}\right) \\ a^{\text{NNLO}}_{s,\text{exp}}(\mu_R^2) &= a^{\text{LO}}_s(\mu_R^2)\left[1 + a^{\text{LO}}_s(\mu_R^2)\left(a^{\text{LO}}_s(\mu_R^2) - a_s(\mu_0^2)\right)(b_2 - b_1^2) \right.\\ & \hspace{60pt} \left. + a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2) b_1 \ln\left(a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2)/a_s(\mu_0^2)\right)\right] \\ - a^{\text{N3LO}}_{s,\text{exp}}(\mu_R^2) &= a^{\text{NNLO}}_s(\mu_R^2) + \frac{a^{\text{LO}}_s(\mu_R^2)^4}{2 b_0^3 \left(1 + b_0 L_{\mu} a^{\text{LO}}_s(\mu_R^2)\right)^4} \left\{ \right. \\ - & -2 b_1^3 L_{\text{LO}}^3 + 5 b_1^3 L_{\alpha_s}^2 + 2 b_1^3 L_{\alpha_s}^3 + b_1^3 L_{\text{LO}}^2 \left(5 + 6 L_{\alpha_s} \right) \\ - & + 2 b_0 b_1 L_{\alpha_s} \left[ b_2 + 2 \left(b_1^2 - b_0 b_2 \right) L_{\mu} a^{\text{LO}}_s(\mu_R^2) \right] \\ - & - b_0^2 L_{\mu} a^{\text{LO}}_s(\mu_R^2) \left[ -2 b_1 b_2 + 2 b_0 b_3 + \left( b_1^3 - 2 b_0 b_1 b_2 + b_0^2 b_3 \right) L_{\mu} a^{\text{LO}}_s(\mu_R^2) \right] \\ - & - 2 b_1 L_{\text{LO}} \left[ 5 b_1^2 L_{\alpha_s} + 3 b_1^2 L_{\alpha_s}^2 + b_0 \left[b_2 + 2 \left(b_1^2 - b_0 b_2\right) L_{\mu} a^{\text{LO}}_s(\mu_R^2)\right] \right ] \\ + a^{\text{N3LO}}_{s,\text{exp}}(\mu_R^2) &= a^{\text{NNLO}}_s(\mu_R^2) + \frac{a^{\text{LO}}_s(\mu_R^2)^4}{2 b_0^3} \left\{ \right. \\ + & -2 b_1^3 L_{0}^3 + 5 b_1^3 L_{\text{LO}}^2 + 2 b_1^3 L_{\text{LO}}^3 + b_1^3 L_{0}^2 \left(5 + 6 L_{\text{LO}} \right) \\ + & + 2 b_0 b_1 L_{\text{LO}} \left[ b_2 + 2 \left(b_1^2 - b_0 b_2 \right) L_{\mu} a_s(\mu_0^2) \right] \\ + & - b_0^2 L_{\mu} a_s(\mu_0^2) \left[ -2 b_1 b_2 + 2 b_0 b_3 + \left( b_1^3 - 2 b_0 b_1 b_2 + b_0^2 b_3 \right) L_{\mu} a_s(\mu_0^2) \right] \\ + & - 2 b_1 L_{0} \left[ 5 b_1^2 L_{\text{LO}} + 3 b_1^2 L_{\text{LO}}^2 + b_0 \left[b_2 + 2 \left(b_1^2 - b_0 b_2\right) L_{\mu} a_s(\mu_0^2)\right] \right ] \\ & \left. \right\} being: .. math :: L_{\mu} &= \ln(\mu_R^2/\mu_0^2) \\ + L_{0} &= \ln(a_s(\mu_0^2)) \\ L_{\text{LO}} &= \ln(a^{\text{LO}}_s(\mu_R^2)) \\ - L_{\alpha_s} &= \ln \left( \frac{a^{\text{LO}}_s(\mu_R^2)}{(1 + b_0 L_{\mu} a^{\text{LO}}_s(\mu_R^2)} \right) \\ When the renormalization scale crosses a flavor threshold matching conditions have to be applied :cite:`Schroder:2005hy,Chetyrkin:2005ia`. diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index 1b6120a0c..a4f6a7a59 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -82,37 +82,37 @@ def as_expanded(order, as_ref, nf, scale_from, scale_to): # than the APFEL expansion which is more accurate if order >= 3: b3 = b(3, nf) - log_fact = np.log(as_LO / (1 + beta0 * lmu * as_LO)) + log_fact = np.log(as_LO) res += ( as_LO ** 4 - / (2 * beta0 ** 3 * (1 + beta0 * lmu * as_LO) ** 4) + / (2 * beta0 ** 3) * ( - -2 * b1 ** 3 * np.log(as_LO) ** 3 + -2 * b1 ** 3 * np.log(as_ref) ** 3 + 5 * b1 ** 3 * log_fact ** 2 + 2 * b1 ** 3 * log_fact ** 3 - + b1 ** 3 * np.log(as_LO) ** 2 * (5 + 6 * log_fact) + + b1 ** 3 * np.log(as_ref) ** 2 * (5 + 6 * log_fact) + 2 * beta0 * b1 * log_fact - * (b2 + 2 * (b1 ** 2 - beta0 * b2) * lmu * as_LO) + * (b2 + 2 * (b1 ** 2 - beta0 * b2) * lmu * as_ref) - beta0 ** 2 * lmu - * as_LO + * as_ref * ( -2 * b1 * b2 + 2 * beta0 * b3 + (b1 ** 3 - 2 * beta0 * b1 * b2 + beta0 ** 2 * b3) * lmu - * as_LO + * as_ref ) - 2 * b1 - * np.log(as_LO) + * np.log(as_ref) * ( 5 * b1 ** 2 * log_fact + 3 * b1 ** 2 * log_fact ** 2 - + beta0 * (b2 + 2 * (b1 ** 2 - beta0 * b2) * lmu * as_LO) + + beta0 * (b2 + 2 * (b1 ** 2 - beta0 * b2) * lmu * as_ref) ) ) ) @@ -510,12 +510,14 @@ def invert_matching_coeffs(c_up): downward matching coefficient matrix """ matching_coeffs_down = np.zeros_like(c_up) - matching_coeffs_down[1, 1] = - c_up[1, 1] - matching_coeffs_down[2, 0] = - c_up[2, 0] - matching_coeffs_down[2, 1] = - c_up[2, 1] + matching_coeffs_down[1, 1] = -c_up[1, 1] + matching_coeffs_down[2, 0] = -c_up[2, 0] + matching_coeffs_down[2, 1] = -c_up[2, 1] matching_coeffs_down[2, 2] = 2.0 * c_up[1, 1] ** 2 - c_up[2, 2] - matching_coeffs_down[3, 0] = - c_up[3, 0] - matching_coeffs_down[3, 1] = 5 *c_up[1, 1] * c_up[2, 0] - c_up[3, 1] + matching_coeffs_down[3, 0] = -c_up[3, 0] + matching_coeffs_down[3, 1] = 5 * c_up[1, 1] * c_up[2, 0] - c_up[3, 1] matching_coeffs_down[3, 2] = 5 * c_up[1, 1] * c_up[2, 1] - c_up[3, 2] - matching_coeffs_down[3, 3] = -5 * c_up[1, 1] ** 3 + 5 * c_up[1, 1] * c_up[2, 2] - c_up[3, 3] + matching_coeffs_down[3, 3] = ( + -5 * c_up[1, 1] ** 3 + 5 * c_up[1, 1] * c_up[2, 2] - c_up[3, 3] + ) return matching_coeffs_down diff --git a/tests/test_strong_coupling.py b/tests/test_strong_coupling.py index 98408bd11..46d6f7d35 100644 --- a/tests/test_strong_coupling.py +++ b/tests/test_strong_coupling.py @@ -174,7 +174,7 @@ def benchmark_expanded_n3lo(self): m2b = 25 m2t = 30625 threshold_list = [m2c, m2b, m2t] - mathematica_val = -0.000101654 + mathematica_val = -0.000169117 # collect my values as_NNLO = StrongCoupling( alphas_ref, From 30b1e161353f5f7d01874406d02489b73c48a811 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Thu, 10 Feb 2022 09:48:24 +0100 Subject: [PATCH 16/24] Fix msbar test with kthr != 1 --- src/eko/msbar_masses.py | 102 +++++++++---------- tests/benchmark_msbar_evolution.py | 152 ++++++++++++++--------------- tests/test_msbar_masses.py | 8 +- 3 files changed, 130 insertions(+), 132 deletions(-) diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py index a57d6e0e3..c4d12d3a4 100644 --- a/src/eko/msbar_masses.py +++ b/src/eko/msbar_masses.py @@ -172,6 +172,56 @@ def msbar_ker_dispatcher(q2_to, q2m_ref, strong_coupling, fact_to_ren, nf): return msbar_ker_exact(a0, a1, order, nf) +def compute_matching_coeffs_up(nf): + """ + |MSbar| matching coefficients :cite:`Liu:2015fxa` at threshold + when moving to a regime with *more* flavors. + + Note :cite:`Liu:2015fxa` (eq 5.) reports the backward relation, + multiplied by a factor of 4 (and 4^2 ...) + + Parameters + ---------- + nf: + number of active flavors in the lower patch + + Returns + ------- + matching_coeffs_up: + downward matching coefficient matrix + """ + matching_coeffs_up = np.zeros((4, 4)) + matching_coeffs_up[2, 0] = -89.0 / 27.0 + matching_coeffs_up[2, 1] = 20.0 / 9.0 + matching_coeffs_up[2, 2] = -4.0 / 3.0 + + matching_coeffs_up[3, 0] = -118.248 - 1.58257 * nf + matching_coeffs_up[3, 1] = 71.7887 + 7.85185 * nf + matching_coeffs_up[3, 2] = -700.0 / 27.0 + matching_coeffs_up[3, 3] = -232.0 / 27.0 + 16.0 / 27.0 * nf + + return matching_coeffs_up + + +def compute_matching_coeffs_down(nf): + """ + |MSbar| matching coefficients :cite:`Liu:2015fxa` (eq 5) at threshold + when moving to a regime with *less* flavors. + + Parameters + ---------- + nf: + number of active flavors in the lower patch + + Returns + ------- + matching_coeffs_down: + downward matching coefficient matrix + """ + _c = compute_matching_coeffs_up(nf) + return invert_matching_coeffs(_c) + + def solve_msbar_mass( m2_ref, q2m_ref, @@ -292,56 +342,6 @@ def evolve_msbar_mass( return m2_ref * ev_mass -def compute_matching_coeffs_up(nf): - """ - |MSbar| matching coefficients :cite:`Liu:2015fxa` at threshold - when moving to a regime with *more* flavors. - - Note :cite:`Liu:2015fxa` (eq 5.) reports the backward relation, - multiplied by a factor of 4 (and 4^2 ...) - - Parameters - ---------- - nf: - number of active flavors in the lower patch - - Returns - ------- - matching_coeffs_up: - downward matching coefficient matrix - """ - matching_coeffs_up = np.zeros((4, 4)) - matching_coeffs_up[2, 0] = -89.0 / 27.0 - matching_coeffs_up[2, 1] = 20.0 / 9.0 - matching_coeffs_up[2, 2] = -4.0 / 3.0 - - matching_coeffs_up[3, 0] = -118.248 - 1.58257 * nf - matching_coeffs_up[3, 1] = 71.7887 + 7.85185 * nf - matching_coeffs_up[3, 2] = -700.0 / 27.0 - matching_coeffs_up[3, 3] = -232.0 / 27.0 + 16.0 / 27.0 * nf - - return matching_coeffs_up - - -def compute_matching_coeffs_down(nf): - """ - |MSbar| matching coefficients :cite:`Liu:2015fxa` (eq 5) at threshold - when moving to a regime with *less* flavors. - - Parameters - ---------- - nf: - number of active flavors in the lower patch - - Returns - ------- - matching_coeffs_down: - downward matching coefficient matrix - """ - _c = compute_matching_coeffs_up(nf) - return invert_matching_coeffs(_c) - - def compute_msbar_mass(theory_card): r""" Compute the |MSbar| masses solving the equation :math:`m_{\bar{MS}}(\mu) = \mu` @@ -447,7 +447,7 @@ def sc(thr_masses): # what if we check just the mass ordering? # Check the msbar ordering if not (masses == np.sort(masses)).all(): - raise ValueError("masses need to be sorted") + raise ValueError("Msbar masses are not to be sorted") # for m2_msbar, hq in zip(masses[:-1], quark_names[4:]): # q2m_ref = np.power(theory_card[f"Qm{hq}"], 2) # m2_ref = np.power(theory_card[f"m{hq}"], 2) diff --git a/tests/benchmark_msbar_evolution.py b/tests/benchmark_msbar_evolution.py index c837b5711..d8021381b 100644 --- a/tests/benchmark_msbar_evolution.py +++ b/tests/benchmark_msbar_evolution.py @@ -14,27 +14,28 @@ use_APFEL = False theory_dict = { - "alphas": 0.1180, - "Qref": 91, - "nfref": 5, - "MaxNfPdf": 6, - "MaxNfAs": 6, - "Q0": 1, - "fact_to_ren_scale_ratio": 1.0, - "mc": 1.5, - "mb": 4.1, - "mt": 175.0, - "kcThr": 1.0, - "kbThr": 1.0, - "ktThr": 1.0, - "HQ": "MSBAR", - "Qmc": 18, - "Qmb": 20, - "Qmt": 175.0, - "PTO": 2, - "ModEv": "EXA", + "alphas": 0.1180, + "Qref": 91, + "nfref": 5, + "MaxNfPdf": 6, + "MaxNfAs": 6, + "Q0": 1, + "fact_to_ren_scale_ratio": 1.0, + "mc": 1.5, + "mb": 4.1, + "mt": 175.0, + "kcThr": 1.0, + "kbThr": 1.0, + "ktThr": 1.0, + "HQ": "MSBAR", + "Qmc": 18, + "Qmb": 20, + "Qmt": 175.0, + "PTO": 2, + "ModEv": "EXA", } + class BenchmarkMSbar: def benchmark_APFEL_msbar_evolution(self): Q2s = np.power([1, 96, 150], 2) @@ -144,7 +145,7 @@ def benchmark_APFEL_msbar_evolution(self): ) # check myself to APFEL np.testing.assert_allclose( - apfel_vals, np.sqrt(np.array(my_vals)), rtol=5e-3 + apfel_vals, np.sqrt(np.array(my_vals)), rtol=2e-3 ) def benchmark_APFEL_msbar_solution(self): @@ -187,61 +188,58 @@ def benchmark_APFEL_msbar_solution(self): apfel_vals, np.sqrt(np.array(my_masses)), rtol=5e-5 ) - # With this test you can see that in Apfel - # the solution value of mb is affected by "kbThr", while - # in EKO this doesn't happend, since mb is searched - # with an Nf=5 larger range. + def benchmark_msbar_solution_kthr(self): + """ + With this test you can see that in EKO + the solution value of mb is not affected by "kbThr", + since mb is searched with an Nf=5 larger range. + While in Apfel this doesn't happend + """ + theory_dict.update( + { + "mc": 1.5, + "mb": 4.1, + "mt": 175.0, + "kcThr": 1.2, + "kbThr": 1.8, + "ktThr": 1.0, + "Qmc": 18, + "Qmb": 20, + "Qmt": 175.0, + "PTO": 0, + } + ) + my_masses_thr = compute_msbar_mass(theory_dict) + apfel_masses_thr = [1.9891, 4.5102, 175.0000] + theory_dict.update( + { + "kcThr": 1.0, + "kbThr": 1.0, + } + ) + my_masses_plain = compute_msbar_mass(theory_dict) + apfel_masses_plain = ([1.9855, 4.8062, 175.0000],) + + # Eko bottom mass is the same + np.testing.assert_allclose(my_masses_thr[1], my_masses_plain[1]) + # Eko charm mass is not the same + np.testing.assert_raises( + AssertionError, + np.testing.assert_allclose, + my_masses_thr[0], + my_masses_plain[0], + ) - # def benchmark_APFEL_msbar_solution_kthr(self): - # theory_dict.update({ - # "mc": 1.5, - # "mb": 4.1, - # "mt": 175.0, - # "kcThr": 1.2, - # "kbThr": 1.8, - # "ktThr": 1.0, - # "Qmc": 18, - # "Qmb": 20, - # "Qmt": 175.0, - # }) - # apfel_vals_dict = { - # "EXA": { - # 0: np.array( - # [1.9891,4.5102,175.0000], - # ), - # 1: np.array( - # [2.1316, 4.6047, 175.0000], - # ), - # 2: np.array([2.1574, 4.6153, 175.0000]), - # }, - # } - # # collect my values - # for order in [0, 1, 2]: - # theory_dict["PTO"] = order - # my_masses = compute_msbar_mass(theory_dict) - # # get APFEL numbers - if available else use cache - # apfel_vals = apfel_vals_dict[theory_dict["ModEv"]][order] - # if use_APFEL: - # # run apfel - # apfel.CleanUp() - # apfel.SetTheory("QCD") - # apfel.SetPerturbativeOrder(order) - # apfel.SetAlphaEvolution("exact") - # apfel.SetAlphaQCDRef(theory_dict["alphas"], theory_dict["Qref"]) - # apfel.SetVFNS() - # apfel.SetMSbarMasses( - # theory_dict["mc"], theory_dict["mb"], theory_dict["mt"] - # ) - # apfel.SetMassScaleReference( - # theory_dict["Qmc"], theory_dict["Qmb"], theory_dict["Qmt"] - # ) - # apfel.SetMassMatchingScales( - # theory_dict["kcThr"], theory_dict["kbThr"], theory_dict["ktThr"] - # ) - # apfel.SetRenFacRatio(theory_dict["fact_to_ren_scale_ratio"]) - # apfel.InitializeAPFEL() - # apfel.EnableWelcomeMessage(1) - # #check myself to APFEL - # np.testing.assert_allclose( - # apfel_vals, np.sqrt(np.array(my_masses)), rtol=5e-5 - # ) + # Apfel bottom masses are not the same + np.testing.assert_raises( + AssertionError, + np.testing.assert_allclose, + apfel_masses_thr[0], + apfel_masses_plain[0], + ) + np.testing.assert_raises( + AssertionError, + np.testing.assert_allclose, + apfel_masses_thr[0], + apfel_masses_plain[0], + ) diff --git a/tests/test_msbar_masses.py b/tests/test_msbar_masses.py index a0093a884..32a1a9f4d 100644 --- a/tests/test_msbar_masses.py +++ b/tests/test_msbar_masses.py @@ -56,7 +56,7 @@ def test_compute_msbar_mass(self): q2_to=m2_computed[nf - 3], ) ) - np.testing.assert_allclose(m2_computed, m2_test, rtol=9e-3 if order==3 else 5e-3) + np.testing.assert_allclose(m2_computed, m2_test, rtol=4e-3) def test_compute_msbar_mass_VFNS(self): # test the solution now with some initial contition @@ -64,7 +64,7 @@ def test_compute_msbar_mass_VFNS(self): theory_dict.update( { "ModEv": "TRN", - "PTO": 2, + "PTO": 3, "mc": 2.0, "mb": 4.0, "Qmc": 80.0, @@ -88,7 +88,7 @@ def test_compute_msbar_mass_VFNS(self): q2_to=m2_computed[nf - 3], ) ) - np.testing.assert_allclose(m2_computed, m2_test, rtol=5e-3) + np.testing.assert_allclose(m2_computed, m2_test, rtol=3e-3) def test_errors(self): @@ -103,7 +103,7 @@ def test_errors(self): # ) # ) # compute_msbar_mass(theory_dict) - with pytest.raises(ValueError, match="masses need to be sorted"): + with pytest.raises(ValueError, match="Msbar masses are not to be sorted"): theory_dict.update( dict( mc=1.1, From 829afe1dbcb39252c3f690ccf43459c02ed270eb Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Thu, 10 Feb 2022 11:56:20 +0100 Subject: [PATCH 17/24] Removing nf annotations --- src/eko/anomalous_dimensions/harmonics.py | 2 +- src/eko/anomalous_dimensions/lo.py | 6 +++--- src/eko/anomalous_dimensions/nlo.py | 14 +++++++------- src/eko/anomalous_dimensions/nnlo.py | 16 ++++++++-------- src/eko/beta.py | 8 ++++---- src/eko/gamma.py | 8 ++++---- 6 files changed, 27 insertions(+), 27 deletions(-) diff --git a/src/eko/anomalous_dimensions/harmonics.py b/src/eko/anomalous_dimensions/harmonics.py index 6270f0f0b..75163b781 100644 --- a/src/eko/anomalous_dimensions/harmonics.py +++ b/src/eko/anomalous_dimensions/harmonics.py @@ -17,7 +17,7 @@ @nb.njit("c16(c16,u1)", cache=True) -def cern_polygamma(Z, K: int): # pylint: disable=all +def cern_polygamma(Z, K): # pylint: disable=all """ Computes the polygamma functions :math:`\\psi_k(z)`. diff --git a/src/eko/anomalous_dimensions/lo.py b/src/eko/anomalous_dimensions/lo.py index fc2bf9e71..92dfe081c 100644 --- a/src/eko/anomalous_dimensions/lo.py +++ b/src/eko/anomalous_dimensions/lo.py @@ -32,7 +32,7 @@ def gamma_ns_0(N, s1): @nb.njit("c16(c16,u1)", cache=True) -def gamma_qg_0(N, nf: int): +def gamma_qg_0(N, nf): """ Computes the leading-order quark-gluon anomalous dimension @@ -78,7 +78,7 @@ def gamma_gq_0(N): @nb.njit("c16(c16,c16,u1)", cache=True) -def gamma_gg_0(N, s1, nf: int): +def gamma_gg_0(N, s1, nf): """ Computes the leading-order gluon-gluon anomalous dimension @@ -104,7 +104,7 @@ def gamma_gg_0(N, s1, nf: int): @nb.njit("c16[:,:](c16,c16,u1)", cache=True) -def gamma_singlet_0(N, s1, nf: int): +def gamma_singlet_0(N, s1, nf): r""" Computes the leading-order singlet anomalous dimension matrix diff --git a/src/eko/anomalous_dimensions/nlo.py b/src/eko/anomalous_dimensions/nlo.py index e4ec996c9..3466be105 100644 --- a/src/eko/anomalous_dimensions/nlo.py +++ b/src/eko/anomalous_dimensions/nlo.py @@ -15,7 +15,7 @@ @nb.njit("c16(c16,u1)", cache=True) -def gamma_nsm_1(n, nf: int): +def gamma_nsm_1(n, nf): """ Computes the |NLO| valence-like non-singlet anomalous dimension. @@ -58,7 +58,7 @@ def gamma_nsm_1(n, nf: int): @nb.njit("c16(c16,u1)", cache=True) -def gamma_nsp_1(n, nf: int): +def gamma_nsp_1(n, nf): """ Computes the |NLO| singlet-like non-singlet anomalous dimension. @@ -99,7 +99,7 @@ def gamma_nsp_1(n, nf: int): @nb.njit("c16(c16,u1)", cache=True) -def gamma_ps_1(n, nf: int): +def gamma_ps_1(n, nf): """ Computes the |NLO| pure-singlet quark-quark anomalous dimension. @@ -126,7 +126,7 @@ def gamma_ps_1(n, nf: int): @nb.njit("c16(c16,u1)", cache=True) -def gamma_qg_1(n, nf: int): +def gamma_qg_1(n, nf): """ Computes the |NLO| quark-gluon singlet anomalous dimension. @@ -159,7 +159,7 @@ def gamma_qg_1(n, nf: int): @nb.njit("c16(c16,u1)", cache=True) -def gamma_gq_1(n, nf: int): +def gamma_gq_1(n, nf): """ Computes the |NLO| gluon-quark singlet anomalous dimension. @@ -195,7 +195,7 @@ def gamma_gq_1(n, nf: int): @nb.njit("c16(c16,u1)", cache=True) -def gamma_gg_1(n, nf: int): +def gamma_gg_1(n, nf): """ Computes the |NLO| gluon-gluon singlet anomalous dimension. @@ -233,7 +233,7 @@ def gamma_gg_1(n, nf: int): @nb.njit("c16[:,:](c16,u1)", cache=True) -def gamma_singlet_1(N, nf: int): +def gamma_singlet_1(N, nf): r""" Computes the next-leading-order singlet anomalous dimension matrix diff --git a/src/eko/anomalous_dimensions/nnlo.py b/src/eko/anomalous_dimensions/nnlo.py index 5cf5f1f52..c99da71fe 100644 --- a/src/eko/anomalous_dimensions/nnlo.py +++ b/src/eko/anomalous_dimensions/nnlo.py @@ -17,7 +17,7 @@ @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_nsm_2(n, nf: int, sx): +def gamma_nsm_2(n, nf, sx): """ Computes the |NNLO| valence-like non-singlet anomalous dimension. @@ -94,7 +94,7 @@ def gamma_nsm_2(n, nf: int, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_nsp_2(n, nf: int, sx): +def gamma_nsp_2(n, nf, sx): """ Computes the |NNLO| singlet-like non-singlet anomalous dimension. @@ -171,7 +171,7 @@ def gamma_nsp_2(n, nf: int, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_nsv_2(n, nf: int, sx): +def gamma_nsv_2(n, nf, sx): """ Computes the |NNLO| valence non-singlet anomalous dimension. @@ -226,7 +226,7 @@ def gamma_nsv_2(n, nf: int, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_ps_2(n, nf: int, sx): +def gamma_ps_2(n, nf, sx): """ Computes the |NNLO| pure-singlet quark-quark anomalous dimension. @@ -298,7 +298,7 @@ def gamma_ps_2(n, nf: int, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_qg_2(n, nf: int, sx): +def gamma_qg_2(n, nf, sx): """ Computes the |NNLO| quark-gluon singlet anomalous dimension. @@ -372,7 +372,7 @@ def gamma_qg_2(n, nf: int, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_gq_2(n, nf: int, sx): +def gamma_gq_2(n, nf, sx): """ Computes the |NNLO| gluon-quark singlet anomalous dimension. @@ -462,7 +462,7 @@ def gamma_gq_2(n, nf: int, sx): @nb.njit("c16(c16,u1,c16[:])", cache=True) -def gamma_gg_2(n, nf: int, sx): +def gamma_gg_2(n, nf, sx): """ Computes the |NNLO| gluon-gluon singlet anomalous dimension. @@ -550,7 +550,7 @@ def gamma_gg_2(n, nf: int, sx): @nb.njit("c16[:,:](c16,u1,c16[:])", cache=True) -def gamma_singlet_2(N, nf: int, sx): +def gamma_singlet_2(N, nf, sx): r""" Computes the |NNLO| singlet anomalous dimension matrix diff --git a/src/eko/beta.py b/src/eko/beta.py index b7ed021ad..445c44f47 100644 --- a/src/eko/beta.py +++ b/src/eko/beta.py @@ -12,7 +12,7 @@ @nb.njit("f8(u1)", cache=True) -def beta_0(nf: int): +def beta_0(nf): """ Computes the first coefficient of the QCD beta function. @@ -33,7 +33,7 @@ def beta_0(nf: int): @nb.njit("f8(u1)", cache=True) -def beta_1(nf: int): +def beta_1(nf): """ Computes the second coefficient of the QCD beta function. @@ -58,7 +58,7 @@ def beta_1(nf: int): @nb.njit("f8(u1)", cache=True) -def beta_2(nf: int): +def beta_2(nf): """ Computes the third coefficient of the QCD beta function @@ -87,7 +87,7 @@ def beta_2(nf: int): @nb.njit("f8(u1)", cache=True) -def beta_3(nf: int): +def beta_3(nf): """ Computes the fourth coefficient of the QCD beta function diff --git a/src/eko/gamma.py b/src/eko/gamma.py index aa7e1346d..2abed82b5 100644 --- a/src/eko/gamma.py +++ b/src/eko/gamma.py @@ -25,7 +25,7 @@ def gamma_0(): @nb.njit("f8(u1)", cache=True) -def gamma_1(nf: int): +def gamma_1(nf): """ Computes the second coefficient of the QCD gamma function. @@ -45,7 +45,7 @@ def gamma_1(nf: int): @nb.njit("f8(u1)", cache=True) -def gamma_2(nf: int): +def gamma_2(nf): """ Computes the third coefficient of the QCD gamma function. @@ -65,7 +65,7 @@ def gamma_2(nf: int): @nb.njit("f8(u1)", cache=True) -def gamma_3(nf: int): +def gamma_3(nf): """ Computes the fourth coefficient of the QCD gamma function. @@ -98,7 +98,7 @@ def gamma_3(nf: int): @nb.njit("f8(u1,u1)", cache=True) -def gamma(order, nf: int): +def gamma(order, nf): """ Compute the value of a gamma coefficient From 8709be92e5a9f380572d606a8787df15f8fe986b Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Thu, 10 Feb 2022 11:57:56 +0100 Subject: [PATCH 18/24] remove _c form coeff mathing --- src/eko/msbar_masses.py | 4 ++-- src/eko/strong_coupling.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py index c4d12d3a4..a2314b64b 100644 --- a/src/eko/msbar_masses.py +++ b/src/eko/msbar_masses.py @@ -218,8 +218,8 @@ def compute_matching_coeffs_down(nf): matching_coeffs_down: downward matching coefficient matrix """ - _c = compute_matching_coeffs_up(nf) - return invert_matching_coeffs(_c) + c_up = compute_matching_coeffs_up(nf) + return invert_matching_coeffs(c_up) def solve_msbar_mass( diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index a4f6a7a59..c7983f519 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -477,8 +477,8 @@ def compute_matching_coeffs_down(mass_scheme, nf): matching_coeffs_down: downward matching coefficient matrix """ - _c = compute_matching_coeffs_up(mass_scheme, nf) - return invert_matching_coeffs(_c) + c_up = compute_matching_coeffs_up(mass_scheme, nf) + return invert_matching_coeffs(c_up) def invert_matching_coeffs(c_up): From eaad75ca4f7ecde6302340d28b8787ec26671ee4 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Thu, 10 Feb 2022 12:45:47 +0100 Subject: [PATCH 19/24] Make msbar function names more explicit --- src/eko/msbar_masses.py | 34 ++++++++++++++---------------- src/eko/runner.py | 4 ++-- tests/benchmark_msbar_evolution.py | 10 ++++----- tests/test_msbar_masses.py | 20 +++++++++--------- 4 files changed, 33 insertions(+), 35 deletions(-) diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py index a2314b64b..e57b7a726 100644 --- a/src/eko/msbar_masses.py +++ b/src/eko/msbar_masses.py @@ -12,7 +12,7 @@ from .strong_coupling import StrongCoupling, invert_matching_coeffs -def msbar_ker_exact(a0, a1, order, nf): +def ker_exact(a0, a1, order, nf): r""" Exact |MSbar| |RGE| kernel @@ -29,12 +29,12 @@ def msbar_ker_exact(a0, a1, order, nf): Returns ------- - ker: float + ker_exact: float Exact |MSbar| kernel: .. math:: - k_{exact} = e^{- \int_{a_s(\mu_{h,0}^2)}^{a_s(\mu^2)} \gamma_m(a_s) / \beta(a_s) da_s} - """ # pylint:disable=line-too-long + k_{exact} = e^{-\int_{a_s(\mu_{h,0}^2)}^{a_s(\mu^2)}\gamma_m(a_s)/ \beta(a_s)da_s} + """ b_vec = [beta(0, nf)] g_vec = [gamma(0, nf)] if order >= 1: @@ -69,7 +69,7 @@ def integrand(a, b_vec, g_vec): @nb.njit("f8(f8,f8,u1,u1)", cache=True) -def msbar_ker_expanded(a0, a1, order, nf): +def ker_expanded(a0, a1, order, nf): r""" Expanded |MSbar| |RGE| kernel @@ -86,7 +86,7 @@ def msbar_ker_expanded(a0, a1, order, nf): Returns ------- - ker: float + ker_expanded: float Expanded |MSbar| kernel: .. math:: @@ -140,7 +140,7 @@ def msbar_ker_expanded(a0, a1, order, nf): return ev_mass * num / den -def msbar_ker_dispatcher(q2_to, q2m_ref, strong_coupling, fact_to_ren, nf): +def ker_dispatcher(q2_to, q2m_ref, strong_coupling, fact_to_ren, nf): r""" Select the |MSbar| kernel and compute the strong coupling values @@ -168,8 +168,8 @@ def msbar_ker_dispatcher(q2_to, q2m_ref, strong_coupling, fact_to_ren, nf): method = strong_coupling.method order = strong_coupling.order if method == "expanded": - return msbar_ker_expanded(a0, a1, order, nf) - return msbar_ker_exact(a0, a1, order, nf) + return ker_expanded(a0, a1, order, nf) + return ker_exact(a0, a1, order, nf) def compute_matching_coeffs_up(nf): @@ -222,7 +222,7 @@ def compute_matching_coeffs_down(nf): return invert_matching_coeffs(c_up) -def solve_msbar_mass( +def solve( m2_ref, q2m_ref, strong_coupling, @@ -256,8 +256,7 @@ def solve_msbar_mass( def rge(m2, q2m_ref, strong_coupling, fact_to_ren, nf_ref): return ( m2_ref - * msbar_ker_dispatcher(m2, q2m_ref, strong_coupling, fact_to_ren, nf_ref) - ** 2 + * ker_dispatcher(m2, q2m_ref, strong_coupling, fact_to_ren, nf_ref) ** 2 - m2 ) @@ -267,7 +266,7 @@ def rge(m2, q2m_ref, strong_coupling, fact_to_ren, nf_ref): return float(msbar_mass) -def evolve_msbar_mass( +def evolve( m2_ref, q2m_ref, strong_coupling, @@ -336,13 +335,12 @@ def evolve_msbar_mass( fact += as_thr ** pto * L ** l * m_coeffs[pto, l] ev_mass *= ( fact - * msbar_ker_dispatcher(q2_final, q2_init, strong_coupling, fact_to_ren, nf) - ** 2 + * ker_dispatcher(q2_final, q2_init, strong_coupling, fact_to_ren, nf) ** 2 ) return m2_ref * ev_mass -def compute_msbar_mass(theory_card): +def compute(theory_card): r""" Compute the |MSbar| masses solving the equation :math:`m_{\bar{MS}}(\mu) = \mu` for each heavy quark and consistent boundary contitions. @@ -425,7 +423,7 @@ def sc(thr_masses): # len(masses[q2m_ref > masses]) + 3 is the nf at the given reference scale if nf_target != len(masses[q2m_ref > masses]) + 3: q2_to = masses[q_idx + shift] - m2_ref = evolve_msbar_mass( + m2_ref = evolve( m2_ref, q2m_ref, sc(masses), @@ -435,7 +433,7 @@ def sc(thr_masses): q2m_ref = q2_to # now solve the RGE - masses[q_idx] = solve_msbar_mass( + masses[q_idx] = solve( m2_ref, q2m_ref, sc(masses), diff --git a/src/eko/runner.py b/src/eko/runner.py index 2ee4ed554..863c9ab22 100644 --- a/src/eko/runner.py +++ b/src/eko/runner.py @@ -9,8 +9,8 @@ from . import basis_rotation as br from . import interpolation +from . import msbar_masses from .evolution_operator.grid import OperatorGrid -from .msbar_masses import compute_msbar_mass from .output import Output from .strong_coupling import StrongCoupling from .thresholds import ThresholdsAtlas @@ -49,7 +49,7 @@ def __init__(self, theory_card, operators_card): # setup the Threshold path, compute masses if necessary masses = None if theory_card["HQ"] == "MSBAR": - masses = compute_msbar_mass(theory_card) + masses = msbar_masses.compute(theory_card) tc = ThresholdsAtlas.from_dict(theory_card, masses=masses) self.out["q2_ref"] = float(tc.q2_ref) diff --git a/tests/benchmark_msbar_evolution.py b/tests/benchmark_msbar_evolution.py index d8021381b..bff7bebe9 100644 --- a/tests/benchmark_msbar_evolution.py +++ b/tests/benchmark_msbar_evolution.py @@ -2,7 +2,7 @@ """This module benchmarks MSbar mass evolution against APFEL.""" import numpy as np -from eko.msbar_masses import compute_msbar_mass, evolve_msbar_mass, compute_msbar_mass +from eko import msbar_masses from eko.strong_coupling import StrongCoupling # try to load APFEL - if not available, we'll use the cached values @@ -109,7 +109,7 @@ def benchmark_APFEL_msbar_evolution(self): my_masses = [] for n in [3, 4, 5]: my_masses.append( - evolve_msbar_mass( + msbar_masses.evolve( m2[n - 3], Q2m[n - 3], strong_coupling=as_VFNS, @@ -163,7 +163,7 @@ def benchmark_APFEL_msbar_solution(self): # collect my values for order in [0, 1, 2]: theory_dict["PTO"] = order - my_masses = compute_msbar_mass(theory_dict) + my_masses = msbar_masses.compute(theory_dict) # get APFEL numbers - if available else use cache apfel_vals = apfel_vals_dict[theory_dict["ModEv"]][order] if use_APFEL: @@ -209,7 +209,7 @@ def benchmark_msbar_solution_kthr(self): "PTO": 0, } ) - my_masses_thr = compute_msbar_mass(theory_dict) + my_masses_thr = msbar_masses.compute(theory_dict) apfel_masses_thr = [1.9891, 4.5102, 175.0000] theory_dict.update( { @@ -217,7 +217,7 @@ def benchmark_msbar_solution_kthr(self): "kbThr": 1.0, } ) - my_masses_plain = compute_msbar_mass(theory_dict) + my_masses_plain = msbar_masses.compute(theory_dict) apfel_masses_plain = ([1.9855, 4.8062, 175.0000],) # Eko bottom mass is the same diff --git a/tests/test_msbar_masses.py b/tests/test_msbar_masses.py index 32a1a9f4d..097177a30 100644 --- a/tests/test_msbar_masses.py +++ b/tests/test_msbar_masses.py @@ -5,8 +5,8 @@ import numpy as np import pytest +from eko import msbar_masses from eko.evolution_operator.flavors import quark_names -from eko.msbar_masses import compute_msbar_mass, evolve_msbar_mass from eko.strong_coupling import StrongCoupling theory_dict = { @@ -41,14 +41,14 @@ def test_compute_msbar_mass(self): strong_coupling = StrongCoupling.from_dict(theory_dict) # compute the scale such msbar(m) = m - m2_computed = compute_msbar_mass(theory_dict) + m2_computed = msbar_masses.compute(theory_dict) m2_test = [] for nf in [3, 4, 5]: # compute msbar( m ) m2_ref = theory_dict[f"m{quark_names[nf]}"] ** 2 Q2m_ref = theory_dict[f"Qm{quark_names[nf]}"] ** 2 m2_test.append( - evolve_msbar_mass( + msbar_masses.evolve( m2_ref, Q2m_ref, strong_coupling=strong_coupling, @@ -73,14 +73,14 @@ def test_compute_msbar_mass_VFNS(self): ) strong_coupling = StrongCoupling.from_dict(theory_dict) # compute the scale such msbar(m) = m - m2_computed = compute_msbar_mass(theory_dict) + m2_computed = msbar_masses.compute(theory_dict) m2_test = [] for nf in [3, 4, 5]: # compute msbar( m ) m2_ref = theory_dict[f"m{quark_names[nf]}"] ** 2 Q2m_ref = theory_dict[f"Qm{quark_names[nf]}"] ** 2 m2_test.append( - evolve_msbar_mass( + msbar_masses.evolve( m2_ref, Q2m_ref, strong_coupling=strong_coupling, @@ -112,26 +112,26 @@ def test_errors(self): Qmb=1.0, ) ) - compute_msbar_mass(theory_dict) + msbar_masses.compute(theory_dict) # test forward conditions on alphas_ref with pytest.raises(ValueError, match="should be lower than"): theory_dict.update(dict(Qmb=91.0001)) - compute_msbar_mass(theory_dict) + msbar_masses.compute(theory_dict) # test backward conditions on alphas_ref with pytest.raises(ValueError, match="should be greater than"): theory_dict.update(dict(Qmt=89.9999)) - compute_msbar_mass(theory_dict) + msbar_masses.compute(theory_dict) theory_dict.update(dict(Qmb=4.0, Qmt=175)) # test forward conditions on masses with pytest.raises(ValueError, match="should be lower than m"): theory_dict.update(dict(mt=174)) - compute_msbar_mass(theory_dict) + msbar_masses.compute(theory_dict) # test backward conditions on masses with pytest.raises(ValueError, match="should be greater than m"): theory_dict.update(dict(mb=4.1, mt=176)) - compute_msbar_mass(theory_dict) + msbar_masses.compute(theory_dict) From e4a80bcdfb6461632e89c47272b0370c1e35d3fe Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Thu, 10 Feb 2022 12:54:24 +0100 Subject: [PATCH 20/24] More pylint fixes and docs spelling --- benchmarks/apfel_bench.py | 4 ++-- doc/source/theory/pQCD.rst | 2 +- src/eko/msbar_masses.py | 4 ++-- src/eko/output.py | 10 +++++----- src/eko/strong_coupling.py | 15 ++++++++------- 5 files changed, 18 insertions(+), 17 deletions(-) diff --git a/benchmarks/apfel_bench.py b/benchmarks/apfel_bench.py index f7e57ca3f..5e73f7e31 100644 --- a/benchmarks/apfel_bench.py +++ b/benchmarks/apfel_bench.py @@ -93,8 +93,8 @@ def benchmark_msbar(self, pto): MSbar heavy quark mass scheme when passing kthr != 1 both apfel and eko use ``kThr * msbar``, as thr scale, where ``msbar`` is the usual ms_bar solution. - However apfel and eko mange the alpha_s thr differently - (apfel uses the given mass parameters as thr multipied by the kthr), + However apfel and eko manage the alpha_s thr differently: + apfel uses the given mass parameters as thr multiplied by the kthr, while in eko only the already computed thr matters, so the benchmark is not a proper comparison with this option allowed. """ diff --git a/doc/source/theory/pQCD.rst b/doc/source/theory/pQCD.rst index 967f7018e..7f892f947 100644 --- a/doc/source/theory/pQCD.rst +++ b/doc/source/theory/pQCD.rst @@ -43,7 +43,7 @@ being: When the renormalization scale crosses a flavor threshold matching conditions have to be applied :cite:`Schroder:2005hy,Chetyrkin:2005ia`. -In particular the matching when passing form :math:`n_f` to :math:`n_f-1` schemes +In particular, the matching involved in the change from :math:`n_f` to :math:`n_f-1` schemes is presented in equation 3.1 of :cite:`Schroder:2005hy` for |MSbar| masses, while the same expression for POLE masses is reported in Appendix A. diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py index e57b7a726..2f07749e0 100644 --- a/src/eko/msbar_masses.py +++ b/src/eko/msbar_masses.py @@ -245,7 +245,7 @@ def solve( nf_ref: int number of active flavours at the scale q2m_ref, where the solution is searched fact_to_ren: float - :math:`\mu_F^2/\muR^2` + :math:`\mu_F^2/\mu_R^2` Returns ------- @@ -287,7 +287,7 @@ def evolve( Instance of :class:`~eko.strong_coupling.StrongCoupling` able to generate a_s for any q fact_to_ren: float - :math:`\mu_F^2/\muR^2` + :math:`\mu_F^2/\mu_R^2` q2_to: float, optional scale at which the mass is computed diff --git a/src/eko/output.py b/src/eko/output.py index f6b7aa963..28c004926 100644 --- a/src/eko/output.py +++ b/src/eko/output.py @@ -11,7 +11,7 @@ import numpy as np import yaml -from lz4 import frame +from lz4 import frame as lz4_frame from . import basis_rotation as br from . import interpolation, version @@ -308,7 +308,7 @@ def get_raw(self, binarize=True, skip_q2_grid=False): out["Q2grid"][q2][k] = float(v) continue if binarize: - out["Q2grid"][q2][k] = frame.compress(v.tobytes()) + out["Q2grid"][q2][k] = lz4_frame.compress(v.tobytes()) else: out["Q2grid"][q2][k] = v.tolist() else: @@ -393,7 +393,7 @@ def dump_tar(self, tarname): stream = io.BytesIO() np.save(stream, operator) stream.seek(0) - with frame.open( + with lz4_frame.open( (tmpdir / kind).with_suffix(".npy.lz4"), "wb" ) as fo: fo.write(stream.read()) @@ -436,7 +436,7 @@ def load_yaml(cls, stream, skip_q2_grid=False): elif isinstance(v, list): v = np.array(v) elif isinstance(v, bytes): - v = np.frombuffer(frame.decompress(v)) + v = np.frombuffer(lz4_frame.decompress(v)) v = v.reshape(len_tpids, len_tgrid, len_ipids, len_igrid) op[k] = v return cls(obj) @@ -500,7 +500,7 @@ def load_tar(cls, tarname): grids = {} for fp in innerdir.glob("*.npy.lz4"): - with frame.open(fp, "rb") as fd: + with lz4_frame.open(fp, "rb") as fd: stream = io.BytesIO(fd.read()) stream.seek(0) grids[pathlib.Path(fp.stem).stem] = np.load(stream) diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index c7983f519..dc77eabbc 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -13,7 +13,8 @@ import scipy from . import constants, thresholds -from .beta import b, beta +from .beta import b as beta_b +from .beta import beta logger = logging.getLogger(__name__) @@ -67,12 +68,12 @@ def as_expanded(order, as_ref, nf, scale_from, scale_to): res = as_LO # NLO if order >= 1: - b1 = b(1, nf) + b1 = beta_b(1, nf) as_NLO = as_LO * (1 - b1 * as_LO * np.log(den)) res = as_NLO # NNLO if order >= 2: - b2 = b(2, nf) + b2 = beta_b(2, nf) res = as_LO * ( 1.0 + as_LO * (as_LO - as_ref) * (b2 - b1 ** 2) @@ -81,7 +82,7 @@ def as_expanded(order, as_ref, nf, scale_from, scale_to): # N3LO expansion is taken from Luca Rottoli and it's simpler # than the APFEL expansion which is more accurate if order >= 3: - b3 = b(3, nf) + b3 = beta_b(3, nf) log_fact = np.log(as_LO) res += ( as_LO ** 4 @@ -295,13 +296,13 @@ def compute_exact(self, as_ref, nf, scale_from, scale_to): b_vec = [1] # NLO if self.order >= 1: - b_vec.append(b(1, nf)) + b_vec.append(beta_b(1, nf)) # NNLO if self.order >= 2: - b_vec.append(b(2, nf)) + b_vec.append(beta_b(2, nf)) # N3LO if self.order >= 3: - b_vec.append(b(3, nf)) + b_vec.append(beta_b(3, nf)) # integration kernel def rge(_t, a, b_vec): return -(a ** 2) * np.sum([a ** k * b for k, b in enumerate(b_vec)]) From 5dd80cfcb11c5d88a89a5e75eb21a703d3ac28c0 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Thu, 10 Feb 2022 13:06:04 +0100 Subject: [PATCH 21/24] Erasing msbar complex test and docs --- doc/source/theory/pQCD.rst | 7 ++++--- src/eko/msbar_masses.py | 18 ------------------ tests/benchmark_strong_coupling.py | 4 ++-- tests/test_interpolation.py | 8 ++++---- tests/test_msbar_masses.py | 10 ---------- 5 files changed, 10 insertions(+), 37 deletions(-) diff --git a/doc/source/theory/pQCD.rst b/doc/source/theory/pQCD.rst index 7f892f947..e73cf3c38 100644 --- a/doc/source/theory/pQCD.rst +++ b/doc/source/theory/pQCD.rst @@ -174,10 +174,11 @@ To find coherent solutions and perform the mass running in the correct patches i is necessary to always start computing the mass scales closer to :math:`\mu_{ref}`. Eventually, to ensure that the threshold values are properly set, we add a -consistency check, asserting: +consistency check, asserting that the :math:`m_{\overline{MS},h}` are properly sorted. -.. math :: - m_{\overline{MS},h} (m_h) \leq m_{\overline{MS},h+1} (m_h) +Note that also for |MSbar| mass running when the heavy threshold scales are +crossed we need to apply non trivial matching from order +:math:`\mathcal{O}(a_s^2)` as described here :cite:`Liu:2015fxa`. We provide the following as an illustrative example of how this procedure works: when the strong coupling is given with boundary condition :math:`\alpha_s(\mu_{ref}=91, n_{f_{ref}}=5)` diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py index 2f07749e0..cb9231481 100644 --- a/src/eko/msbar_masses.py +++ b/src/eko/msbar_masses.py @@ -441,25 +441,7 @@ def sc(thr_masses): fact_to_ren, ) - # TODO: this test seems to be quite hard to contradict - # what if we check just the mass ordering? # Check the msbar ordering if not (masses == np.sort(masses)).all(): raise ValueError("Msbar masses are not to be sorted") - # for m2_msbar, hq in zip(masses[:-1], quark_names[4:]): - # q2m_ref = np.power(theory_card[f"Qm{hq}"], 2) - # m2_ref = np.power(theory_card[f"m{hq}"], 2) - # # check that m_msbar_hq < msbar_hq+1 (m_msbar_hq) - # m2_test = evolve_msbar_mass( - # m2_ref, - # q2m_ref, - # sc(masses), - # fact_to_ren, - # m2_msbar, - # ) - # if m2_msbar >= m2_test: - # raise ValueError( - # "The MSBAR masses do not preserve the correct ordering,\ - # check the initial reference values" - # ) return masses diff --git a/tests/benchmark_strong_coupling.py b/tests/benchmark_strong_coupling.py index e42a73b64..8c1848102 100644 --- a/tests/benchmark_strong_coupling.py +++ b/tests/benchmark_strong_coupling.py @@ -170,7 +170,7 @@ def benchmark_pegasus_ffns(self): 0.02224042559454458, 0.014988806602725375, 0.01140950716164355, - 0.009245168439298405 + 0.009245168439298405, ] ), } @@ -688,7 +688,7 @@ def benchmark_lhapdf_exact(self): 0.025955428065113105, 0.015862752802768762, 0.011614793662449371, - 0.009214009540864635 + 0.009214009540864635, ] ), } diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index c4e06b7b6..a58c5e60d 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -240,11 +240,11 @@ def test_log_eval_N(self): assert_almost_equal(act_c, res_c) # Full -> \tilde p_0(N) = exp(-N)(exp(N)-1-N)/N^2 # MMa: Integrate[x^(n-1) (-Log[x]),{x,1/E,1}] - p0Nref_full = lambda N, lnx: ((np.exp(N) - 1 - N) / N**2) * np.exp( + p0Nref_full = lambda N, lnx: ((np.exp(N) - 1 - N) / N ** 2) * np.exp( -N * (lnx + 1) ) # partial = lower bound is neglected; - p0Nref_partial = lambda N, lnx: (1 / N**2) * np.exp(-N * lnx) + p0Nref_partial = lambda N, lnx: (1 / N ** 2) * np.exp(-N * lnx) p1N = inter_N[1] assert len(p1N.areas) == 1 p1_cs_ref = [1, 1] @@ -252,8 +252,8 @@ def test_log_eval_N(self): assert_almost_equal(act_c, res_c) # p_1(x) = 1+\ln(x) -> \tilde p_1(N) = (exp(-N)-1+N)/N^2 # MMa: Integrate[x^(n-1) (1+Log[x]),{x,1/E,1}] - p1Nref_full = lambda N, lnx: ((np.exp(-N) - 1 + N) / N**2) * np.exp(-N * lnx) - p1Nref_partial = lambda N, lnx: (1 / N - 1 / N**2) * np.exp(-N * lnx) + p1Nref_full = lambda N, lnx: ((np.exp(-N) - 1 + N) / N ** 2) * np.exp(-N * lnx) + p1Nref_partial = lambda N, lnx: (1 / N - 1 / N ** 2) * np.exp(-N * lnx) # iterate configurations for N in [1.0, 2.0, complex(1.0, 1.0)]: # check skip diff --git a/tests/test_msbar_masses.py b/tests/test_msbar_masses.py index 097177a30..fbe530534 100644 --- a/tests/test_msbar_masses.py +++ b/tests/test_msbar_masses.py @@ -93,16 +93,6 @@ def test_compute_msbar_mass_VFNS(self): def test_errors(self): # test mass ordering - # with pytest.raises(ValueError, match="do not preserve the correct ordering"): - # theory_dict.update( - # dict( - # mc=1.009, - # mb=1.01, - # Qmc=1.012, - # Qmb=1.01, - # ) - # ) - # compute_msbar_mass(theory_dict) with pytest.raises(ValueError, match="Msbar masses are not to be sorted"): theory_dict.update( dict( From 32698b4043f7a24377f32f97c44828b07d66e650 Mon Sep 17 00:00:00 2001 From: Alessandro Candido Date: Tue, 15 Feb 2022 19:09:34 +0100 Subject: [PATCH 22/24] Run pre-commit on all files --- doc/source/refs.bib | 2 +- doc/source/theory/pQCD.rst | 4 ++-- src/eko/beta.py | 7 ++++--- src/eko/runner.py | 3 +-- src/eko/strong_coupling.py | 24 ++++++++++++------------ tests/test_beta.py | 2 +- tests/test_interpolation.py | 8 ++++---- tests/test_strong_coupling.py | 4 ++-- 8 files changed, 27 insertions(+), 27 deletions(-) diff --git a/doc/source/refs.bib b/doc/source/refs.bib index a46b8394f..30bebb630 100644 --- a/doc/source/refs.bib +++ b/doc/source/refs.bib @@ -445,4 +445,4 @@ @article{Liu:2015fxa volume = "746", pages = "330--334", year = "2015" -} \ No newline at end of file +} diff --git a/doc/source/theory/pQCD.rst b/doc/source/theory/pQCD.rst index e73cf3c38..b0c4078fb 100644 --- a/doc/source/theory/pQCD.rst +++ b/doc/source/theory/pQCD.rst @@ -34,7 +34,7 @@ We implement two different strategies to solve the |RGE|: & - 2 b_1 L_{0} \left[ 5 b_1^2 L_{\text{LO}} + 3 b_1^2 L_{\text{LO}}^2 + b_0 \left[b_2 + 2 \left(b_1^2 - b_0 b_2\right) L_{\mu} a_s(\mu_0^2)\right] \right ] \\ & \left. \right\} -being: +being: .. math :: L_{\mu} &= \ln(\mu_R^2/\mu_0^2) \\ @@ -43,7 +43,7 @@ being: When the renormalization scale crosses a flavor threshold matching conditions have to be applied :cite:`Schroder:2005hy,Chetyrkin:2005ia`. -In particular, the matching involved in the change from :math:`n_f` to :math:`n_f-1` schemes +In particular, the matching involved in the change from :math:`n_f` to :math:`n_f-1` schemes is presented in equation 3.1 of :cite:`Schroder:2005hy` for |MSbar| masses, while the same expression for POLE masses is reported in Appendix A. diff --git a/src/eko/beta.py b/src/eko/beta.py index 445c44f47..dcdeda6de 100644 --- a/src/eko/beta.py +++ b/src/eko/beta.py @@ -105,9 +105,10 @@ def beta_3(nf): """ beta_3 = ( 149753.0 / 6.0 - + 3564.0 * zeta3 + nf * (-1078361.0 / 162.0 - 6508.0 / 27.0 * zeta3) - + nf ** 2 * (50065.0 / 162.0 + 6472.0 / 81.0 * zeta3) - + 1093.0 / 729.0 * nf ** 3 + + 3564.0 * zeta3 + + nf * (-1078361.0 / 162.0 - 6508.0 / 27.0 * zeta3) + + nf**2 * (50065.0 / 162.0 + 6472.0 / 81.0 * zeta3) + + 1093.0 / 729.0 * nf**3 ) return beta_3 diff --git a/src/eko/runner.py b/src/eko/runner.py index 863c9ab22..1e4f59601 100644 --- a/src/eko/runner.py +++ b/src/eko/runner.py @@ -8,8 +8,7 @@ import numpy as np from . import basis_rotation as br -from . import interpolation -from . import msbar_masses +from . import interpolation, msbar_masses from .evolution_operator.grid import OperatorGrid from .output import Output from .strong_coupling import StrongCoupling diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index 18e5fe815..6d8c9efae 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -85,25 +85,25 @@ def as_expanded(order, as_ref, nf, scale_from, scale_to): b3 = beta_b(3, nf) log_fact = np.log(as_LO) res += ( - as_LO ** 4 - / (2 * beta0 ** 3) + as_LO**4 + / (2 * beta0**3) * ( - -2 * b1 ** 3 * np.log(as_ref) ** 3 - + 5 * b1 ** 3 * log_fact ** 2 - + 2 * b1 ** 3 * log_fact ** 3 - + b1 ** 3 * np.log(as_ref) ** 2 * (5 + 6 * log_fact) + -2 * b1**3 * np.log(as_ref) ** 3 + + 5 * b1**3 * log_fact**2 + + 2 * b1**3 * log_fact**3 + + b1**3 * np.log(as_ref) ** 2 * (5 + 6 * log_fact) + 2 * beta0 * b1 * log_fact - * (b2 + 2 * (b1 ** 2 - beta0 * b2) * lmu * as_ref) - - beta0 ** 2 + * (b2 + 2 * (b1**2 - beta0 * b2) * lmu * as_ref) + - beta0**2 * lmu * as_ref * ( -2 * b1 * b2 + 2 * beta0 * b3 - + (b1 ** 3 - 2 * beta0 * b1 * b2 + beta0 ** 2 * b3) + + (b1**3 - 2 * beta0 * b1 * b2 + beta0**2 * b3) * lmu * as_ref ) @@ -111,9 +111,9 @@ def as_expanded(order, as_ref, nf, scale_from, scale_to): * b1 * np.log(as_ref) * ( - 5 * b1 ** 2 * log_fact - + 3 * b1 ** 2 * log_fact ** 2 - + beta0 * (b2 + 2 * (b1 ** 2 - beta0 * b2) * lmu * as_ref) + 5 * b1**2 * log_fact + + 3 * b1**2 * log_fact**2 + + beta0 * (b2 + 2 * (b1**2 - beta0 * b2) * lmu * as_ref) ) ) ) diff --git a/tests/test_beta.py b/tests/test_beta.py index a0c757fa3..b6c61c167 100644 --- a/tests/test_beta.py +++ b/tests/test_beta.py @@ -44,7 +44,7 @@ def test_beta_3(): _flav_test(beta.beta_3) # from hep-ph/9706430 np.testing.assert_allclose( - beta.beta_3(5), 4 ** 4 * (11027.0 / 648.0 * zeta3 - 598391.0 / 373248.0) + beta.beta_3(5), 4**4 * (11027.0 / 648.0 * zeta3 - 598391.0 / 373248.0) ) diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index a58c5e60d..c4e06b7b6 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -240,11 +240,11 @@ def test_log_eval_N(self): assert_almost_equal(act_c, res_c) # Full -> \tilde p_0(N) = exp(-N)(exp(N)-1-N)/N^2 # MMa: Integrate[x^(n-1) (-Log[x]),{x,1/E,1}] - p0Nref_full = lambda N, lnx: ((np.exp(N) - 1 - N) / N ** 2) * np.exp( + p0Nref_full = lambda N, lnx: ((np.exp(N) - 1 - N) / N**2) * np.exp( -N * (lnx + 1) ) # partial = lower bound is neglected; - p0Nref_partial = lambda N, lnx: (1 / N ** 2) * np.exp(-N * lnx) + p0Nref_partial = lambda N, lnx: (1 / N**2) * np.exp(-N * lnx) p1N = inter_N[1] assert len(p1N.areas) == 1 p1_cs_ref = [1, 1] @@ -252,8 +252,8 @@ def test_log_eval_N(self): assert_almost_equal(act_c, res_c) # p_1(x) = 1+\ln(x) -> \tilde p_1(N) = (exp(-N)-1+N)/N^2 # MMa: Integrate[x^(n-1) (1+Log[x]),{x,1/E,1}] - p1Nref_full = lambda N, lnx: ((np.exp(-N) - 1 + N) / N ** 2) * np.exp(-N * lnx) - p1Nref_partial = lambda N, lnx: (1 / N - 1 / N ** 2) * np.exp(-N * lnx) + p1Nref_full = lambda N, lnx: ((np.exp(-N) - 1 + N) / N**2) * np.exp(-N * lnx) + p1Nref_partial = lambda N, lnx: (1 / N - 1 / N**2) * np.exp(-N * lnx) # iterate configurations for N in [1.0, 2.0, complex(1.0, 1.0)]: # check skip diff --git a/tests/test_strong_coupling.py b/tests/test_strong_coupling.py index 8ef173611..4b6e1b88c 100644 --- a/tests/test_strong_coupling.py +++ b/tests/test_strong_coupling.py @@ -166,10 +166,10 @@ def test_exact_LO(self): def benchmark_expanded_n3lo(self): """test N3LO - NNLO expansion with some reference value from Mathematica""" - Q2 = 100 ** 2 + Q2 = 100**2 # use a big alpha_s to enlarge the difference alphas_ref = 0.9 - scale_ref = 90 ** 2 + scale_ref = 90**2 m2c = 2 m2b = 25 m2t = 30625 From 0a9de58214e46f7c029af953e650b4d8fa618f01 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Wed, 16 Feb 2022 11:41:03 +0100 Subject: [PATCH 23/24] typos and fix acknowledge Luca Rottoli --- benchmarks/apfel_bench.py | 2 +- src/eko/strong_coupling.py | 3 +-- tests/benchmark_msbar_evolution.py | 2 +- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/benchmarks/apfel_bench.py b/benchmarks/apfel_bench.py index cf0f12a61..882b982ac 100644 --- a/benchmarks/apfel_bench.py +++ b/benchmarks/apfel_bench.py @@ -180,4 +180,4 @@ def benchmark_sv(self, pto): # obj.benchmark_plain(1) # obj.benchmark_sv(1) # obj.benchmark_kthr(2) - obj.benchmark_msbar(0) + obj.benchmark_msbar(2) diff --git a/src/eko/strong_coupling.py b/src/eko/strong_coupling.py index 6d8c9efae..3e977da7a 100644 --- a/src/eko/strong_coupling.py +++ b/src/eko/strong_coupling.py @@ -79,8 +79,7 @@ def as_expanded(order, as_ref, nf, scale_from, scale_to): + as_LO * (as_LO - as_ref) * (b2 - b1**2) + as_NLO * b1 * np.log(as_NLO / as_ref) ) - # N3LO expansion is taken from Luca Rottoli and it's simpler - # than the APFEL expansion which is more accurate + # N3LO expansion is taken from Luca Rottoli if order >= 3: b3 = beta_b(3, nf) log_fact = np.log(as_LO) diff --git a/tests/benchmark_msbar_evolution.py b/tests/benchmark_msbar_evolution.py index 0b24b4b47..63f16eb63 100644 --- a/tests/benchmark_msbar_evolution.py +++ b/tests/benchmark_msbar_evolution.py @@ -193,7 +193,7 @@ def benchmark_msbar_solution_kthr(self): With this test you can see that in EKO the solution value of mb is not affected by "kbThr", since mb is searched with an Nf=5 larger range. - While in Apfel this doesn't happend + While in Apfel this doesn't happen. """ theory_dict.update( { From 9263ce64008c04ceb63621c191483b761027c242 Mon Sep 17 00:00:00 2001 From: giacomomagni Date: Sun, 20 Feb 2022 16:07:40 +0100 Subject: [PATCH 24/24] revert lz4.frame in output --- src/eko/output.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/eko/output.py b/src/eko/output.py index 28c004926..45c2f387c 100644 --- a/src/eko/output.py +++ b/src/eko/output.py @@ -9,9 +9,9 @@ import tempfile import warnings +import lz4.frame import numpy as np import yaml -from lz4 import frame as lz4_frame from . import basis_rotation as br from . import interpolation, version @@ -308,7 +308,7 @@ def get_raw(self, binarize=True, skip_q2_grid=False): out["Q2grid"][q2][k] = float(v) continue if binarize: - out["Q2grid"][q2][k] = lz4_frame.compress(v.tobytes()) + out["Q2grid"][q2][k] = lz4.frame.compress(v.tobytes()) else: out["Q2grid"][q2][k] = v.tolist() else: @@ -393,7 +393,7 @@ def dump_tar(self, tarname): stream = io.BytesIO() np.save(stream, operator) stream.seek(0) - with lz4_frame.open( + with lz4.frame.open( (tmpdir / kind).with_suffix(".npy.lz4"), "wb" ) as fo: fo.write(stream.read()) @@ -436,7 +436,7 @@ def load_yaml(cls, stream, skip_q2_grid=False): elif isinstance(v, list): v = np.array(v) elif isinstance(v, bytes): - v = np.frombuffer(lz4_frame.decompress(v)) + v = np.frombuffer(lz4.frame.decompress(v)) v = v.reshape(len_tpids, len_tgrid, len_ipids, len_igrid) op[k] = v return cls(obj) @@ -500,7 +500,7 @@ def load_tar(cls, tarname): grids = {} for fp in innerdir.glob("*.npy.lz4"): - with lz4_frame.open(fp, "rb") as fd: + with lz4.frame.open(fp, "rb") as fd: stream = io.BytesIO(fd.read()) stream.seek(0) grids[pathlib.Path(fp.stem).stem] = np.load(stream)