Pass nl as an argument in qed beta functions

benchmarks/eko/benchmark_alphaem.py

@@ -0,0 +1,144 @@
+"""This module benchmarks alpha_em against alphaQED23 and validphys.
+
+alphaQED23 can be obtained from http://www-com.physik.hu-berlin.de/~fjeger/software.html .
+"""
+
+import numpy as np
+import pytest
+
+from eko.couplings import Couplings
+from eko.quantities.couplings import CouplingEvolutionMethod, CouplingsInfo
+from eko.quantities.heavy_quarks import QuarkMassScheme
+
+
+@pytest.mark.isolated
+class BenchmarkCouplings:
+    def test_alphaQED_high(self):
+        """testing aginst alphaQED23 for high Q values"""
+        alphaQED23 = np.array(
+            [
+                0.0075683,
+                0.0075867,
+                0.0076054,
+                0.0076245,
+                0.0076437,
+                0.0076631,
+                0.0076827,
+                0.0077024,
+                0.0077222,
+                0.0077421,
+                0.0077621,
+            ]
+        )
+        qvec = np.geomspace(10, 100, 11)
+        couplings = CouplingsInfo.from_dict(
+            dict(
+                alphas=0.118,
+                alphaem=7.7553e-03,
+                scale=91.2,
+                num_flavs_ref=5,
+                max_num_flavs=5,
+                em_running=True,
+            )
+        )
+        eko_alpha = Couplings(
+            couplings,
+            (3, 2),
+            method=CouplingEvolutionMethod.EXACT,
+            masses=[m**2 for m in [1.51, 4.92, 172.5]],
+            hqm_scheme=QuarkMassScheme.POLE,
+            thresholds_ratios=[1.0, 1.0, np.inf],
+        )
+        alpha_eko = np.array([eko_alpha.a_em(q**2) * 4 * np.pi for q in qvec])
+
+        np.testing.assert_allclose(alphaQED23, alpha_eko, rtol=1.8e-4)
+
+    def test_alphaQED_low(self):
+        """testing aginst alphaQED23 for low Q values: they are close but not identical"""
+        alphaQED23 = np.array(
+            [
+                0.0074192,
+                0.007431,
+                0.0074434,
+                0.0074565,
+                0.0074702,
+                0.0074847,
+                0.0075,
+                0.0075161,
+                0.0075329,
+                0.0075503,
+                0.0075683,
+            ]
+        )
+        qvec = np.geomspace(1, 10, 11)
+        couplings = CouplingsInfo.from_dict(
+            dict(
+                alphas=0.118,
+                alphaem=7.7553e-03,
+                scale=91.2,
+                num_flavs_ref=5,
+                max_num_flavs=5,
+                em_running=True,
+            )
+        )
+        eko_alpha = Couplings(
+            couplings,
+            (3, 2),
+            method=CouplingEvolutionMethod.EXACT,
+            masses=[m**2 for m in [1.51, 4.92, 172.5]],
+            hqm_scheme=QuarkMassScheme.POLE,
+            thresholds_ratios=[1.0, 1.0, np.inf],
+        )
+        alpha_eko = np.array([eko_alpha.a_em(q**2) * 4 * np.pi for q in qvec])
+
+        np.testing.assert_allclose(alphaQED23, alpha_eko, rtol=3.2e-3)
+
+    def test_validphys(self):
+        """testing aginst validphys"""
+        alpha_vp = np.array(
+            [
+                0.007408255774054356,
+                0.007425240094018394,
+                0.007449051094996458,
+                0.007476301027742958,
+                0.007503751810862984,
+                0.007532299008699658,
+                0.007561621958607614,
+                0.007591174885612722,
+                0.007620960508577136,
+                0.00765098158940323,
+                0.007681240933888789,
+                0.007711741392602655,
+                0.007742485861781425,
+                0.007773477284247778,
+                0.007804718650351058,
+                0.007836212998930739,
+                0.00786796341830342,
+                0.007899973047274033,
+                0.007932245076171957,
+                0.00796478274791273,
+            ]
+        )
+        qvec = np.geomspace(1, 1000, 20)
+        couplings = CouplingsInfo.from_dict(
+            dict(
+                alphas=0.118,
+                alphaem=7.7553e-03,
+                scale=91.2,
+                num_flavs_ref=5,
+                max_num_flavs=5,
+                em_running=True,
+            )
+        )
+        eko_alpha = Couplings(
+            couplings,
+            (3, 2),
+            method=CouplingEvolutionMethod.EXACT,
+            masses=[m**2 for m in [1.51, 4.92, 172.5]],
+            hqm_scheme=QuarkMassScheme.POLE,
+            thresholds_ratios=[1.0, 1.0, np.inf],
+        )
+        eko_alpha.decoupled_running = True
+        alpha_eko = np.array([eko_alpha.a_em(q**2) * 4 * np.pi for q in qvec])
+
+        np.testing.assert_allclose(alpha_vp, alpha_eko, rtol=5e-6)

