Factorize masses computation

src/eko/io/runcards.py

@@ -13,9 +13,15 @@
 import numpy.typing as npt
 
 from .. import basis_rotation as br
-from .. import interpolation
+from .. import interpolation, msbar_masses
 from .. import version as vmod
-from ..quantities.heavy_quarks import HeavyQuarkMasses, MatchingScales, QuarkMassScheme
+from ..couplings import couplings_mod_ev
+from ..quantities.heavy_quarks import (
+    HeavyQuarkMasses,
+    MatchingRatios,
+    MatchingScales,
+    QuarkMassScheme,
+)
 from .dictlike import DictLike
 from .types import (
     CouplingsRef,
@@ -52,11 +58,35 @@ class TheoryCard(DictLike):
     """List of heavy quark masses."""
     quark_masses_scheme: QuarkMassScheme
     """Scheme used to specify heavy quark masses."""
-    matching: MatchingScales
+    matching: MatchingRatios
     """Matching scale of heavy quark masses"""
     xif: float
     """Ratio between factorization scale and process scale."""
 
+    @property
+    def matching_scales(self) -> MatchingScales:
+        """Compute matching scales."""
+        return np.power(list(iter(self.matching)), 2.0) * np.power(
+            list(iter(self.quark_masses)), 2.0
+        )
+
+
+def masses(theory: TheoryCard, evmeth: EvolutionMethod):
+    """Compute masses in the chosen scheme."""
+    if theory.quark_masses_scheme is QuarkMassScheme.MSBAR:
+        return msbar_masses.compute(
+            theory.quark_masses,
+            theory.couplings,
+            theory.order,
+            couplings_mod_ev(evmeth),
+            np.power(list(iter(theory.matching)), 2.0),
+            xif2=theory.xif**2,
+        ).tolist()
+    if theory.quark_masses_scheme is QuarkMassScheme.POLE:
+        return [mq.value**2 for mq in theory.quark_masses]
+
+    raise ValueError(f"Unknown mass scheme '{theory.quark_masses_scheme}'")
+
 
 @dataclass
 class Debug(DictLike):
@@ -80,6 +110,10 @@ class Configs(DictLike):
     """
     ev_op_iterations: int
     """Number of intervals in which to break the global path."""
+    scvar_method: Optional[ScaleVariationsMethod]
+    """"""
+    inversion_method: Optional[InversionMethod]
+    """Which method to use for backward matching conditions."""
     interpolation_polynomial_degree: int
     """Degree of elements of the intepolation polynomial basis."""
     interpolation_is_log: bool
@@ -90,10 +124,6 @@ class Configs(DictLike):
     """If `true` do polarized evolution."""
     time_like: bool
     """If `true` do time-like evolution."""
-    scvar_method: Optional[ScaleVariationsMethod]
-    """"""
-    inversion_method: Optional[InversionMethod]
-    """Which method to use for backward matching conditions."""
     n_integration_cores: int = 1
     """Number of cores used to parallelize integration."""
 

src/eko/quantities/heavy_quarks.py

@@ -3,7 +3,7 @@
 import enum
 
 from ..io.dictlike import DictLike
-from ..io.types import LinearScale, reference_running
+from ..io.types import LinearScale, SquaredScale, reference_running
 
 FLAVORS = "cbt"
 
@@ -40,7 +40,7 @@ def __getitem__(self, key: int):
 QuarkMass = LinearScale
 QuarkMassRef = reference_running(QuarkMass)
 HeavyQuarkMasses = _heavy_quark(QuarkMassRef)
-MatchingScale = LinearScale
+MatchingScale = SquaredScale
 MatchingScales = _heavy_quark(MatchingScale)
 MatchingRatio = float
 MatchingRatios = _heavy_quark(MatchingRatio)

src/eko/runner/legacy.py

@@ -5,12 +5,11 @@
 
 import numpy as np
 
-from .. import interpolation, msbar_masses
+from .. import interpolation
 from ..couplings import Couplings, couplings_mod_ev
 from ..evolution_operator.grid import OperatorGrid
 from ..io import EKO, Operator, runcards
 from ..io.types import RawCard
-from ..quantities.heavy_quarks import QuarkMassScheme
 from ..thresholds import ThresholdsAtlas
 from . import commons
 
@@ -62,20 +61,7 @@ def __init__(
         )
 
         # setup the Threshold path, compute masses if necessary
-        masses = None
-        if new_theory.quark_masses_scheme is QuarkMassScheme.MSBAR:
-            masses = msbar_masses.compute(
-                new_theory.quark_masses,
-                new_theory.couplings,
-                new_theory.order,
-                couplings_mod_ev(new_operator.configs.evolution_method),
-                np.power(list(iter(new_theory.matching)), 2.0),
-                xif2=new_theory.xif**2,
-            ).tolist()
-        elif new_theory.quark_masses_scheme is QuarkMassScheme.POLE:
-            masses = [mq.value**2 for mq in new_theory.quark_masses]
-        else:
-            raise ValueError(f"Unknown mass scheme '{new_theory.quark_masses_scheme}'")
+        masses = runcards.masses(new_theory, new_operator.configs.evolution_method)
 
         # call explicitly iter to explain the static analyzer that is an
         # iterable

src/eko/runner/recipes.py

@@ -3,8 +3,11 @@
 from dataclasses import dataclass
 
 from .. import EKO
+from .. import scale_variations as sv
+from ..io import runcards
 from ..io.dictlike import DictLike
 from ..io.types import SquaredScale
+from ..thresholds import ThresholdsAtlas
 
 
 @dataclass
@@ -33,3 +36,19 @@ class MatchingRecipe(Recipe):
 def create(eko: EKO):
     """Create all associated recipes."""
     _ = eko.theory_card.matching
+
+    masses = runcards.masses(
+        eko.theory_card, eko.operator_card.configs.evolution_method
+    )
+
+    tc = ThresholdsAtlas(
+        masses=masses,
+        q2_ref=eko.operator_card.mu20,
+        nf_ref=eko.theory_card.num_flavs_init,
+        thresholds_ratios=None,
+        max_nf=eko.theory_card.num_flavs_max_pdf,
+    )
+
+    for mu2 in eko.mu2grid:
+        expanded = eko.operator_card.configs.scvar_method is sv.Modes.expanded
+        mu2f = mu2 * eko.theory_card.xif**2 if expanded else mu2