Next runcards update

benchmarks/eko/benchmark_alphaem.py

@@ -35,9 +35,7 @@ def test_alphaQED_high(self):
             dict(
                 alphas=0.118,
                 alphaem=7.7553e-03,
-                scale=91.2,
-                num_flavs_ref=5,
-                max_num_flavs=5,
+                ref=(91.2, 5),
                 em_running=True,
             )
         )
@@ -75,9 +73,7 @@ def test_alphaQED_low(self):
             dict(
                 alphas=0.118,
                 alphaem=7.7553e-03,
-                scale=91.2,
-                num_flavs_ref=5,
-                max_num_flavs=5,
+                ref=(91.2, 5),
                 em_running=True,
             )
         )
@@ -124,9 +120,7 @@ def test_validphys(self):
             dict(
                 alphas=0.118,
                 alphaem=7.7553e-03,
-                scale=91.2,
-                num_flavs_ref=5,
-                max_num_flavs=5,
+                ref=(91.2, 5),
                 em_running=True,
             )
         )

benchmarks/eko/benchmark_evol_to_unity.py

@@ -19,16 +19,12 @@ def update_cards(theory: TheoryCard, operator: OperatorCard):
     theory.couplings = CouplingsInfo(
         alphas=0.35,
         alphaem=0.007496,
-        scale=float(np.sqrt(2)),
-        max_num_flavs=6,
-        num_flavs_ref=None,
+        ref=(float(np.sqrt(2)), 4),
     )
-    theory.heavy.num_flavs_init = 4
-    theory.heavy.intrinsic_flavors = [4, 5]
     theory.heavy.masses.c.value = 1.0
     theory.heavy.masses.b.value = 4.75
     theory.heavy.masses.t.value = 173.0
-    operator.mu0 = float(np.sqrt(2))
+    operator.init = (float(np.sqrt(2)), 4)
     operator.mugrid = [(10, 5)]
     operator.xgrid = XGrid(np.linspace(1e-1, 1, 30))
     operator.configs.interpolation_polynomial_degree = 1

benchmarks/eko/benchmark_inverse_matching.py

