diff --git a/benchmarks/eko/benchmark_alphaem.py b/benchmarks/eko/benchmark_alphaem.py index 16e56ac4a..9cc05536c 100644 --- a/benchmarks/eko/benchmark_alphaem.py +++ b/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, ) ) diff --git a/benchmarks/eko/benchmark_evol_to_unity.py b/benchmarks/eko/benchmark_evol_to_unity.py index 41bc4b4d0..c466d00d8 100644 --- a/benchmarks/eko/benchmark_evol_to_unity.py +++ b/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 diff --git a/benchmarks/eko/benchmark_inverse_matching.py b/benchmarks/eko/benchmark_inverse_matching.py index d9c16b0b2..d63b776f0 100644 --- a/benchmarks/eko/benchmark_inverse_matching.py +++ b/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( diff --git a/benchmarks/eko/benchmark_msbar_evolution.py b/benchmarks/eko/benchmark_msbar_evolution.py index b9e3bd252..49dc0e952 100644 --- a/benchmarks/eko/benchmark_msbar_evolution.py +++ b/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( diff --git a/benchmarks/eko/benchmark_strong_coupling.py b/benchmarks/eko/benchmark_strong_coupling.py index d13a264bd..d60bd28b0 100644 --- a/benchmarks/eko/benchmark_strong_coupling.py +++ b/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( diff --git a/benchmarks/ekobox/benchmark_evol_pdf.py b/benchmarks/ekobox/benchmark_evol_pdf.py index 3a9dc0f4a..a0ddcaac4 100644 --- a/benchmarks/ekobox/benchmark_evol_pdf.py +++ b/benchmarks/ekobox/benchmark_evol_pdf.py @@ -21,15 +21,13 @@ def benchmark_evolve_single_member( theory = theory_card theory.order = (1, 0) theory.couplings.alphas = 0.118000 - theory.couplings.scale = 91.1876 + theory.couplings.ref = (91.1876, 5) theory.couplings.alphaem = 0.007496 - theory.couplings.max_num_flavs = 3 - theory.heavy.num_flavs_max_pdf = 3 theory.heavy.masses.c.value = 1.3 theory.heavy.masses.b.value = 4.75 theory.heavy.masses.t.value = 172 op = operator_card - op.mu0 = 5.0 + op.init = (5.0, 4) op.mugrid = mugrid # lhapdf import (maybe i have to dump with a x*), do plots) with lhapdf_path(test_pdf): @@ -49,7 +47,7 @@ def benchmark_evolve_single_member( ev_pdf = lhapdf.mkPDF("EvPDF", 0) assert info["XMin"] == op.xgrid.raw[0] assert info["SetDesc"] == "MyEvolvedPDF" - assert info["MZ"] == theory.couplings.scale + assert info["MZ"] == theory.couplings.ref[0] assert info["Debug"] == "Debug" xgrid = op.xgrid.raw for idx, mu2 in enumerate(op.mu2grid): @@ -76,7 +74,7 @@ def benchmark_evolve_more_members( theory = theory_card theory.order = (1, 0) op = operator_card - op.mu0 = 1.0 + op.init = (1.0, 3) op.mugrid = [(10.0, 5), (100.0, 5)] op.xgrid = XGrid([1e-7, 0.01, 0.1, 0.2, 0.3]) with lhapdf_path(test_pdf): diff --git a/benchmarks/lha_paper_bench.py b/benchmarks/lha_paper_bench.py index 137d8b3d7..bdc48ce59 100644 --- a/benchmarks/lha_paper_bench.py +++ b/benchmarks/lha_paper_bench.py @@ -13,19 +13,15 @@ register(__file__) +_sqrt2 = float(np.sqrt(2.0)) + base_theory = { "ModEv": "EXA", - "Q0": np.sqrt( - 2.0 - ), # Eq. (30) :cite:`Giele:2002hx`, Eq. (4.53) :cite:`Dittmar:2005ed` - "mc": np.sqrt( - 2.0 - ), # Eq. (34) :cite:`Giele:2002hx`, Eq. (4.56) :cite:`Dittmar:2005ed` + "Q0": _sqrt2, # Eq. (30) :cite:`Giele:2002hx`, Eq. (4.53) :cite:`Dittmar:2005ed` + "mc": _sqrt2, # Eq. (34) :cite:`Giele:2002hx`, Eq. (4.56) :cite:`Dittmar:2005ed` "mb": 4.5, "mt": 175, - "Qref": np.sqrt( - 2.0 - ), # Eq. (32) :cite:`Giele:2002hx`,Eq. (4.53) :cite:`Dittmar:2005ed` + "Qref": _sqrt2, # Eq. (32) :cite:`Giele:2002hx`,Eq. (4.53) :cite:`Dittmar:2005ed` "alphas": 0.35, # Eq. (4.55) :cite:`Dittmar:2005ed` "alphaqed": 0.007496, "QED": 0, @@ -79,11 +75,11 @@ def sv_theories(self, pto): """ low = self.theory.copy() low["PTO"] = pto - low["XIF"] = np.sqrt(1.0 / 2.0) + low["XIF"] = 1.0 / _sqrt2 low["ModSV"] = "exponentiated" high = self.theory.copy() high["PTO"] = pto - high["XIF"] = np.sqrt(2.0) + high["XIF"] = _sqrt2 high["ModSV"] = "exponentiated" return [high, low] @@ -302,7 +298,7 @@ def benchmark_sv(self, pto): sv_theory["kcThr"] = 1.0 + 1e-15 sv_theory["nfref"] = 4 sv_theory["EScaleVar"] = 0 - low["XIR"] = np.sqrt(2.0) - high["XIR"] = np.sqrt(0.5) + low["XIR"] = _sqrt2 + high["XIR"] = 1.0 / _sqrt2 self.run_lha([low, high]) diff --git a/doc/source/overview/tutorials/alpha_s.ipynb b/doc/source/overview/tutorials/alpha_s.ipynb index a9d75c0ca..f9cc5a240 100644 --- a/doc/source/overview/tutorials/alpha_s.ipynb +++ b/doc/source/overview/tutorials/alpha_s.ipynb @@ -39,12 +39,10 @@ "from eko.quantities.heavy_quarks import MatchingScales, QuarkMassScheme\n", "\n", "# set the (alpha_s, alpha_em) reference values\n", - "couplings_ref = CouplingsInfo(\n", - " alphas=0.118, alphaem=0.007496252, scale=91.0, num_flavs_ref=None, max_num_flavs=5\n", - ")\n", + "couplings_ref = CouplingsInfo(alphas=0.118, alphaem=0.007496252, ref=(91.0, 5))\n", "\n", "# set heavy quark masses and their threshold ratios\n", - "heavy_quark_masses = np.power([1.51, 4.92, 172.0],2)\n", + "heavy_quark_masses = np.power([1.51, 4.92, 172.0], 2)\n", "thresholds_ratios = np.array([1.0, 1.0, 1.0])\n", "\n", "# set (QCD,QED) perturbative order\n", @@ -88,9 +86,9 @@ } ], "source": [ - "target_scale = 10.0 ** 2\n", + "target_scale = 10.0**2\n", "a_s = sc.a_s(target_scale)\n", - "print(\"The value of alpha_s at Q^2=100 GeV^2 is: \", 4. * np.pi * a_s)" + "print(\"The value of alpha_s at Q^2=100 GeV^2 is: \", 4.0 * np.pi * a_s)" ] }, { @@ -105,13 +103,10 @@ } ], "metadata": { - "interpreter": { - "hash": "0a84ba3ac8703c04e87bc503a7d00188dfd591ad56130da93c406115a1e4a408" - }, "kernelspec": { - "display_name": "eko-KkPVjVhh-py3.10", + "display_name": "Python 3 (ipykernel)", "language": "python", - "name": "eko-kkpvjvhh-py3.10" + "name": "python3" }, "language_info": { "codemirror_mode": { diff --git a/doc/source/overview/tutorials/dglap.ipynb b/doc/source/overview/tutorials/dglap.ipynb index 0e271a00c..98cb4a6e1 100644 --- a/doc/source/overview/tutorials/dglap.ipynb +++ b/doc/source/overview/tutorials/dglap.ipynb @@ -68,7 +68,7 @@ "th_card = example.theory()\n", "op_card = example.operator()\n", "# here we replace the grid with a very minimal one, to speed up the example\n", - "op_card.xgrid = [1e-3, 1e-2, 1e-1, 5e-1, 1.]" + "op_card.xgrid = [1e-3, 1e-2, 1e-1, 5e-1, 1.0]" ] }, { @@ -91,14 +91,9 @@ "{'order': [1, 0],\n", " 'couplings': {'alphas': 0.118,\n", " 'alphaem': 0.007496252,\n", - " 'scale': 91.2,\n", - " 'max_num_flavs': 6,\n", - " 'num_flavs_ref': 5,\n", + " 'ref': [91.2, 5],\n", " 'em_running': False},\n", - " 'heavy': {'num_flavs_init': 4,\n", - " 'num_flavs_max_pdf': 6,\n", - " 'intrinsic_flavors': [4],\n", - " 'masses': [[2.0, nan], [4.5, nan], [173.07, nan]],\n", + " 'heavy': {'masses': [[2.0, nan], [4.5, nan], [173.07, nan]],\n", " 'masses_scheme': 'pole',\n", " 'matching_ratios': [1.0, 1.0, 1.0]},\n", " 'xif': 1.0,\n", @@ -123,7 +118,7 @@ { "data": { "text/plain": [ - "{'mu0': 1.65,\n", + "{'init': [1.65, 4],\n", " 'mugrid': [[100.0, 5]],\n", " 'xgrid': [0.001, 0.01, 0.1, 0.5, 1.0],\n", " 'configs': {'evolution_method': 'iterate-exact',\n", @@ -174,7 +169,7 @@ "id": "880aadcf-8f87-4918-a0bc-09581d0d3579", "metadata": {}, "source": [ - "The actual result is a complicate EKO object, we will discuss it in a separate tutorial.\n", + "The actual result is a complicate EKO object, which we will discuss it in a separate tutorial.\n", "\n", "You have just run your first DGLAP calculation!" ] diff --git a/doc/source/overview/tutorials/output.ipynb b/doc/source/overview/tutorials/output.ipynb index 47345ac35..1cfe9f071 100644 --- a/doc/source/overview/tutorials/output.ipynb +++ b/doc/source/overview/tutorials/output.ipynb @@ -34,7 +34,7 @@ "id": "2f8f0666-c6ab-40f6-86f5-15773f205b51", "metadata": {}, "source": [ - "We can access the operator, by using the `open` method (similar to python's `open`):" + "We can access the operator, by using the `read` method:" ] }, { @@ -76,8 +76,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "TheoryCard(order=(1, 0), couplings=CouplingsInfo(alphas=0.118, alphaem=0.007496252, scale=91.2, max_num_flavs=6, num_flavs_ref=5, em_running=False), heavy=HeavyInfo(num_flavs_init=4, num_flavs_max_pdf=6, intrinsic_flavors=[4, 5, 6], masses=[[2.0, nan], [4.5, nan], [173.07, nan]], masses_scheme=, matching_ratios=[1.0, 1.0, 1.0]), xif=1.0, n3lo_ad_variation=(0, 0, 0, 0))\n", - "OperatorCard(mu0=1.65, mugrid=[(100.0, 5)], xgrid=, configs=Configs(evolution_method=, ev_op_max_order=(10, 0), ev_op_iterations=10, scvar_method=None, inversion_method=None, interpolation_polynomial_degree=4, interpolation_is_log=True, polarized=False, time_like=False, n_integration_cores=0), debug=Debug(skip_singlet=False, skip_non_singlet=False), eko_version='0.0.0')\n" + "TheoryCard(order=(1, 0), couplings=CouplingsInfo(alphas=0.118, alphaem=0.007496252, ref=(91.2, 5), em_running=False), heavy=HeavyInfo(masses=[[2.0, nan], [4.5, nan], [173.07, nan]], masses_scheme=, matching_ratios=[1.0, 1.0, 1.0]), xif=1.0, n3lo_ad_variation=(0, 0, 0, 0))\n", + "OperatorCard(init=(1.65, 4), mugrid=[(100.0, 5)], xgrid=, configs=Configs(evolution_method=, ev_op_max_order=(10, 0), ev_op_iterations=10, scvar_method=None, inversion_method=None, interpolation_polynomial_degree=4, interpolation_is_log=True, polarized=False, time_like=False, n_integration_cores=0), debug=Debug(skip_singlet=False, skip_non_singlet=False), eko_version='0.0.0')\n" ] } ], @@ -143,7 +143,7 @@ ], "source": [ "with eko.EKO.read(\"./myeko.tar\") as evolution_operator:\n", - " with evolution_operator.operator((10000.,5)) as op:\n", + " with evolution_operator.operator((10000.0, 5)) as op:\n", " print(f\"operator: {op.operator.shape}\")\n", " print(f\"error: {op.error.shape}\")" ] diff --git a/doc/source/overview/tutorials/pdf.ipynb b/doc/source/overview/tutorials/pdf.ipynb index 4607ef891..bfc6ad3ad 100644 --- a/doc/source/overview/tutorials/pdf.ipynb +++ b/doc/source/overview/tutorials/pdf.ipynb @@ -13,9 +13,9 @@ "id": "3e0dcd0f", "metadata": {}, "source": [ - "## Method 1: Using apply_pdf\n", + "## Method 1: Using `apply_pdf`\n", "\n", - "In this first part we will compute the eko and subsequently apply the initial PDF \"manually\" calling a dedicated function. " + "In this first part, we compute the eko and subsequently apply the initial PDF \"manually\" calling a dedicated function. " ] }, { @@ -52,6 +52,7 @@ "import pathlib\n", "import eko\n", "from banana import toy\n", + "\n", "pdf = toy.mkPDF(\"\", 0)" ] }, @@ -71,6 +72,7 @@ "outputs": [], "source": [ "from ekobox.apply import apply_pdf\n", + "\n", "with eko.EKO.read(\"./myeko.tar\") as evolution_operator:\n", " evolved_pdfs = apply_pdf(evolution_operator, pdf)" ] @@ -92,7 +94,7 @@ { "data": { "text/plain": [ - "dict_keys([10000.0])" + "dict_keys([(10000.0, 5)])" ] }, "execution_count": 3, @@ -134,7 +136,7 @@ } ], "source": [ - "evolved_pdfs[10000.0][\"pdfs\"][21]" + "evolved_pdfs[(10000.0, 5)][\"pdfs\"][21]" ] }, { @@ -150,7 +152,7 @@ "id": "e925d2c9", "metadata": {}, "source": [ - "## Method 2: Using evolve_pdfs\n", + "## Method 2: Using `evolve_pdfs`\n", "\n", "In this second part we illustrate how to get (and install) directly a LHAPDF set evolved with eko. " ] @@ -171,6 +173,7 @@ "outputs": [], "source": [ "from banana import toy\n", + "\n", "pdf = toy.mkPDF(\"\", 0)" ] }, @@ -195,10 +198,10 @@ "th_card = example.theory()\n", "op_card = example.operator()\n", "# here we replace the grid with a very minimal one, to speed up the example\n", - "op_card.xgrid = eko.interpolation.XGrid([1e-3, 1e-2, 1e-1, 5e-1, 1.])\n", - "op_card.mugrid = [(10.,5), (100.,5)]\n", + "op_card.xgrid = eko.interpolation.XGrid([1e-3, 1e-2, 1e-1, 5e-1, 1.0])\n", + "op_card.mugrid = [(10.0, 5), (100.0, 5)]\n", "# set QCD LO evolution\n", - "th_card.orders = (1,0)" + "th_card.orders = (1, 0)" ] }, { @@ -236,16 +239,10 @@ ], "source": [ "from ekobox.evol_pdf import evolve_pdfs\n", + "\n", "path = pathlib.Path(\"./myeko2.tar\")\n", "path.unlink(missing_ok=True)\n", - "evolve_pdfs(\n", - " [pdf],\n", - " th_card,\n", - " op_card,\n", - " install=True,\n", - " name=\"Evolved_PDF\",\n", - " store_path=path\n", - ")" + "evolve_pdfs([pdf], th_card, op_card, install=True, name=\"Evolved_PDF\", store_path=path)" ] }, { @@ -273,6 +270,7 @@ ], "source": [ "import lhapdf\n", + "\n", "evolved_pdf = lhapdf.mkPDF(\"Evolved_PDF\", 0)" ] }, @@ -300,12 +298,12 @@ } ], "source": [ - "pid = 21 # gluon pid\n", - "Q2 = 89.10 # Q^2 in Gev^2\n", - "x = 0.01 # momentum fraction \n", + "pid = 21 # gluon pid\n", + "Q2 = 89.10 # Q^2 in Gev^2\n", + "x = 0.01 # momentum fraction\n", "\n", "# check that the particle is present\n", - "print(\"has gluon?\",evolved_pdf.hasFlavor(pid))\n", + "print(\"has gluon?\", evolved_pdf.hasFlavor(pid))\n", "# now do the lookup\n", "xg = evolved_pdf.xfxQ2(pid, x, Q2)\n", "print(f\"xg(x={x}, Q2={Q2}) = {xg}\")" @@ -352,28 +350,28 @@ "import lhapdf\n", "from ekobox.cards import example\n", "from eko.interpolation import make_grid\n", - "from eko.quantities.heavy_quarks import QuarkMassRef,HeavyQuarks\n", + "from eko.quantities.heavy_quarks import QuarkMassRef, HeavyQuarks\n", "\n", "# get the PDF object\n", "ct14llo = lhapdf.mkPDF(\"CT14llo\")\n", "\n", "# setup the operator card\n", "op_card = example.operator()\n", - "op_card.xgrid = eko.interpolation.XGrid(make_grid(30, 30)) # x grid\n", - "op_card.mugrid = [(float(q),5) for q in np.geomspace(5., 100, 5)] # Q2 grid\n", - "op_card.mu0 = 1.295000 # starting point for the evolution \n", + "op_card.xgrid = eko.interpolation.XGrid(make_grid(30, 30)) # x grid\n", + "op_card.mugrid = [(float(q), 5) for q in np.geomspace(5.0, 100, 5)] # Q2 grid\n", + "op_card.init = (1.295000, 3) # starting point for the evolution\n", "\n", "# setup the theory card - this can be mostly inferred from the PDF's .info file\n", - "\n", "th_card = example.theory()\n", - "th_card.orders = (1,0) # QCD LO\n", - "th_card.heavy.masses = HeavyQuarks([QuarkMassRef([1.3,nan]), QuarkMassRef([4.75,nan]), QuarkMassRef([172.,nan])]) # quark mass\n", - "th_card.couplings.alphas = 0.130000 # reference value of alpha_s\n", - "th_card.couplings.scale = 91.1876 # the reference scale at which alpha_s is provided\n", - "th_card.couplings.num_flavs_ref = 5 # the number of flavors active at the alpha_s reference scale\n", - "th_card.couplings.max_num_flavs = 5 # the maximum number of flavors active in the alpha_s evolution\n", - "th_card.couplings.num_flavs_init = 3 # the number of flavors active at the reference scale\n", - "th_card.num_flavs_max_pdf = 5 # the maximum number of flavors active in the pdf evolution." + "th_card.orders = (1, 0) # QCD LO\n", + "th_card.heavy.masses = HeavyQuarks(\n", + " [QuarkMassRef([1.3, nan]), QuarkMassRef([4.75, nan]), QuarkMassRef([172.0, nan])]\n", + ") # quark mass\n", + "th_card.couplings.alphas = 0.130000 # reference value of alpha_s\n", + "th_card.couplings.ref = (\n", + " 91.1876,\n", + " 5,\n", + ") # the reference scale together with the number of flavors at which alpha_s is provided" ] }, { @@ -407,15 +405,11 @@ ], "source": [ "from ekobox.evol_pdf import evolve_pdfs\n", + "\n", "path = pathlib.Path(\"./myeko_ct14llo.tar\")\n", "path.unlink(missing_ok=True)\n", "evolve_pdfs(\n", - " [ct14llo],\n", - " th_card,\n", - " op_card,\n", - " install=True,\n", - " name=\"my_ct14llo\",\n", - " store_path=path\n", + " [ct14llo], th_card, op_card, install=True, name=\"my_ct14llo\", store_path=path\n", ")" ] }, @@ -439,33 +433,33 @@ "output_type": "stream", "text": [ "LHAPDF 6.4.0 loading /home/felix/local/share/LHAPDF/my_ct14llo/my_ct14llo_0000.dat\n", + "my_ct14llo PDF set, member #0, version 1\n", " x Q2 ct14llo my_ct14llo relative_diff\n", - "0 0.000010 25.000000 7.635785e+01 7.630461e+01 0.000697\n", - "1 0.000173 25.000000 3.194273e+01 3.192092e+01 0.000683\n", - "2 0.003000 25.000000 1.081843e+01 1.081086e+01 0.000701\n", - "3 0.051962 25.000000 1.958956e+00 1.958632e+00 0.000166\n", - "4 0.900000 25.000000 1.922415e-05 1.955026e-05 -0.016963\n", - "5 0.000010 111.803399 1.333957e+02 1.332985e+02 0.000729\n", - "6 0.000173 111.803399 4.777286e+01 4.773664e+01 0.000758\n", - "7 0.003000 111.803399 1.341028e+01 1.339967e+01 0.000791\n", - "8 0.051962 111.803399 1.978216e+00 1.978130e+00 0.000044\n", - "9 0.900000 111.803399 6.644805e-06 6.753652e-06 -0.016381\n", - "10 0.000010 500.000000 1.967032e+02 1.965456e+02 0.000801\n", - "11 0.000173 500.000000 6.291393e+01 6.286095e+01 0.000842\n", - "12 0.003000 500.000000 1.542347e+01 1.540996e+01 0.000876\n", - "13 0.051962 500.000000 1.947465e+00 1.947391e+00 0.000038\n", - "14 0.900000 500.000000 2.929060e-06 2.977306e-06 -0.016471\n", - "15 0.000010 2236.067977 2.633266e+02 2.631109e+02 0.000819\n", - "16 0.000173 2236.067977 7.708540e+01 7.701938e+01 0.000856\n", - "17 0.003000 2236.067977 1.700410e+01 1.698928e+01 0.000872\n", - "18 0.051962 2236.067977 1.893923e+00 1.893971e+00 -0.000025\n", - "19 0.900000 2236.067977 1.544450e-06 1.570997e-06 -0.017189\n", - "20 0.000010 10000.000000 3.314097e+02 3.311351e+02 0.000829\n", - "21 0.000173 10000.000000 9.023010e+01 9.015279e+01 0.000857\n", - "22 0.003000 10000.000000 1.825934e+01 1.824402e+01 0.000839\n", - "23 0.051962 10000.000000 1.830992e+00 1.831183e+00 -0.000104\n", - "24 0.900000 10000.000000 9.288458e-07 9.447927e-07 -0.017169\n", - "my_ct14llo PDF set, member #0, version 1\n" + "0 0.000010 25.000000 7.635785e+01 7.630719e+01 0.000663\n", + "1 0.000173 25.000000 3.194273e+01 3.192239e+01 0.000637\n", + "2 0.003000 25.000000 1.081843e+01 1.081160e+01 0.000632\n", + "3 0.051962 25.000000 1.958956e+00 1.958820e+00 0.000069\n", + "4 0.900000 25.000000 1.922415e-05 1.955440e-05 -0.017179\n", + "5 0.000010 111.803399 1.333957e+02 1.333028e+02 0.000697\n", + "6 0.000173 111.803399 4.777286e+01 4.773855e+01 0.000718\n", + "7 0.003000 111.803399 1.341028e+01 1.340044e+01 0.000734\n", + "8 0.051962 111.803399 1.978216e+00 1.978292e+00 -0.000038\n", + "9 0.900000 111.803399 6.644805e-06 6.756354e-06 -0.016787\n", + "10 0.000010 500.000000 1.967032e+02 1.965517e+02 0.000770\n", + "11 0.000173 500.000000 6.291393e+01 6.286327e+01 0.000805\n", + "12 0.003000 500.000000 1.542347e+01 1.541073e+01 0.000826\n", + "13 0.051962 500.000000 1.947465e+00 1.947532e+00 -0.000034\n", + "14 0.900000 500.000000 2.929060e-06 2.979511e-06 -0.017224\n", + "15 0.000010 2236.067977 2.633266e+02 2.631189e+02 0.000789\n", + "16 0.000173 2236.067977 7.708540e+01 7.702204e+01 0.000822\n", + "17 0.003000 2236.067977 1.700410e+01 1.699004e+01 0.000827\n", + "18 0.051962 2236.067977 1.893923e+00 1.894094e+00 -0.000090\n", + "19 0.900000 2236.067977 1.544450e-06 1.572860e-06 -0.018395\n", + "20 0.000010 10000.000000 3.314097e+02 3.311450e+02 0.000799\n", + "21 0.000173 10000.000000 9.023010e+01 9.015576e+01 0.000824\n", + "22 0.003000 10000.000000 1.825934e+01 1.824477e+01 0.000798\n", + "23 0.051962 10000.000000 1.830992e+00 1.831291e+00 -0.000163\n", + "24 0.900000 10000.000000 9.288458e-07 9.463689e-07 -0.018866\n" ] } ], @@ -475,15 +469,15 @@ "# load evolved pdf\n", "my_ct14llo = lhapdf.mkPDF(\"my_ct14llo\", 0)\n", "\n", - "pid = 21 # gluon pid\n", + "pid = 21 # gluon pid\n", "\n", "# collect data\n", - "log = {\"x\": [], \"Q2\" : [], \"ct14llo\": [], \"my_ct14llo\": [], \"relative_diff\": []} \n", - "for q in np.geomspace(5., 100, 5):\n", - " q2 = q**2.\n", + "log = {\"x\": [], \"Q2\": [], \"ct14llo\": [], \"my_ct14llo\": [], \"relative_diff\": []}\n", + "for q in np.geomspace(5.0, 100, 5):\n", + " q2 = q**2.0\n", " for x in np.geomspace(1e-5, 0.9, 5):\n", " value = ct14llo.xfxQ2(pid, x, q2)\n", - " my_value = my_ct14llo.xfxQ2(pid, x, q2)\n", + " my_value = my_ct14llo.xfxQ2(pid, x, q2)\n", " log[\"x\"].append(x)\n", " log[\"Q2\"].append(q2)\n", " log[\"ct14llo\"].append(value)\n", diff --git a/extras/n3lo_bench/splitting_function_utils.py b/extras/n3lo_bench/splitting_function_utils.py index 6b32c86cc..66a516a5e 100644 --- a/extras/n3lo_bench/splitting_function_utils.py +++ b/extras/n3lo_bench/splitting_function_utils.py @@ -103,9 +103,7 @@ def compute_a_s(q2=None, xif2=1.0, nf=None, order=(4, 0)): ref = CouplingsInfo( alphas=0.1181, alphaem=0.007496, - scale=91.00, - max_num_flavs=6, - num_flavs_ref=5, + ref=(91.00, 5), ) sc = Couplings( couplings=ref, diff --git a/src/eko/couplings.py b/src/eko/couplings.py index a0826c09a..502695b32 100644 --- a/src/eko/couplings.py +++ b/src/eko/couplings.py @@ -441,7 +441,7 @@ def assert_positive(name, var): assert_positive("alpha_s_ref", couplings.alphas) assert_positive("alpha_em_ref", couplings.alphaem) - assert_positive("scale_ref", couplings.scale) + assert_positive("scale_ref", couplings.ref[0]) if order[0] not in [1, 2, 3, 4]: raise NotImplementedError( "QCD order has to be an integer between 1 (LO) and 4 (N3LO)" @@ -455,7 +455,7 @@ def assert_positive(name, var): raise ValueError(f"Unknown method {method.value}") self.method = method.value - nf_ref = couplings.num_flavs_ref + nf_ref = couplings.ref[1] scheme_name = hqm_scheme.name self.alphaem_running = couplings.em_running self.decoupled_running = False @@ -464,7 +464,7 @@ def assert_positive(name, var): self.a_ref = np.array(couplings.values) / 4.0 / np.pi # convert to a_s and a_em matching_scales = (np.array(masses) * np.array(thresholds_ratios)).tolist() self.thresholds_ratios = list(thresholds_ratios) - self.atlas = matchings.Atlas(matching_scales, (couplings.scale**2, nf_ref)) + self.atlas = matchings.Atlas(matching_scales, (couplings.ref[0] ** 2, nf_ref)) self.hqm_scheme = scheme_name logger.info( "Strong Coupling: a_s(µ_R^2=%f)%s=%f=%f/(4π)", diff --git a/src/eko/evolution_operator/grid.py b/src/eko/evolution_operator/grid.py index 6a9a3d9db..969750389 100644 --- a/src/eko/evolution_operator/grid.py +++ b/src/eko/evolution_operator/grid.py @@ -61,7 +61,6 @@ def __init__( masses: List[float], mass_scheme, thresholds_ratios: List[float], - intrinsic_flavors: list, xif: float, n3lo_ad_variation: tuple, matching_order: Order, @@ -75,7 +74,6 @@ def __init__( # check config: Dict[str, Any] = {} config["order"] = order - config["intrinsic_range"] = intrinsic_flavors config["xif2"] = xif**2 config["HQ"] = mass_scheme config["ModSV"] = configs.scvar_method @@ -175,16 +173,15 @@ def generate(self, q2: EPoint) -> OpDict: operator.compute() is_downward = is_downward_path(path) - intrinsic_range = [4, 5, 6] if is_downward else self.config["intrinsic_range"] qed = self.config["order"][1] > 0 final_op = physical.PhysicalOperator.ad_to_evol_map( - operator.op_members, operator.nf, operator.q2_to, intrinsic_range, qed + operator.op_members, operator.nf, operator.q2_to, qed ) # and multiply the lower ones from the right for op in reversed(list(thr_ops)): phys_op = physical.PhysicalOperator.ad_to_evol_map( - op.op_members, op.nf, op.q2_to, intrinsic_range, qed + op.op_members, op.nf, op.q2_to, qed ) # join with the basis rotation, since matching requires c+ (or likewise) @@ -193,7 +190,6 @@ def generate(self, q2: EPoint) -> OpDict: self._matching_operators[op.q2_to], nf_match, op.q2_to, - intrinsic_range=intrinsic_range, qed=qed, ) if is_downward: diff --git a/src/eko/evolution_operator/matching_condition.py b/src/eko/evolution_operator/matching_condition.py index 8df4775ff..e84ef4d2e 100644 --- a/src/eko/evolution_operator/matching_condition.py +++ b/src/eko/evolution_operator/matching_condition.py @@ -20,8 +20,7 @@ def split_ad_to_evol_map( op_members, nf, q2_thr, - intrinsic_range, - qed, + qed=False, ): """ Create the instance from the |OME|. @@ -35,8 +34,6 @@ def split_ad_to_evol_map( number of active flavors *below* the threshold q2_thr: float threshold value - intrinsic_range : list - list of intrinsic quark pids qed : bool activate qed """ @@ -78,27 +75,26 @@ def split_ad_to_evol_map( ) # intrinsic matching - if len(intrinsic_range) != 0: - op_id = member.OpMember.id_like(op_members[(200, 200)]) - for intr_fl in intrinsic_range: - ihq = br.quark_names[intr_fl - 1] # find name - if intr_fl > nf + 1: - # keep the higher quarks as they are - m[f"{ihq}+.{ihq}+"] = op_id.copy() - m[f"{ihq}-.{ihq}-"] = op_id.copy() - elif intr_fl == nf + 1: - # match the missing contribution from h+ and h- - m.update( - { - f"{ihq}+.{ihq}+": op_members[ - (br.matching_hplus_pid, br.matching_hplus_pid) - ], - f"S.{ihq}+": op_members[(100, br.matching_hplus_pid)], - f"g.{ihq}+": op_members[(21, br.matching_hplus_pid)], - f"{ihq}-.{ihq}-": op_members[ - (br.matching_hminus_pid, br.matching_hminus_pid) - ], - # f"V.{ihq}-": op_members[(200, br.matching_hminus_pid)], - } - ) + op_id = member.OpMember.id_like(op_members[(200, 200)]) + for intr_fl in [4, 5, 6]: + ihq = br.quark_names[intr_fl - 1] # find name + if intr_fl > nf + 1: + # keep the higher quarks as they are + m[f"{ihq}+.{ihq}+"] = op_id.copy() + m[f"{ihq}-.{ihq}-"] = op_id.copy() + elif intr_fl == nf + 1: + # match the missing contribution from h+ and h- + m.update( + { + f"{ihq}+.{ihq}+": op_members[ + (br.matching_hplus_pid, br.matching_hplus_pid) + ], + f"S.{ihq}+": op_members[(100, br.matching_hplus_pid)], + f"g.{ihq}+": op_members[(21, br.matching_hplus_pid)], + f"{ihq}-.{ihq}-": op_members[ + (br.matching_hminus_pid, br.matching_hminus_pid) + ], + # f"V.{ihq}-": op_members[(200, br.matching_hminus_pid)], + } + ) return cls.promote_names(m, q2_thr) diff --git a/src/eko/evolution_operator/operator_matrix_element.py b/src/eko/evolution_operator/operator_matrix_element.py index fc29aafc6..022f80730 100644 --- a/src/eko/evolution_operator/operator_matrix_element.py +++ b/src/eko/evolution_operator/operator_matrix_element.py @@ -247,10 +247,6 @@ def __init__(self, config, managers, nf, q2, is_backward, L, is_msbar): self.backward_method = matching_method( config["backward_inversion"] if is_backward else None ) - if is_backward: - self.is_intrinsic = True - else: - self.is_intrinsic = bool(len(config["intrinsic_range"]) != 0) self.L = L self.is_msbar = is_msbar # Note for the moment only QCD matching is implemented @@ -271,11 +267,10 @@ def labels(self): logger.warning("%s: skipping non-singlet sector", self.log_label) else: labels.append((200, 200)) - if self.is_intrinsic or self.backward_method is not MatchingMethods.FORWARD: - # intrinsic labels, which are not zero at NLO - labels.append((br.matching_hminus_pid, br.matching_hminus_pid)) - # These contributions are always 0 for the moment - # labels.extend([(200, br.matching_hminus_pid), (br.matching_hminus_pid, 200)]) + # intrinsic labels, which are not zero at NLO + labels.append((br.matching_hminus_pid, br.matching_hminus_pid)) + # These contributions are always 0 for the moment + # labels.extend([(200, br.matching_hminus_pid), (br.matching_hminus_pid, 200)]) # same for singlet if self.config["debug_skip_singlet"]: logger.warning("%s: skipping singlet sector", self.log_label) @@ -287,14 +282,13 @@ def labels(self): (br.matching_hplus_pid, 100), ] ) - if self.is_intrinsic or self.backward_method is not MatchingMethods.FORWARD: - labels.extend( - [ - (21, br.matching_hplus_pid), - (100, br.matching_hplus_pid), - (br.matching_hplus_pid, br.matching_hplus_pid), - ] - ) + labels.extend( + [ + (21, br.matching_hplus_pid), + (100, br.matching_hplus_pid), + (br.matching_hplus_pid, br.matching_hplus_pid), + ] + ) return labels def quad_ker(self, label, logx, areas): diff --git a/src/eko/evolution_operator/physical.py b/src/eko/evolution_operator/physical.py index ca01c4ea4..f8a3545d8 100644 --- a/src/eko/evolution_operator/physical.py +++ b/src/eko/evolution_operator/physical.py @@ -23,7 +23,7 @@ class PhysicalOperator(member.OperatorBase): """ @classmethod - def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range, qed): + def ad_to_evol_map(cls, op_members, nf, q2_final, qed=False): """ Obtain map between the 3-dimensional anomalous dimension basis and the 4-dimensional evolution basis. @@ -35,8 +35,6 @@ def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range, qed): operator members in anomalous dimension basis nf : int number of active light flavors - intrinsic_range : sequence - intrinsic heavy flavors qed : bool activate qed @@ -94,15 +92,14 @@ def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range, qed): m["Tu8.Tu8"] = op_members[(br.non_singlet_pids_map["ns+u"], 0)] m["Vu8.Vu8"] = op_members[(br.non_singlet_pids_map["ns-u"], 0)] # deal with intrinsic heavy quark PDFs - if intrinsic_range is not None: - hqfl = "cbt" - op_id = member.OpMember.id_like(op_members[(21, 21)]) - for intr_fl in intrinsic_range: - if intr_fl <= nf: # light quarks are not intrinsic - continue - hq = hqfl[intr_fl - 4] # find name - # intrinsic means no evolution, i.e. they are evolving with the identity - m[f"{hq}+.{hq}+"] = op_id.copy() - m[f"{hq}-.{hq}-"] = op_id.copy() + hqfl = "cbt" + op_id = member.OpMember.id_like(op_members[(21, 21)]) + for intr_fl in [4, 5, 6]: + if intr_fl <= nf: # light quarks are not intrinsic + continue + hq = hqfl[intr_fl - 4] # find name + # intrinsic means no evolution, i.e. they are evolving with the identity + m[f"{hq}+.{hq}+"] = op_id.copy() + m[f"{hq}-.{hq}-"] = op_id.copy() # map key to MemberName return cls.promote_names(m, q2_final) diff --git a/src/eko/io/legacy.py b/src/eko/io/legacy.py index 0458b6a78..5b15f154d 100644 --- a/src/eko/io/legacy.py +++ b/src/eko/io/legacy.py @@ -100,9 +100,6 @@ class PseudoTheory(DictLike): def from_old(cls, old: RawCard): """Load from old metadata.""" heavy = HeavyInfo( - num_flavs_init=4, - num_flavs_max_pdf=5, - intrinsic_flavors=[], masses=HeavyQuarkMasses( [ ReferenceRunning([_MC, np.inf]), @@ -125,7 +122,7 @@ class PseudoOperator(DictLike): """ - mu20: float + init: EPoint evolgrid: List[EPoint] xgrid: XGrid configs: dict @@ -133,7 +130,7 @@ class PseudoOperator(DictLike): @classmethod def from_old(cls, old: RawCard): """Load from old metadata.""" - mu20 = float(old["q2_ref"]) + mu0 = float(np.sqrt(float(old["q2_ref"]))) mu2list = old.get("Q2grid") if mu2list is None: mu2list = old["mu2grid"] @@ -149,7 +146,7 @@ def from_old(cls, old: RawCard): interpolation_is_log=old.get("interpolation_is_log"), ) - return cls(mu20=mu20, evolgrid=evolgrid, xgrid=xgrid, configs=configs) + return cls(init=(mu0, 4), evolgrid=evolgrid, xgrid=xgrid, configs=configs) ARRAY_SUFFIX = ".npy.lz4" diff --git a/src/eko/io/runcards.py b/src/eko/io/runcards.py index e48437dc3..e9fac3893 100644 --- a/src/eko/io/runcards.py +++ b/src/eko/io/runcards.py @@ -105,7 +105,7 @@ class Configs(DictLike): class OperatorCard(DictLike): """Operator Card info.""" - mu0: float + init: EPoint """Initial scale.""" mugrid: List[EPoint] xgrid: interpolation.XGrid @@ -126,7 +126,7 @@ class OperatorCard(DictLike): @property def mu20(self): """Squared value of initial scale.""" - return self.mu0**2 + return self.init[0] ** 2 @property def mu2grid(self) -> npt.NDArray: @@ -194,25 +194,12 @@ def new_theory(self): alphas=old["alphas"], alphaem=alphaem, em_running=em_running, - scale=old["Qref"], - num_flavs_ref=old["nfref"], - max_num_flavs=old["MaxNfAs"], + ref=(old["Qref"], old["nfref"]), ) new["heavy"] = { - "num_flavs_init": ( - nf_default(old["Q0"] ** 2.0, default_atlas(ms, ks)) - if old["nf0"] is None - else old["nf0"] - ), - "num_flavs_max_pdf": old["MaxNfPdf"], "matching_ratios": self.heavies("k%sThr", old), "masses_scheme": old["HQ"], } - intrinsic = [] - for idx, q in enumerate(hq.FLAVORS): - if old.get(f"i{q}".upper()) == 1: - intrinsic.append(idx + 4) - new["heavy"]["intrinsic_flavors"] = intrinsic if old["HQ"] == "POLE": new["heavy"]["masses"] = [[m, nan] for m in ms] elif old["HQ"] == "MSBAR": @@ -236,17 +223,22 @@ def new_operator(self): old_th = self.theory new = {} - new["mu0"] = old_th["Q0"] + ms = self.heavies("m%s", old_th) + ks = self.heavies("k%sThr", old_th) + new["init"] = ( + old_th["Q0"], + ( + nf_default(old_th["Q0"] ** 2.0, default_atlas(ms, ks)) + if old_th["nf0"] is None + else old_th["nf0"] + ), + ) if "mugrid" in old: mugrid = old["mugrid"] else: mu2grid = old["Q2grid"] if "Q2grid" in old else old["mu2grid"] mugrid = np.sqrt(mu2grid).tolist() - new["mugrid"] = flavored_mugrid( - mugrid, - list(self.heavies("m%s", old_th)), - list(self.heavies("k%sThr", old_th)), - ) + new["mugrid"] = flavored_mugrid(mugrid, list(ms), list(ks)) new["configs"] = {} evmod = old_th["ModEv"] @@ -278,27 +270,6 @@ def new_operator(self): return OperatorCard.from_dict(new) -def update(theory: Union[RawCard, TheoryCard], operator: Union[RawCard, OperatorCard]): - """Update legacy runcards. - - This function is mainly defined for compatibility with the old interface. - Prefer direct usage of :class:`Legacy` in new code. - - Consecutive applications of this function yield identical results:: - - cards = update(theory, operator) - assert update(*cards) == cards - """ - if isinstance(theory, TheoryCard) or isinstance(operator, OperatorCard): - # if one is not a dict, both have to be new cards - assert isinstance(theory, TheoryCard) - assert isinstance(operator, OperatorCard) - return theory, operator - - cards = Legacy(theory, operator) - return cards.new_theory, cards.new_operator - - def default_atlas(masses: list, matching_ratios: list): r"""Create default landscape. diff --git a/src/eko/io/struct.py b/src/eko/io/struct.py index 1520e1003..fe3bf7c6e 100644 --- a/src/eko/io/struct.py +++ b/src/eko/io/struct.py @@ -547,7 +547,7 @@ def build(self) -> EKO: self.access.open = True metadata = Metadata( _path=self.path, - origin=(self.operator.mu20, self.theory.heavy.num_flavs_init), + origin=(self.operator.init[0] ** 2, self.operator.init[1]), xgrid=self.operator.xgrid, ) InternalPaths(self.path).bootstrap( diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py index de605764b..63e1b72d1 100644 --- a/src/eko/msbar_masses.py +++ b/src/eko/msbar_masses.py @@ -378,8 +378,8 @@ def compute( """ # TODO: sketch in the docs how the MSbar computation works with a figure. - mu2_ref = couplings.scale**2 - nf_ref: FlavorsNumber = couplings.num_flavs_ref + mu2_ref = couplings.ref[0] ** 2 + nf_ref: FlavorsNumber = couplings.ref[1] masses = np.concatenate((np.zeros(nf_ref - 3), np.full(6 - nf_ref, np.inf))) def sc(thr_masses): diff --git a/src/eko/quantities/couplings.py b/src/eko/quantities/couplings.py index b223ce2a7..cbeb40a62 100644 --- a/src/eko/quantities/couplings.py +++ b/src/eko/quantities/couplings.py @@ -4,7 +4,8 @@ import enum from ..io.dictlike import DictLike -from ..io.types import FlavorsNumber, LinearScale, ReferenceRunning, Scalar +from ..io.types import EvolutionPoint as EPoint +from ..io.types import LinearScale, ReferenceRunning, Scalar Coupling = Scalar CouplingRef = ReferenceRunning[Coupling] @@ -20,15 +21,7 @@ class CouplingsInfo(DictLike): alphas: Coupling alphaem: Coupling - scale: LinearScale - max_num_flavs: FlavorsNumber - num_flavs_ref: FlavorsNumber - r"""Number of active flavors at strong coupling reference scale. - - I.e. :math:`n_{f,\text{ref}}(\mu^2_{\text{ref}})`, formerly called - ``nfref``. - - """ + ref: EPoint em_running: bool = False @property diff --git a/src/eko/quantities/heavy_quarks.py b/src/eko/quantities/heavy_quarks.py index 8db084ed4..7328f0f12 100644 --- a/src/eko/quantities/heavy_quarks.py +++ b/src/eko/quantities/heavy_quarks.py @@ -7,13 +7,7 @@ import numpy as np from ..io.dictlike import DictLike -from ..io.types import ( - FlavorsNumber, - IntrinsicFlavors, - LinearScale, - ReferenceRunning, - SquaredScale, -) +from ..io.types import FlavorsNumber, LinearScale, ReferenceRunning, SquaredScale FLAVORS = "cbt" @@ -91,16 +85,6 @@ class HeavyInfo(DictLike): """ - num_flavs_init: FlavorsNumber - r"""Number of active flavors at fitting scale. - - I.e. :math:`n_{f,\text{ref}}(\mu^2_0)`, formerly called ``nf0``. - - """ - num_flavs_max_pdf: FlavorsNumber - """Maximum number of quark PDFs.""" - intrinsic_flavors: IntrinsicFlavors - """List of intrinsic quark PDFs.""" masses: HeavyQuarkMasses """List of heavy quark masses.""" masses_scheme: QuarkMassScheme diff --git a/src/eko/runner/commons.py b/src/eko/runner/commons.py index 543fb927c..af547b9a6 100644 --- a/src/eko/runner/commons.py +++ b/src/eko/runner/commons.py @@ -33,7 +33,7 @@ def atlas(theory: TheoryCard, operator: OperatorCard) -> Atlas: # TODO: cache result masses = runcards.masses(theory, operator.configs.evolution_method) matching_scales = np.power(theory.heavy.matching_ratios, 2.0) * np.array(masses) - return Atlas(matching_scales.tolist(), (operator.mu20, theory.heavy.num_flavs_init)) + return Atlas(matching_scales.tolist(), (operator.mu20, operator.init[1])) def couplings(theory: TheoryCard, operator: OperatorCard) -> Couplings: diff --git a/src/eko/runner/legacy.py b/src/eko/runner/legacy.py index 9ce925a9f..b6e9bf1c4 100644 --- a/src/eko/runner/legacy.py +++ b/src/eko/runner/legacy.py @@ -4,6 +4,8 @@ import os from typing import Union +from ekomark.data import update_runcards + from ..evolution_operator.grid import OperatorGrid from ..io import EKO, Operator, runcards from ..io.types import RawCard @@ -44,8 +46,7 @@ def __init__( path where to store the computed operator """ - new_theory, new_operator = runcards.update(theory_card, operators_card) - new_theory.heavy.intrinsic_flavors = [4, 5, 6] + new_theory, new_operator = update_runcards(theory_card, operators_card) # Store inputs self.path = path @@ -70,7 +71,6 @@ def __init__( masses=masses, mass_scheme=new_theory.heavy.masses_scheme.value, thresholds_ratios=new_theory.heavy.squared_ratios, - intrinsic_flavors=new_theory.heavy.intrinsic_flavors, xif=new_theory.xif, configs=new_operator.configs, debug=new_operator.debug, diff --git a/src/eko/runner/managed.py b/src/eko/runner/managed.py index 486c8035c..df5223700 100644 --- a/src/eko/runner/managed.py +++ b/src/eko/runner/managed.py @@ -21,8 +21,6 @@ def solve(theory: TheoryCard, operator: OperatorCard, path: Path): """Solve DGLAP equations in terms of evolution kernel operators (EKO).""" - theory.heavy.intrinsic_flavors = [4, 5, 6] - with EKO.create(path) as builder: eko = builder.load_cards(theory, operator).build() # pylint: disable=E1101 diff --git a/src/eko/runner/parts.py b/src/eko/runner/parts.py index 798353f71..18d22f377 100644 --- a/src/eko/runner/parts.py +++ b/src/eko/runner/parts.py @@ -39,31 +39,6 @@ def managers(eko: EKO) -> Managers: ) -def blowup_info(eko: EKO) -> dict: - """Prepare common information to blow up to flavor basis. - - Note - ---- - ``intrinsic_range`` is a fully deprecated feature, here and anywhere else, - since a full range is already always used for backward evolution, and it is - not harmful to use it also for forward. - - Indeed, the only feature of non-intrinsic evolution is to absorb a - non-trivial boundary condition when an intrinsic PDF is defined. - But to achieve this, is sufficient to not specify any intrinsic boundary - condition at all, while if something is there, it is intuitive enough that - it will be consistently evolved. - - Moreover, since two different behavior are applied for the forward and - backward evolution, the intrinsic range is a "non-local" function, since it - does not depend only on the evolution segment, but also on the previous - evolution history (to determine if evolution is backward in flavor, - irrespectively of happening for an increasing or decreasing interval in - scale at fixed flavor). - """ - return dict(intrinsic_range=[4, 5, 6], qed=eko.theory_card.order[1] > 0) - - def evolve_configs(eko: EKO) -> dict: """Create configs for :class:`Operator`. @@ -99,10 +74,10 @@ def evolve(eko: EKO, recipe: Evolution) -> Operator: ) op.compute() - binfo = blowup_info(eko) + qed = eko.theory_card.order[1] > 0 res, err = physical.PhysicalOperator.ad_to_evol_map( - op.op_members, op.nf, op.q2_to, **binfo - ).to_flavor_basis_tensor(qed=binfo["qed"]) + op.op_members, op.nf, op.q2_to, qed + ).to_flavor_basis_tensor(qed) return Operator(res, err) @@ -119,7 +94,6 @@ def matching_configs(eko: EKO) -> dict: return dict( **evolve_configs(eko), backward_inversion=ocard.configs.inversion_method, - intrinsic_range=tcard.heavy.intrinsic_flavors, ) @@ -146,10 +120,9 @@ def match(eko: EKO, recipe: Matching) -> Operator: eko.theory_card.heavy.masses_scheme is QuarkMassScheme.MSBAR, ) op.compute() - - binfo = blowup_info(eko) + qed = eko.theory_card.order[1] > 0 res, err = matching_condition.MatchingCondition.split_ad_to_evol_map( - op.op_members, op.nf, recipe.scale, **binfo - ).to_flavor_basis_tensor(qed=binfo["qed"]) + op.op_members, op.nf, recipe.scale, qed + ).to_flavor_basis_tensor(qed) return Operator(res, err) diff --git a/src/ekobox/cards.py b/src/ekobox/cards.py index 9476bb77a..86233c0df 100644 --- a/src/ekobox/cards.py +++ b/src/ekobox/cards.py @@ -14,14 +14,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=[4], masses=[ReferenceRunning([mq, nan]) for mq in (2.0, 4.5, 173.07)], masses_scheme="POLE", matching_ratios=[1.0, 1.0, 1.0], @@ -33,7 +28,7 @@ ) _operator = dict( - mu0=1.65, + init=(1.65, 4), mugrid=[(100.0, 5)], xgrid=np.geomspace(1e-7, 1.0, 50).tolist(), configs=dict( diff --git a/src/ekobox/cli/runcards.py b/src/ekobox/cli/runcards.py index 227e012f4..99db353bd 100644 --- a/src/ekobox/cli/runcards.py +++ b/src/ekobox/cli/runcards.py @@ -36,7 +36,7 @@ def sub_example(destination: pathlib.Path): theory.order = (1, 0) cards.dump(theory.raw, path=destination / "theory.yaml") operator = cards.example.operator() - operator.mu0 = 1.65 + operator.init = (1.65, 4) operator.mugrid = [(np.sqrt(1e5), 5)] cards.dump(operator.raw, path=destination / "operator.yaml") _logger.info(f"Runcards generated to '{destination}'") diff --git a/src/ekobox/info_file.py b/src/ekobox/info_file.py index 56439599d..16dde0197 100644 --- a/src/ekobox/info_file.py +++ b/src/ekobox/info_file.py @@ -38,7 +38,9 @@ def build( info file in lhapdf format """ template_info = copy.deepcopy(load.template_info) - template_info["SetDesc"] = "Evolved PDF from " + str(operators_card.mu0) + " GeV" + template_info["SetDesc"] = ( + "Evolved PDF from " + str(operators_card.init[0]) + " GeV" + ) template_info["Authors"] = "" template_info["FlavorScheme"] = "variable" template_info.update(info_update) @@ -51,7 +53,7 @@ def build( template_info["OrderQCD"] = theory_card.order[0] - 1 template_info["QMin"] = round(math.sqrt(operators_card.mu2grid[0]), 4) template_info["QMax"] = round(math.sqrt(operators_card.mu2grid[-1]), 4) - template_info["MZ"] = theory_card.couplings.scale + template_info["MZ"] = theory_card.couplings.ref[0] template_info["MUp"] = 0.0 template_info["MDown"] = 0.0 template_info["MStrange"] = 0.0 diff --git a/src/ekobox/utils.py b/src/ekobox/utils.py index 091cde15e..10ee536a3 100644 --- a/src/ekobox/utils.py +++ b/src/ekobox/utils.py @@ -39,7 +39,7 @@ def ekos_product( # another kind of output which includes the theory and operator runcards) ep_match = eko_ini.approx( - (eko_fin.operator_card.mu0**2, eko_fin.theory_card.heavy.num_flavs_init), + (eko_fin.operator_card.init[0] ** 2, eko_fin.operator_card.init[1]), rtol=rtol, atol=atol, ) diff --git a/src/ekomark/data/__init__.py b/src/ekomark/data/__init__.py index 797409204..89fb6377b 100644 --- a/src/ekomark/data/__init__.py +++ b/src/ekomark/data/__init__.py @@ -1 +1,30 @@ """EKO database configuration.""" + +from typing import Union + +from eko.io.runcards import Legacy, OperatorCard, TheoryCard +from eko.io.types import RawCard + + +def update_runcards( + theory: Union[RawCard, TheoryCard], operator: Union[RawCard, OperatorCard] +): + """Update legacy runcards. + + This function is mainly defined for compatibility with the old interface. + Prefer direct usage of :class:`Legacy` in new code. + + Consecutive applications of this function yield identical results:: + + cards = update(theory, operator) + assert update(*cards) == cards + + """ + if isinstance(theory, TheoryCard) or isinstance(operator, OperatorCard): + # if one is not a dict, both have to be new cards + assert isinstance(theory, TheoryCard) + assert isinstance(operator, OperatorCard) + return theory, operator + + cards = Legacy(theory, operator) + return cards.new_theory, cards.new_operator diff --git a/tests/eko/evolution_operator/test_grid.py b/tests/eko/evolution_operator/test_grid.py index 804388db6..01754578a 100644 --- a/tests/eko/evolution_operator/test_grid.py +++ b/tests/eko/evolution_operator/test_grid.py @@ -72,7 +72,7 @@ def test_mod_expanded(theory_card, theory_ffns, operator_card, tmp_path: pathlib else: theory = theory_card theory.order = (1, 0) - theory.heavy.num_flavs_init = nf0 + operator_card.init = (operator_card.init[0], nf0) path.unlink(missing_ok=True) opgrid = legacy.Runner(theory, operator_card, path=path).op_grid opg = opgrid.compute() diff --git a/tests/eko/evolution_operator/test_matching_condition.py b/tests/eko/evolution_operator/test_matching_condition.py index c854fb90c..bb6aeebbc 100644 --- a/tests/eko/evolution_operator/test_matching_condition.py +++ b/tests/eko/evolution_operator/test_matching_condition.py @@ -22,12 +22,6 @@ def mkOME(self): (br.matching_hplus_pid, 21), (200, 200), (br.matching_hminus_pid, 200), - ]: - ome.update({key: mkOM(self.shape)}) - return ome - - def update_intrinsic_OME(self, ome): - for key in [ (br.matching_hplus_pid, br.matching_hplus_pid), (br.matching_hminus_pid, br.matching_hminus_pid), (200, br.matching_hminus_pid), @@ -35,10 +29,11 @@ def update_intrinsic_OME(self, ome): (21, br.matching_hplus_pid), ]: ome.update({key: mkOM(self.shape)}) + return ome def test_split_ad_to_evol_map(self): ome = self.mkOME() - a = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [], False) + a = MatchingCondition.split_ad_to_evol_map(ome, 3, 1) triv_keys = [ "V.V", "T3.T3", @@ -55,6 +50,14 @@ def test_split_ad_to_evol_map(self): "c+.S", "c+.g", # "c-.V", + "S.c+", + "g.c+", + "c+.c+", + "c-.c-", + "b+.b+", + "b-.b-", + "t+.t+", + "t-.t-", ] assert sorted(str(k) for k in a.op_members.keys()) == sorted( [*triv_keys, *keys3] @@ -63,36 +66,8 @@ def test_split_ad_to_evol_map(self): a.op_members[member.MemberName("V.V")].value, ome[(200, 200)].value, ) - # # if alpha is zero, nothing non-trivial should happen - b = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [], False) - assert sorted(str(k) for k in b.op_members.keys()) == sorted( - [*triv_keys, *keys3] - ) - # assert_almost_equal( - # b.op_members[member.MemberName("V.V")].value, - # np.eye(self.shape[0]), - # ) - # nf=3 + IC - self.update_intrinsic_OME(ome) - c = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [4], False) - assert sorted(str(k) for k in c.op_members.keys()) == sorted( - [*triv_keys, *keys3, "S.c+", "g.c+", "c+.c+", "c-.c-"] - ) - assert_almost_equal( - c.op_members[member.MemberName("V.V")].value, - b.op_members[member.MemberName("V.V")].value, - ) - # nf=3 + IB - d = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [5], False) - assert sorted(str(k) for k in d.op_members.keys()) == sorted( - [*triv_keys, *keys3, "b+.b+", "b-.b-"] - ) - assert_almost_equal( - d.op_members[member.MemberName("b+.b+")].value, - np.eye(self.shape[0]), - ) # nf=4 + IB - d = MatchingCondition.split_ad_to_evol_map(ome, 4, 1, [5], False) + d = MatchingCondition.split_ad_to_evol_map(ome, 4, 1) assert sorted(str(k) for k in d.op_members.keys()) == sorted( [ *triv_keys, @@ -106,6 +81,8 @@ def test_split_ad_to_evol_map(self): "b+.b+", # "b-.V", "b-.b-", + "t+.t+", + "t-.t-", ] ) assert_almost_equal( @@ -127,7 +104,7 @@ def test_split_ad_to_evol_map(self): def test_split_ad_to_evol_map_qed(self): ome = self.mkOME() - a = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [], qed=True) + a = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, qed=True) triv_keys = [ "ph.ph", "S.S", @@ -144,7 +121,15 @@ def test_split_ad_to_evol_map_qed(self): keys3 = [ "c+.S", "c+.g", + "S.c+", + "g.c+", + "c+.c+", + "c-.c-", # "c-.V", + "b+.b+", + "b-.b-", + "t+.t+", + "t-.t-", ] assert sorted(str(k) for k in a.op_members.keys()) == sorted( [*triv_keys, *keys3] @@ -153,36 +138,8 @@ def test_split_ad_to_evol_map_qed(self): a.op_members[member.MemberName("V.V")].value, ome[(200, 200)].value, ) - # # if alpha is zero, nothing non-trivial should happen - b = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [], qed=True) - assert sorted(str(k) for k in b.op_members.keys()) == sorted( - [*triv_keys, *keys3] - ) - # assert_almost_equal( - # b.op_members[member.MemberName("V.V")].value, - # np.eye(self.shape[0]), - # ) - # nf=3 + IC - self.update_intrinsic_OME(ome) - c = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [4], qed=True) - assert sorted(str(k) for k in c.op_members.keys()) == sorted( - [*triv_keys, *keys3, "S.c+", "g.c+", "c+.c+", "c-.c-"] - ) - assert_almost_equal( - c.op_members[member.MemberName("V.V")].value, - b.op_members[member.MemberName("V.V")].value, - ) - # nf=3 + IB - d = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [5], qed=True) - assert sorted(str(k) for k in d.op_members.keys()) == sorted( - [*triv_keys, *keys3, "b+.b+", "b-.b-"] - ) - assert_almost_equal( - d.op_members[member.MemberName("b+.b+")].value, - np.eye(self.shape[0]), - ) # nf=4 + IB - d = MatchingCondition.split_ad_to_evol_map(ome, 4, 1, [5], qed=True) + d = MatchingCondition.split_ad_to_evol_map(ome, 4, 1, qed=True) assert sorted(str(k) for k in d.op_members.keys()) == sorted( [ *triv_keys, @@ -196,6 +153,8 @@ def test_split_ad_to_evol_map_qed(self): "b+.b+", # "b-.V", "b-.b-", + "t+.t+", + "t-.t-", ] ) assert_almost_equal( diff --git a/tests/eko/evolution_operator/test_ome.py b/tests/eko/evolution_operator/test_ome.py index 553e5823c..e807a587b 100644 --- a/tests/eko/evolution_operator/test_ome.py +++ b/tests/eko/evolution_operator/test_ome.py @@ -339,6 +339,7 @@ def test_labels(self, theory_ffns, operator_card, tmp_path: pathlib.Path): path = tmp_path / "eko.tar" for skip_singlet in [True, False]: for skip_ns in [True, False]: + operator_card.configs.inversion_method = InversionMethod.EXACT operator_card.debug.skip_singlet = skip_singlet operator_card.debug.skip_non_singlet = skip_ns path.unlink(missing_ok=True) diff --git a/tests/eko/evolution_operator/test_physical.py b/tests/eko/evolution_operator/test_physical.py index e9882de5b..0ba41d91f 100644 --- a/tests/eko/evolution_operator/test_physical.py +++ b/tests/eko/evolution_operator/test_physical.py @@ -221,27 +221,21 @@ def mk_op_members(shape=(2, 2), qed=False): return om -def get_ad_to_evol_map(nf, intrinsic_range=None, qed=False): +def get_ad_to_evol_map(nf, qed=False): oms = mk_op_members(qed=qed) - m = PhysicalOperator.ad_to_evol_map(oms, nf, 1, intrinsic_range, qed) + m = PhysicalOperator.ad_to_evol_map(oms, nf, 1, qed) return sorted(map(str, m.op_members.keys())) def test_ad_to_evol_map(): triv_ops = ("S.S", "S.g", "g.S", "g.g", "V.V", "V3.V3", "T3.T3", "V8.V8", "T8.T8") # nf=3 - assert sorted(triv_ops) == get_ad_to_evol_map(3) - # nf=3 + IC - assert sorted([*triv_ops, "c+.c+", "c-.c-"]) == get_ad_to_evol_map(3, [4]) - # nf=3 + IC + IB assert sorted( - [*triv_ops, "c+.c+", "c-.c-", "b+.b+", "b-.b-"] - ) == get_ad_to_evol_map(3, [4, 5]) - # nf=4 + IC(non-existant) + IB - ks = sorted([*triv_ops, "V15.V15", "T15.T15", "b+.b+", "b-.b-"]) - assert ks == get_ad_to_evol_map(4, [4, 5]) - # nf=4 + IB - assert ks == get_ad_to_evol_map(4, [5]) + [*triv_ops, "c+.c+", "c-.c-", "b+.b+", "b-.b-", "t+.t+", "t-.t-"] + ) == get_ad_to_evol_map(3) + # nf=4 + ks = sorted([*triv_ops, "V15.V15", "T15.T15", "b+.b+", "b-.b-", "t+.t+", "t-.t-"]) + assert ks == get_ad_to_evol_map(4) # nf=6 assert sorted( [*triv_ops, "T15.T15", "V15.V15", "T24.T24", "V24.V24", "T35.T35", "V35.V35"] @@ -274,18 +268,12 @@ def test_ad_to_evol_map_qed(): "Td3.Td3", ) # nf=3 - assert sorted(triv_ops) == get_ad_to_evol_map(3, qed=True) - # nf=3 + IC - assert sorted([*triv_ops, "c+.c+", "c-.c-"]) == get_ad_to_evol_map(3, [4], qed=True) - # nf=3 + IC + IB assert sorted( - [*triv_ops, "c+.c+", "c-.c-", "b+.b+", "b-.b-"] - ) == get_ad_to_evol_map(3, [4, 5], qed=True) - # nf=4 + IC(non-existant) + IB - ks = sorted([*triv_ops, "Vu3.Vu3", "Tu3.Tu3", "b+.b+", "b-.b-"]) - assert ks == get_ad_to_evol_map(4, [4, 5], qed=True) - # nf=4 + IB - assert ks == get_ad_to_evol_map(4, [5], qed=True) + [*triv_ops, "c+.c+", "c-.c-", "b+.b+", "b-.b-", "t+.t+", "t-.t-"] + ) == get_ad_to_evol_map(3, True) + # nf=4 + ks = sorted([*triv_ops, "Vu3.Vu3", "Tu3.Tu3", "b+.b+", "b-.b-", "t+.t+", "t-.t-"]) + assert ks == get_ad_to_evol_map(4, True) # nf=6 assert sorted( [*triv_ops, "Tu3.Tu3", "Vu3.Vu3", "Td8.Td8", "Vd8.Vd8", "Tu8.Tu8", "Vu8.Vu8"] diff --git a/tests/eko/io/test_runcards.py b/tests/eko/io/test_runcards.py index 4b4996336..4fe9a8e4b 100644 --- a/tests/eko/io/test_runcards.py +++ b/tests/eko/io/test_runcards.py @@ -1,13 +1,10 @@ -import copy from io import StringIO import numpy as np import pytest import yaml -from banana.data.theories import default_card as theory_card from eko.io import runcards as rc -from ekomark.data.operators import default_card as operator_card def test_flavored_mu2grid(): @@ -30,48 +27,6 @@ def test_flavored_mu2grid(): assert t == tuple(l) -def test_runcards_opcard(): - # check conversion - tc = copy.deepcopy(theory_card) - oc = copy.deepcopy(operator_card) - tc["Q0"] = 2.0 - # mugrid - oc["mugrid"] = [2.0, 10.0] - _nt, no = rc.update(tc, oc) - assert isinstance(no, rc.OperatorCard) - assert len(no.evolgrid) == len(oc["mugrid"]) - assert len(no.mu2grid) == len(no.evolgrid) - assert no.evolgrid[0][-1] == 4 - assert no.evolgrid[1][-1] == 5 - np.testing.assert_allclose(no.mu0, tc["Q0"]) - np.testing.assert_allclose(no.mu20, tc["Q0"] ** 2.0) - assert len(no.pids) == 14 - check_dumpable(no) - del oc["mugrid"] - # or mu2grid - oc["mu2grid"] = [9.0, 30.0, 32.0] - _nt, no = rc.update(tc, oc) - assert isinstance(no, rc.OperatorCard) - assert len(no.evolgrid) == len(oc["mu2grid"]) - assert len(no.mu2grid) == len(no.evolgrid) - assert no.evolgrid[0][-1] == 4 - assert no.evolgrid[1][-1] == 5 - assert no.evolgrid[2][-1] == 5 - check_dumpable(no) - del oc["mu2grid"] - # or Q2grid - oc["Q2grid"] = [15.0, 130.0, 140.0, 1e5] - _nt, no = rc.update(tc, oc) - assert isinstance(no, rc.OperatorCard) - assert len(no.evolgrid) == len(oc["Q2grid"]) - assert len(no.mu2grid) == len(no.evolgrid) - assert no.evolgrid[0][-1] == 4 - assert no.evolgrid[1][-1] == 5 - assert no.evolgrid[2][-1] == 5 - assert no.evolgrid[3][-1] == 6 - check_dumpable(no) - - def check_dumpable(no): """Check we can write and read to yaml.""" so = StringIO() @@ -80,50 +35,6 @@ def check_dumpable(no): noo = yaml.safe_load(so) -def test_runcards_ekomark(): - # check conversion - tc = copy.deepcopy(theory_card) - oc = copy.deepcopy(operator_card) - nt, no = rc.update(tc, oc) - assert isinstance(nt, rc.TheoryCard) - assert isinstance(no, rc.OperatorCard) - # check is idempotent - nnt, nno = rc.update(nt, no) - assert nnt == nt - assert nno == no - - -def test_runcards_quarkmass(): - tc = copy.deepcopy(theory_card) - tc["nfref"] = 5 - tc["IC"] = 1 - oc = copy.deepcopy(operator_card) - nt, no = rc.update(tc, oc) - assert nt.heavy.intrinsic_flavors == [4] - for _, scale in nt.heavy.masses: - assert np.isnan(scale) - m2s = rc.masses(nt, no.configs.evolution_method) - raw = rc.Legacy.heavies("m%s", tc) - raw2 = np.power(raw, 2.0) - np.testing.assert_allclose(m2s, raw2) - tc["HQ"] = "MSBAR" - tc["Qmc"] = raw[0] * 1.1 - tc["Qmb"] = raw[1] * 1.1 - tc["Qmt"] = raw[2] * 0.9 - nt, no = rc.update(tc, oc) - for _, scale in nt.heavy.masses: - assert not np.isnan(scale) - m2s = rc.masses(nt, no.configs.evolution_method) - for m1, m2 in zip(m2s, raw2): - assert not np.isclose(m1, m2) - tc["HQ"] = "Blub" - with pytest.raises(ValueError, match="mass scheme"): - _nt, _no = rc.update(tc, oc) - nt.heavy.masses_scheme = "Bla" - with pytest.raises(ValueError, match="mass scheme"): - _ms = rc.masses(nt, no.configs.evolution_method) - - def test_legacy_fallback(): assert rc.Legacy.fallback(1, 2, 3) == 1 assert rc.Legacy.fallback(None, 2, 3) == 2 diff --git a/tests/eko/kernels/test_kernels_QEDns.py b/tests/eko/kernels/test_kernels_QEDns.py index 9fce3097b..bfa655a81 100644 --- a/tests/eko/kernels/test_kernels_QEDns.py +++ b/tests/eko/kernels/test_kernels_QEDns.py @@ -120,9 +120,7 @@ def test_zero_true_gamma(): dict( alphas=alpharef[0], alphaem=alpharef[1], - scale=muref, - num_flavs_ref=5, - max_num_flavs=6, + ref=(muref, 5), ) ) evmod = CouplingEvolutionMethod.EXACT diff --git a/tests/eko/quantities/test_couplings.py b/tests/eko/quantities/test_couplings.py index 00ccae9d1..e4d194d5b 100644 --- a/tests/eko/quantities/test_couplings.py +++ b/tests/eko/quantities/test_couplings.py @@ -3,7 +3,7 @@ def test_couplings_ref(): scale = 90.0 - d = dict(alphas=0.1, alphaem=0.01, scale=scale, max_num_flavs=6, num_flavs_ref=5) + d = dict(alphas=0.1, alphaem=0.01, ref=(scale, 5)) couplings = CouplingsInfo.from_dict(d) - assert couplings.scale == scale + assert couplings.ref[0] == scale assert not couplings.em_running diff --git a/tests/eko/runner/test_legacy.py b/tests/eko/runner/test_legacy.py index 60395acc5..1828285c3 100644 --- a/tests/eko/runner/test_legacy.py +++ b/tests/eko/runner/test_legacy.py @@ -36,7 +36,7 @@ class FakeEM(enum.Enum): eko.runner.legacy.Runner(theory_card, operator_card, path=path) # MSbar scheme theory_card.heavy.masses_scheme = QuarkMassScheme.MSBAR - theory_card.couplings.num_flavs_ref = 5 + theory_card.couplings.ref = (91.0, 5) theory_card.heavy.masses.c.scale = 2 theory_card.heavy.masses.b.scale = 4.5 theory_card.heavy.masses.t.scale = 173.07 diff --git a/tests/eko/scale_variations/test_diff.py b/tests/eko/scale_variations/test_diff.py index 8ade8bcac..13cc4938f 100644 --- a/tests/eko/scale_variations/test_diff.py +++ b/tests/eko/scale_variations/test_diff.py @@ -27,9 +27,7 @@ def compute_a_s(q2, order): couplings=CouplingsInfo( alphas=0.1181, alphaem=0.007496, - scale=91.00, - max_num_flavs=4, - num_flavs_ref=4, + ref=(91.00, 4), ), order=order, method=CouplingEvolutionMethod.EXPANDED, diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index bfacc1c96..309a87597 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/tests/eko/scale_variations/test_expanded.py @@ -18,9 +18,7 @@ def compute_a_s(q2, order): couplings=CouplingsInfo( alphas=0.1181, alphaem=0.007496, - scale=91.00, - max_num_flavs=4, - num_flavs_ref=4, + ref=(91.00, 4), ), order=order, method=CouplingEvolutionMethod.EXPANDED, diff --git a/tests/eko/test_couplings.py b/tests/eko/test_couplings.py index f8e051e4c..9b68d308b 100644 --- a/tests/eko/test_couplings.py +++ b/tests/eko/test_couplings.py @@ -54,9 +54,7 @@ def test_init(self): dict( alphas=alpharef[0], alphaem=alpharef[1], - scale=muref, - num_flavs_ref=5, - max_num_flavs=6, + ref=(muref, 5), ) ) order = (1, 0) @@ -109,7 +107,7 @@ def test_init(self): ) with pytest.raises(ValueError): coup3 = copy.deepcopy(couplings) - coup3.scale = 0 + coup3.ref = (0.0, 5) Couplings( coup3, order, @@ -161,9 +159,7 @@ def test_ref(self): dict( alphas=alpharef[0], alphaem=alpharef[1], - scale=muref, - num_flavs_ref=nfref, - max_num_flavs=6, + ref=(muref, nfref), ) ) for order_s in [1, 2, 3, 4]: @@ -196,9 +192,7 @@ def test_ref_copy_e841b0dfdee2f31d9ccc1ecee4d9d1a6f6624313(self): dict( alphas=alpharef[0], alphaem=alpharef[1], - scale=muref, - num_flavs_ref=3, # reference nf is needed to force the matching - max_num_flavs=6, + ref=(muref, 3), # reference nf is needed to force the matching ) ) sc = Couplings( @@ -237,9 +231,7 @@ def test_exact(self): dict( alphas=alpharef[0], alphaem=alpharef[1], - scale=muref, - num_flavs_ref=nfref, - max_num_flavs=6, + ref=(muref, nfref), em_running=em_running, ) ) @@ -300,9 +292,7 @@ def benchmark_expanded_n3lo(self): couplings = CouplingsInfo( alphas=alpharef[0], alphaem=alpharef[1], - scale=muref, - num_flavs_ref=5, - max_num_flavs=6, + ref=(muref, 5), ) m2c = 2 m2b = 25 diff --git a/tests/eko/test_msbar_masses.py b/tests/eko/test_msbar_masses.py index 12cb88f91..25718e962 100644 --- a/tests/eko/test_msbar_masses.py +++ b/tests/eko/test_msbar_masses.py @@ -16,9 +16,8 @@ def theory_card(theory_card: TheoryCard): th = theory_card th.order = (3, 0) th.couplings.alphas = 0.1180 - th.couplings.scale = 91.0 th.couplings.alphaem = 0.00781 - th.couplings.num_flavs_ref = 5 + th.couplings.ref = (91.0, 5) th.heavy.masses = HeavyQuarkMasses( [QuarkMassRef(val) for val in [(2.0, 2.1), (4.0, 4.1), (175.0, 174.9)]] ) diff --git a/tests/eko/test_quantities.py b/tests/eko/test_quantities.py index 0c2628b30..d3c412a25 100644 --- a/tests/eko/test_quantities.py +++ b/tests/eko/test_quantities.py @@ -33,9 +33,6 @@ def test_HeavyQuarks(): def test_HeavyInfo(): i = hq.HeavyInfo( - num_flavs_init=4, - num_flavs_max_pdf=6, - intrinsic_flavors=[4, 5], masses=hq.HeavyQuarkMasses( [ hq.QuarkMassRef([2.0, nan]), diff --git a/tests/ekobox/test_cards.py b/tests/ekobox/test_cards.py index 1a35e137e..00ef0e4c3 100644 --- a/tests/ekobox/test_cards.py +++ b/tests/ekobox/test_cards.py @@ -10,14 +10,14 @@ def test_generate_ocard(): mu0 = 1.65 mugrid = [(10.0, 6), (100.0, 5)] op = cards.example.operator() - op.mu0 = mu0 + op.init = (mu0, 4) op.mugrid = mugrid assert pytest.approx(op.mugrid) == mugrid assert pytest.approx(op.mu2grid) == np.array([mu**2 for mu, _ in mugrid]) assert op.configs.interpolation_polynomial_degree == 4 mugrid1 = [100.0, 5] op = cards.example.operator() - op.mu0 = mu0 + op.init = (mu0, 4) op.mugrid = mugrid1 op.configs.interpolation_polynomial_degree = 2 op.configs.interpolation_is_log = False @@ -35,7 +35,7 @@ def test_dump_load_op_card(tmp_path, cd): path1 = tmp_path / "debug_op.yaml" path2 = tmp_path / "debug_op_two.yaml" op = cards.example.operator() - op.mu0 = mu0 + op.init = (mu0, 4) cards.dump(op.raw, path1) cards.dump(op.raw, path2) op_loaded = cards.load(path1) @@ -52,7 +52,7 @@ def test_generate_theory_card(): assert theory.order[0] == 2 rawt = cards.example.raw_theory() assert isinstance(rawt, dict) - assert theory.heavy.num_flavs_init == rawt["heavy"]["num_flavs_init"] + assert theory.heavy.masses.c[0] == rawt["heavy"]["masses"][0][0] def containsnan(obj) -> bool: diff --git a/tests/ekobox/test_evol_pdf.py b/tests/ekobox/test_evol_pdf.py index 141787b49..7d788277b 100644 --- a/tests/ekobox/test_evol_pdf.py +++ b/tests/ekobox/test_evol_pdf.py @@ -12,7 +12,7 @@ def init_cards(): op = cards.example.operator() - op.mu0 = 1.65 + op.init = (1.65, 4) op.xgrid = XGrid([0.1, 0.5, 1.0]) op.configs.interpolation_polynomial_degree = 1 theory = cards.example.theory() diff --git a/tests/ekobox/test_info_file.py b/tests/ekobox/test_info_file.py index acf3f7272..62df7d739 100644 --- a/tests/ekobox/test_info_file.py +++ b/tests/ekobox/test_info_file.py @@ -10,7 +10,7 @@ def test_build(): theory.order = (2, 0) theory.couplings.alphas = 0.2 op = cards.example.operator() - op.mu0 = 1.0 + op.init = (1.0, 3) op.mugrid = [(10.0, 5), (100.0, 5)] info = info_file.build( theory, op, 4, info_update={"SetDesc": "Prova", "NewArg": 15.3, "MTop": 1.0} diff --git a/tests/ekobox/test_utils.py b/tests/ekobox/test_utils.py index 557db679c..1233d1d91 100644 --- a/tests/ekobox/test_utils.py +++ b/tests/ekobox/test_utils.py @@ -16,10 +16,9 @@ def test_ekos_product(tmp_path): theory = cards.example.theory() theory.order = (1, 0) - theory.heavy.num_flavs_init = 5 op1 = cards.example.operator() - op1.mu0 = mu01 + op1.init = (mu01, 5) op1.mugrid = mugrid1 op1.xgrid = xgrid op1.configs.interpolation_polynomial_degree = 1 @@ -28,13 +27,13 @@ def test_ekos_product(tmp_path): mugrid2 = [(8.0, 5), (10.0, 5), (12.0, 5)] op2 = cards.example.operator() - op2.mu0 = mu0 + op2.init = (mu0, 5) op2.mugrid = mugrid2 op2.xgrid = xgrid op2.configs.interpolation_polynomial_degree = 1 op_err = copy.deepcopy(op2) - op_err.mu0 = mu01 + op_err.init = (mu01, 5) mu2first = (mugrid2[0][0] ** 2, mugrid2[0][1]) diff --git a/tests/ekomark/data/__init__.py b/tests/ekomark/data/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/ekomark/data/test_init.py b/tests/ekomark/data/test_init.py new file mode 100644 index 000000000..278a06181 --- /dev/null +++ b/tests/ekomark/data/test_init.py @@ -0,0 +1,96 @@ +import copy + +import numpy as np +import pytest +from banana.data.theories import default_card as theory_card + +from eko.io import runcards as rc +from ekomark.data import update_runcards +from ekomark.data.operators import default_card as operator_card + +from ...eko.io.test_runcards import check_dumpable + + +def test_runcards_opcard(): + # check conversion + tc = copy.deepcopy(theory_card) + oc = copy.deepcopy(operator_card) + tc["Q0"] = 2.0 + # mugrid + oc["mugrid"] = [2.0, 10.0] + _nt, no = update_runcards(tc, oc) + assert isinstance(no, rc.OperatorCard) + assert len(no.evolgrid) == len(oc["mugrid"]) + assert len(no.mu2grid) == len(no.evolgrid) + assert no.evolgrid[0][-1] == 4 + assert no.evolgrid[1][-1] == 5 + np.testing.assert_allclose(no.init[0], tc["Q0"]) + np.testing.assert_allclose(no.mu20, tc["Q0"] ** 2.0) + assert len(no.pids) == 14 + check_dumpable(no) + del oc["mugrid"] + # or mu2grid + oc["mu2grid"] = [9.0, 30.0, 32.0] + _nt, no = update_runcards(tc, oc) + assert isinstance(no, rc.OperatorCard) + assert len(no.evolgrid) == len(oc["mu2grid"]) + assert len(no.mu2grid) == len(no.evolgrid) + assert no.evolgrid[0][-1] == 4 + assert no.evolgrid[1][-1] == 5 + assert no.evolgrid[2][-1] == 5 + check_dumpable(no) + del oc["mu2grid"] + # or Q2grid + oc["Q2grid"] = [15.0, 130.0, 140.0, 1e5] + _nt, no = update_runcards(tc, oc) + assert isinstance(no, rc.OperatorCard) + assert len(no.evolgrid) == len(oc["Q2grid"]) + assert len(no.mu2grid) == len(no.evolgrid) + assert no.evolgrid[0][-1] == 4 + assert no.evolgrid[1][-1] == 5 + assert no.evolgrid[2][-1] == 5 + assert no.evolgrid[3][-1] == 6 + check_dumpable(no) + + +def test_runcards_ekomark(): + # check conversion + tc = copy.deepcopy(theory_card) + oc = copy.deepcopy(operator_card) + nt, no = update_runcards(tc, oc) + assert isinstance(nt, rc.TheoryCard) + assert isinstance(no, rc.OperatorCard) + # check is idempotent + nnt, nno = update_runcards(nt, no) + assert nnt == nt + assert nno == no + + +def test_runcards_quarkmass(): + tc = copy.deepcopy(theory_card) + tc["nfref"] = 5 + tc["IC"] = 1 + oc = copy.deepcopy(operator_card) + nt, no = update_runcards(tc, oc) + for _, scale in nt.heavy.masses: + assert np.isnan(scale) + m2s = rc.masses(nt, no.configs.evolution_method) + raw = rc.Legacy.heavies("m%s", tc) + raw2 = np.power(raw, 2.0) + np.testing.assert_allclose(m2s, raw2) + tc["HQ"] = "MSBAR" + tc["Qmc"] = raw[0] * 1.1 + tc["Qmb"] = raw[1] * 1.1 + tc["Qmt"] = raw[2] * 0.9 + nt, no = update_runcards(tc, oc) + for _, scale in nt.heavy.masses: + assert not np.isnan(scale) + m2s = rc.masses(nt, no.configs.evolution_method) + for m1, m2 in zip(m2s, raw2): + assert not np.isclose(m1, m2) + tc["HQ"] = "Blub" + with pytest.raises(ValueError, match="mass scheme"): + _nt, _no = update_runcards(tc, oc) + nt.heavy.masses_scheme = "Bla" + with pytest.raises(ValueError, match="mass scheme"): + _ms = rc.masses(nt, no.configs.evolution_method)