src/eko/beta.py

@@ -32,7 +32,7 @@ def beta_qcd_as2(nf):
 
 
 @nb.njit(cache=True)
-def beta_qed_aem2(nf):
+def beta_qed_aem2(nf, nl):
     r"""Compute the first coefficient of the QED beta function.
 
     Implements Eq. (7) of :cite:`Surguladze:1996hx`.
@@ -41,6 +41,8 @@ def beta_qed_aem2(nf):
     ----------
         nf : int
             number of active flavors
+        nl : int
+            number of leptons partecipating to alphaem running
 
     Returns
     -------
@@ -50,7 +52,6 @@ def beta_qed_aem2(nf):
     """
     nu = constants.uplike_flavors(nf)
     nd = nf - nu
-    nl = 3  # TODO : pass nl as an argument??
     beta_qed_aem2 = (
         -4.0 / 3 * (nl + constants.NC * (nu * constants.eu2 + nd * constants.ed2))
     )
@@ -83,7 +84,7 @@ def beta_qcd_as3(nf):
 
 
 @nb.njit(cache=True)
-def beta_qed_aem3(nf):
+def beta_qed_aem3(nf, nl):
     r"""Compute the second coefficient of the QED beta function.
 
     Implements Eq. (7) of :cite:`Surguladze:1996hx`.
@@ -92,6 +93,8 @@ def beta_qed_aem3(nf):
     ----------
         nf : int
             number of active flavors
+        nl : int
+            number of leptons partecipating to alphaem running
 
     Returns
     -------
@@ -101,7 +104,6 @@ def beta_qed_aem3(nf):
     """
     nu = constants.uplike_flavors(nf)
     nd = nf - nu
-    nl = 3  # TODO : pass nl as an argument??
     beta_qed_aem3 = -4.0 * (
         nl + constants.NC * (nu * constants.eu2**2 + nd * constants.ed2**2)
     )
@@ -246,7 +248,7 @@ def beta_qcd(k, nf):
 
 
 @nb.njit(cache=True)
-def beta_qed(k, nf):
+def beta_qed(k, nf, nl):
     r"""Compute value of a beta_qed coefficients.
 
     Parameters
@@ -255,6 +257,8 @@ def beta_qed(k, nf):
             perturbative order
         nf : int
             number of active flavors
+        nl : int
+            number of leptons partecipating to alphaem running
 
     Returns
     -------
@@ -264,9 +268,9 @@ def beta_qed(k, nf):
     """
     beta_ = 0
     if k == (0, 2):
-        beta_ = beta_qed_aem2(nf)
+        beta_ = beta_qed_aem2(nf, nl)
     elif k == (0, 3):
-        beta_ = beta_qed_aem3(nf)
+        beta_ = beta_qed_aem3(nf, nl)
     elif k == (1, 2):
         beta_ = beta_qed_aem2as1(nf)
     else:
@@ -295,7 +299,7 @@ def b_qcd(k, nf):
 
 
 @nb.njit(cache=True)
-def b_qed(k, nf):
+def b_qed(k, nf, nl):
     r"""Compute b_qed coefficient.
 
     Parameters
@@ -304,11 +308,13 @@ def b_qed(k, nf):
             perturbative order
         nf : int
             number of active flavors
+        nl : int
+            number of leptons partecipating to alphaem running
 
     Returns
     -------
         b_qed : float
             b_qed_k(nf)
 
     """
-    return beta_qed(k, nf) / beta_qed((0, 2), nf)
+    return beta_qed(k, nf, nl) / beta_qed((0, 2), nf, nl)

src/eko/constants.py

@@ -25,6 +25,9 @@
 Defaults to :math:`\frac{N_C^2-1}{2N_C} = 4/3`.
 """
 
+MTAU = 1.777
+"""Mass of the tau."""
+
 eu2 = 4.0 / 9
 """Up quarks charge squared."""
 
@@ -84,7 +87,7 @@ def uplike_flavors(nf):
     nu : int
 
     """
-    if nf not in range(2, 6 + 1):
+    if nf > 6:
         raise NotImplementedError("Selected nf is not implemented")
     nu = nf // 2
     return nu