@@ -20,14 +20,9 @@
     couplings=dict(
         alphas=0.118,
         alphaem=0.007496252,
-        scale=91.2,
-        num_flavs_ref=5,
-        max_num_flavs=6,
+        ref=(91.2, 5),
     ),
     heavy=dict(
-        num_flavs_init=4,
-        num_flavs_max_pdf=6,
-        intrinsic_flavors=[],
         masses=[ReferenceRunning([mq, np.nan]) for mq in (MC, 4.92, 172.5)],
         masses_scheme="POLE",
         matching_ratios=[1.0, 1.0, np.inf],
@@ -40,7 +35,7 @@
 
 # operator settings
 op_raw = dict(
-    mu0=1.65,
+    init=(1.65, 4),
     xgrid=[0.0001, 0.001, 0.01, 0.1, 1],
     mugrid=[(MC, 3), (MC, 4)],
     configs=dict(

benchmarks/eko/benchmark_msbar_evolution.py

@@ -21,9 +21,8 @@
 
 def update_theory(theory: TheoryCard):
     theory.order = (3, 0)
-    theory.couplings.scale = 91
+    theory.couplings.ref = (91, 5)
     theory.couplings.alphaem = 0.007496
-    theory.couplings.num_flavs_ref = 5
     theory.heavy.masses_scheme = QuarkMassScheme.MSBAR
     theory.heavy.masses.c = QuarkMassRef([1.5, 18])
     theory.heavy.masses.b = QuarkMassRef([4.1, 20])
@@ -150,7 +149,7 @@ def benchmark_APFEL_msbar_evolution(
                     apfel.SetTheory("QCD")
                     apfel.SetPerturbativeOrder(order - 1)
                     apfel.SetAlphaEvolution(method)
-                    apfel.SetAlphaQCDRef(coupl.alphas, coupl.scale)
+                    apfel.SetAlphaQCDRef(coupl.alphas, coupl.ref[0])
                     apfel.SetVFNS()
                     apfel.SetMSbarMasses(
                         qmasses.c.value, qmasses.b.value, qmasses.t.value
@@ -218,7 +217,7 @@ def benchmark_APFEL_msbar_solution(
                 apfel.SetTheory("QCD")
                 apfel.SetPerturbativeOrder(order - 1)
                 apfel.SetAlphaEvolution("exact")
-                apfel.SetAlphaQCDRef(coupl.alphas, coupl.scale)
+                apfel.SetAlphaQCDRef(coupl.alphas, coupl.ref[0])
                 apfel.SetVFNS()
                 apfel.SetMSbarMasses(qmasses.c.value, qmasses.b.value, qmasses.t.value)
                 apfel.SetMassScaleReference(

benchmarks/eko/benchmark_strong_coupling.py

@@ -8,7 +8,7 @@
 from eko.couplings import Couplings
 from eko.io.runcards import TheoryCard
 from eko.quantities.couplings import CouplingEvolutionMethod, CouplingsInfo
-from eko.quantities.heavy_quarks import QuarkMassScheme
+from eko.quantities.heavy_quarks import FlavorsNumber, QuarkMassScheme
 
 # try to load LHAPDF - if not available, we'll use the cached values
 try:
@@ -35,16 +35,11 @@
     use_PEGASUS = False
 
 
-def ref_couplings(
-    ref_values,
-    ref_scale: float,
-) -> CouplingsInfo:
+def ref_couplings(ref_values, ref_scale: float, ref_nf: FlavorsNumber) -> CouplingsInfo:
     return CouplingsInfo(
         alphas=ref_values[0],
         alphaem=ref_values[1],
-        scale=ref_scale,
-        max_num_flavs=6,
-        num_flavs_ref=None,
+        ref=(ref_scale, ref_nf),
     )
 
 
@@ -54,11 +49,11 @@ def test_a_s(self):
         """Tests the value of alpha_s (for now only at LO)
         for a given set of parameters
         """
-        # TODO @JCM: we need a source for this!
+        # source: JCM
         known_val = 0.0091807954
         ref_mu2 = 90
         ask_q2 = 125
-        ref = ref_couplings([0.1181, 0.007496], np.sqrt(ref_mu2))
+        ref = ref_couplings([0.1181, 0.007496], np.sqrt(ref_mu2), 5)
         as_FFNS_LO = Couplings(
             couplings=ref,
             order=(1, 0),
@@ -80,7 +75,7 @@ def benchmark_LHA_paper(self):
         # LO - FFNS
         # note that the LO-FFNS value reported in :cite:`Giele:2002hx`
         # was corrected in :cite:`Dittmar:2005ed`
-        coupling_ref = ref_couplings([0.35, 0.007496], np.sqrt(2))
+        coupling_ref = ref_couplings([0.35, 0.007496], np.sqrt(2), 4)
         as_FFNS_LO = Couplings(
             couplings=coupling_ref,
             order=(1, 0),
@@ -141,7 +136,7 @@ def benchmark_APFEL_ffns(self):
         threshold_holder = matchings.Atlas.ffns(nf, 0.0)
         for order in [1, 2, 3]:
             as_FFNS = Couplings(
-                couplings=ref_couplings(coupling_ref, scale_ref),
+                couplings=ref_couplings(coupling_ref, scale_ref, nf),
                 order=(order, 0),
                 method=CouplingEvolutionMethod.EXPANDED,
                 masses=threshold_holder.walls[1:-1],
@@ -212,8 +207,7 @@ def benchmark_pegasus_ffns(self):
         }
         # collect my values
         threshold_holder = matchings.Atlas.ffns(nf, 0.0)
-        couplings = ref_couplings(coupling_ref, scale_ref)
-        couplings.max_num_flavs = 4
+        couplings = ref_couplings(coupling_ref, scale_ref, nf)
         for order in [1, 2, 3, 4]:
             as_FFNS = Couplings(
                 couplings=couplings,
@@ -259,6 +253,7 @@ def benchmark_APFEL_vfns(self):
         Q2s = [1, 2**2, 3**2, 90**2, 100**2]
         coupling_ref = np.array([0.118, 0.007496])
         scale_ref = 91.0
+        nf_ref = 5
         threshold_list = np.power([2, 4, 175], 2)
         apfel_vals_dict = {
             1: np.array(
@@ -292,7 +287,7 @@ def benchmark_APFEL_vfns(self):
         # collect my values
         for order in [1, 2, 3]:
             as_VFNS = Couplings(
-                couplings=ref_couplings(coupling_ref, scale_ref),
+                couplings=ref_couplings(coupling_ref, scale_ref, nf_ref),
                 order=(order, 0),
                 method=CouplingEvolutionMethod.EXPANDED,
                 masses=threshold_list,
@@ -339,8 +334,9 @@ def benchmark_APFEL_vfns_fact_to_ren(self):
         ]
         coupling_ref = np.array([0.118, 0.007496])
         scale_ref = 91.0
+        nf_refs = (5, 6)
         fact_to_ren_lin_list = [0.567, 2.34]
-        threshold_list = np.power([2, 2 * 4, 2 * 92], 2)
+        masses_list = np.power([2, 2 * 4, 2 * 92], 2)
         apfel_vals_dict_list = [
             {
                 1: np.array(
@@ -431,16 +427,16 @@ def benchmark_APFEL_vfns_fact_to_ren(self):
                 ),
             },
         ]
-        for fact_to_ren_lin, apfel_vals_dict in zip(
-            fact_to_ren_lin_list, apfel_vals_dict_list
+        for fact_to_ren_lin, apfel_vals_dict, nf_ref in zip(
+            fact_to_ren_lin_list, apfel_vals_dict_list, nf_refs
         ):
             # collect my values
             for order in [1, 2, 3]:
                 as_VFNS = Couplings(
-                    couplings=ref_couplings(coupling_ref, scale_ref),
+                    couplings=ref_couplings(coupling_ref, scale_ref, nf_ref),
                     order=(order, 0),
                     method=CouplingEvolutionMethod.EXACT,
-                    masses=threshold_list.tolist(),
+                    masses=masses_list.tolist(),
                     hqm_scheme=QuarkMassScheme.POLE,
                     thresholds_ratios=np.array([1.0, 1.0, 1.0]) / fact_to_ren_lin**2,
                 )
@@ -457,7 +453,7 @@ def benchmark_APFEL_vfns_fact_to_ren(self):
                     apfel.SetAlphaEvolution("exact")
                     apfel.SetAlphaQCDRef(coupling_ref[0], scale_ref)
                     apfel.SetVFNS()
-                    apfel.SetPoleMasses(*np.sqrt(threshold_list))
+                    apfel.SetPoleMasses(*np.sqrt(masses_list))
                     apfel.SetRenFacRatio(1.0 / fact_to_ren_lin)
                     apfel.InitializeAPFEL()
                     # collect a_s
@@ -468,12 +464,18 @@ def benchmark_APFEL_vfns_fact_to_ren(self):
                         )
                     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=2.5e-5)
+                np.testing.assert_allclose(
+                    apfel_vals,
+                    np.array(my_vals),
+                    rtol=2.5e-5,
+                    err_msg=f"{order=},{fact_to_ren_lin=}",
+                )
 
     def benchmark_APFEL_vfns_threshold(self):
         Q2s = np.power([30, 96, 150], 2)
         coupling_ref = np.array([0.118, 0.007496])
         scale_ref = 91.0
+        nf_ref = 4
         threshold_list = np.power([30, 95, 240], 2)
         thresholds_ratios = [2.34**2, 1.0**2, 0.5**2]
         apfel_vals_dict = {
@@ -487,7 +489,7 @@ def benchmark_APFEL_vfns_threshold(self):
         # collect my values
         for order in [2, 3]:
             as_VFNS = Couplings(
-                couplings=ref_couplings(coupling_ref, scale_ref),
+                couplings=ref_couplings(coupling_ref, scale_ref, nf_ref),
                 order=(order, 0),
                 method=CouplingEvolutionMethod.EXPANDED,
                 masses=threshold_list.tolist(),
@@ -496,7 +498,6 @@ def benchmark_APFEL_vfns_threshold(self):
             )
             my_vals = []
             for Q2 in Q2s:
-                print(Q2)
                 my_vals.append(as_VFNS.a(Q2)[0])
             # get APFEL numbers - if available else use cache
             apfel_vals = apfel_vals_dict[order]
@@ -516,7 +517,6 @@ def benchmark_APFEL_vfns_threshold(self):
                 apfel_vals_cur = []
                 for Q2 in Q2s:
                     apfel_vals_cur.append(apfel.AlphaQCD(np.sqrt(Q2)) / (4.0 * np.pi))
-                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))
@@ -525,6 +525,7 @@ def benchmark_APFEL_vfns_msbar(self):
         Q2s = np.power([3, 96, 150], 2)
         coupling_ref = np.array([0.118, 0.007496])
         scale_ref = 91.0
+        nf_ref = 5
         Q2m = np.power([2.0, 4.0, 175], 2)
         m2s = np.power((1.4, 4.0, 175), 2)
         apfel_vals_dict = {
@@ -538,7 +539,7 @@ def benchmark_APFEL_vfns_msbar(self):
         # collect my values
         for order in [2, 3]:
             as_VFNS = Couplings(
-                couplings=ref_couplings(coupling_ref, scale_ref),
+                couplings=ref_couplings(coupling_ref, scale_ref, nf_ref),
                 order=(order, 0),
                 method=CouplingEvolutionMethod.EXPANDED,
                 masses=m2s.tolist(),
@@ -586,7 +587,7 @@ def benchmark_lhapdf_ffns_lo(self):
         # collect my values
         threshold_holder = matchings.Atlas.ffns(nf, 0.0)
         as_FFNS_LO = Couplings(
-            couplings=ref_couplings(coupling_ref, scale_ref),
+            couplings=ref_couplings(coupling_ref, scale_ref, nf),
             order=(1, 0),
             method=CouplingEvolutionMethod.EXACT,
             masses=threshold_holder.walls[1:-1],
@@ -629,8 +630,9 @@ def benchmark_apfel_exact(self):
         Q2s = [1e1, 1e2, 1e3, 1e4]
         coupling_ref = np.array([0.118, 0.007496])
         scale_ref = 90
+        nf = 3
         # collect my values
-        threshold_holder = matchings.Atlas.ffns(3, 0.0)
+        threshold_holder = matchings.Atlas.ffns(nf, 0.0)
         # LHAPDF cache
         apfel_vals_dict = {
             1: np.array(
@@ -660,7 +662,7 @@ def benchmark_apfel_exact(self):
         }
         for order in range(1, 3 + 1):
             sc = Couplings(
-                couplings=ref_couplings(coupling_ref, scale_ref),
+                couplings=ref_couplings(coupling_ref, scale_ref, nf),
                 order=(order, 0),
                 method=CouplingEvolutionMethod.EXACT,
                 masses=threshold_holder.walls[1:-1],
@@ -696,8 +698,9 @@ def benchmark_lhapdf_exact(self):
         Q2s = [1e1, 1e2, 1e3, 1e4]
         coupling_ref = np.array([0.118, 0.007496])
         scale_ref = 90
+        nf = 3
         # collect my values
-        threshold_holder = matchings.Atlas.ffns(3, 0.0)
+        threshold_holder = matchings.Atlas.ffns(nf, 0.0)
         # LHAPDF cache
         lhapdf_vals_dict = {
             1: np.array(
@@ -735,7 +738,7 @@ def benchmark_lhapdf_exact(self):
         }
         for order in range(1, 4 + 1):
             sc = Couplings(
-                couplings=ref_couplings(coupling_ref, scale_ref),
+                couplings=ref_couplings(coupling_ref, scale_ref, nf),
                 order=(order, 0),
                 method=CouplingEvolutionMethod.EXACT,
                 masses=threshold_holder.walls[1:-1],
@@ -766,6 +769,7 @@ def benchmark_lhapdf_zmvfns_lo(self):
         Q2s = [1, 1e1, 1e2, 1e3, 1e4]
         coupling_ref = np.array([0.118, 0.007496])
         scale_ref = 900
+        nf_ref = 5
         m2c = 2.0
         m2b = 25.0
         m2t = 1500.0
@@ -785,7 +789,7 @@ def benchmark_lhapdf_zmvfns_lo(self):
 
         # collect my values
         as_VFNS_LO = Couplings(
-            couplings=ref_couplings(coupling_ref, np.sqrt(scale_ref)),
+            couplings=ref_couplings(coupling_ref, np.sqrt(scale_ref), nf_ref),
             order=(1, 0),
             method=CouplingEvolutionMethod.EXACT,
             masses=threshold_list,
@@ -836,10 +840,8 @@ def benchmark_APFEL_fact_to_ren_lha_settings(self, theory_card: TheoryCard):
         theory = theory_card
         theory.order = (3, 0)
         theory.couplings.alphas = 0.35
-        theory.couplings.scale = float(np.sqrt(2))
         theory.couplings.alphaem = 0.007496
-        theory.couplings.num_flavs_ref = 4
-        theory.heavy.num_flavs_init = 3
+        theory.couplings.ref = (float(np.sqrt(2)), 4)
         theory.xif = np.sqrt(1.0 / 2.0)
         theory.heavy.masses.c.value = np.sqrt(2.0)
         theory.heavy.masses.b.value = 4.5
@@ -881,7 +883,7 @@ def benchmark_APFEL_fact_to_ren_lha_settings(self, theory_card: TheoryCard):
                 apfel.SetTheory("QCD")
                 apfel.SetPerturbativeOrder(theory.order[0] - 1)
                 apfel.SetAlphaEvolution("exact")
-                apfel.SetAlphaQCDRef(theory.couplings.alphas, theory.couplings.scale)
+                apfel.SetAlphaQCDRef(theory.couplings.alphas, theory.couplings.ref[0])
                 apfel.SetVFNS()
                 apfel.SetPoleMasses(qmasses.c.value, qmasses.b.value, qmasses.t.value)
                 apfel.SetMassMatchingScales(