diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index d579fb856..dfd0862f5 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -43,6 +43,13 @@ repos:
         args: ["--add-ignore=D107,D105"]
         additional_dependencies:
           - toml
+  - repo: https://github.com/pre-commit/mirrors-mypy
+    rev: v1.4.1
+    hooks:
+      - id: mypy
+        additional_dependencies: [types-PyYAML]
+        pass_filenames: false
+        args: ["--ignore-missing-imports", "src/"]
   - repo: local
     hooks:
     - id: fmt
diff --git a/src/eko/couplings.py b/src/eko/couplings.py
index 2b261afdb..a0826c09a 100644
--- a/src/eko/couplings.py
+++ b/src/eko/couplings.py
@@ -9,10 +9,11 @@
 """
 
 import logging
-from typing import Iterable, List
+from typing import Dict, Iterable, List, Tuple
 
 import numba as nb
 import numpy as np
+import numpy.typing as npt
 import scipy
 
 from . import constants, matchings
@@ -383,6 +384,10 @@ def couplings_expanded_fixed_alphaem(order, couplings_ref, nf, scale_from, scale
     return np.array([res_as, aem])
 
 
+_CouplingsCacheKey = Tuple[float, float, int, float, float]
+"""Cache key containing (a0, a1, nf, scale_from, scale_to)."""
+
+
 class Couplings:
     r"""Compute the strong and electromagnetic coupling constants :math:`a_s, a_{em}`.
 
@@ -480,7 +485,7 @@ def assert_positive(name, var):
                 self.decoupled_running,
             )
         # cache
-        self.cache = {}
+        self.cache: Dict[_CouplingsCacheKey, npt.NDArray] = {}
 
     @property
     def mu2_ref(self):
diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index ec56b6db6..fe07ade9e 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -8,6 +8,7 @@
 import os
 import time
 from multiprocessing import Pool
+from typing import Dict, Tuple
 
 import numba as nb
 import numpy as np
@@ -20,7 +21,7 @@
 from .. import basis_rotation as br
 from .. import interpolation, mellin
 from .. import scale_variations as sv
-from ..io.types import EvolutionMethod
+from ..io.types import EvolutionMethod, OperatorLabel
 from ..kernels import ev_method
 from ..kernels import non_singlet as ns
 from ..kernels import non_singlet_qed as qed_ns
@@ -602,6 +603,10 @@ def quad_ker_qed(
     return ker
 
 
+OpMembers = Dict[OperatorLabel, OpMember]
+"""Map of all operators."""
+
+
 class Operator(sv.ScaleVariationModeMixin):
     """Internal representation of a single EKO.
 
@@ -627,8 +632,8 @@ class Operator(sv.ScaleVariationModeMixin):
 
     log_label = "Evolution"
     # complete list of possible evolution operators labels
-    full_labels = br.full_labels
-    full_labels_qed = br.full_unified_labels
+    full_labels: Tuple[OperatorLabel, ...] = br.full_labels
+    full_labels_qed: Tuple[OperatorLabel, ...] = br.full_unified_labels
 
     def __init__(
         self, config, managers, segment: Segment, mellin_cut=5e-2, is_threshold=False
@@ -641,9 +646,9 @@ def __init__(
         # TODO make 'cut' external parameter?
         self._mellin_cut = mellin_cut
         self.is_threshold = is_threshold
-        self.op_members = {}
+        self.op_members: OpMembers = {}
         self.order = tuple(config["order"])
-        self.alphaem_running = self.managers["couplings"].alphaem_running
+        self.alphaem_running = self.managers.couplings.alphaem_running
         if self.log_label == "Evolution":
             self.a = self.compute_a()
             self.as_list, self.a_half_list = self.compute_aem_list()
@@ -665,7 +670,7 @@ def xif2(self):
     @property
     def int_disp(self):
         """Return the interpolation dispatcher."""
-        return self.managers["interpol_dispatcher"]
+        return self.managers.interpolator
 
     @property
     def grid_size(self):
@@ -688,7 +693,7 @@ def mu2(self):
 
     def compute_a(self):
         """Return the computed values for :math:`a_s` and :math:`a_{em}`."""
-        coupling = self.managers["couplings"]
+        coupling = self.managers.couplings
         a0 = coupling.a(
             self.mu2[0],
             nf_to=self.nf,
@@ -724,7 +729,7 @@ def compute_aem_list(self):
             as_list = np.array([self.a_s[0], self.a_s[1]])
             a_half = np.zeros((ev_op_iterations, 2))
         else:
-            couplings = self.managers["couplings"]
+            couplings = self.managers.couplings
             mu2_steps = np.geomspace(self.q2_from, self.q2_to, 1 + ev_op_iterations)
             mu2_l = mu2_steps[0]
             as_list = np.array(
diff --git a/src/eko/evolution_operator/__init__.py.patch b/src/eko/evolution_operator/__init__.py.patch
index d8258ee14..995d77b6b 100644
--- a/src/eko/evolution_operator/__init__.py.patch
+++ b/src/eko/evolution_operator/__init__.py.patch
@@ -1,8 +1,8 @@
 diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
-index ec56b6db..374d0d0b 100644
+index fe07ade9..0f58c9e5 100644
 --- a/src/eko/evolution_operator/__init__.py
 +++ b/src/eko/evolution_operator/__init__.py
-@@ -3,15 +3,15 @@ r"""Contains the central operator classes.
+@@ -3,16 +3,16 @@ r"""Contains the central operator classes.
  See :doc:`Operator overview </code/Operators>`.
  """
 
@@ -11,6 +11,7 @@ index ec56b6db..374d0d0b 100644
  import os
  import time
  from multiprocessing import Pool
+ from typing import Dict, Tuple
 
 +import ekors
  import numba as nb
@@ -20,7 +21,7 @@ index ec56b6db..374d0d0b 100644
 
  import ekore.anomalous_dimensions.polarized.space_like as ad_ps
  import ekore.anomalous_dimensions.unpolarized.space_like as ad_us
-@@ -29,92 +29,10 @@ from ..kernels import singlet_qed as qed_s
+@@ -30,92 +30,10 @@ from ..kernels import singlet_qed as qed_s
  from ..kernels import valence_qed as qed_v
  from ..matchings import Segment, lepton_number
  from ..member import OpMember
@@ -114,7 +115,7 @@ index ec56b6db..374d0d0b 100644
  spec = [
      ("is_singlet", nb.boolean),
      ("is_QEDsinglet", nb.boolean),
-@@ -186,422 +104,6 @@ class QuadKerBase:
+@@ -187,421 +105,6 @@ class QuadKerBase:
          return self.path.prefactor * pj * self.path.jac
 
 
@@ -533,11 +534,10 @@ index ec56b6db..374d0d0b 100644
 -            )
 -    return ker
 -
--
- class Operator(sv.ScaleVariationModeMixin):
-     """Internal representation of a single EKO.
 
-@@ -787,50 +289,6 @@ class Operator(sv.ScaleVariationModeMixin):
+ OpMembers = Dict[OperatorLabel, OpMember]
+ """Map of all operators."""
+@@ -792,50 +295,6 @@ class Operator(sv.ScaleVariationModeMixin):
          """Return the evolution method."""
          return ev_method(EvolutionMethod(self.config["method"]))
 
@@ -588,7 +588,7 @@ index ec56b6db..374d0d0b 100644
      def initialize_op_members(self):
          """Init all operators with the identity or zeros."""
          eye = OpMember(
-@@ -853,10 +311,7 @@ class Operator(sv.ScaleVariationModeMixin):
+@@ -858,10 +317,7 @@ class Operator(sv.ScaleVariationModeMixin):
              else:
                  self.op_members[n] = zero.copy()
 
@@ -600,7 +600,7 @@ index ec56b6db..374d0d0b 100644
          """Run the integration for each grid point.
 
          Parameters
-@@ -871,18 +326,56 @@ class Operator(sv.ScaleVariationModeMixin):
+@@ -876,18 +332,56 @@ class Operator(sv.ScaleVariationModeMixin):
          """
          column = []
          k, logx = log_grid
diff --git a/src/eko/evolution_operator/grid.py b/src/eko/evolution_operator/grid.py
index 2539c3fc7..6a9a3d9db 100644
--- a/src/eko/evolution_operator/grid.py
+++ b/src/eko/evolution_operator/grid.py
@@ -6,8 +6,8 @@
 """
 
 import logging
-from dataclasses import astuple
-from typing import Dict, List, Optional
+from dataclasses import dataclass
+from typing import Any, Dict, List, Optional
 
 import numpy as np
 import numpy.typing as npt
@@ -18,9 +18,9 @@
 from ..interpolation import InterpolatorDispatcher
 from ..io.runcards import Configs, Debug
 from ..io.types import EvolutionPoint as EPoint
-from ..io.types import Order
+from ..io.types import Order, SquaredScale
 from ..matchings import Atlas, Segment, flavor_shift, is_downward_path
-from . import Operator, flavors, matching_condition, physical
+from . import Operator, OpMembers, flavors, matching_condition, physical
 from .operator_matrix_element import OperatorMatrixElement
 
 logger = logging.getLogger(__name__)
@@ -29,6 +29,15 @@
 """In particular, only the ``operator`` and ``error`` fields are expected."""
 
 
+@dataclass(frozen=True)
+class Managers:
+    """Set of steering objects."""
+
+    atlas: Atlas
+    couplings: Couplings
+    interpolator: InterpolatorDispatcher
+
+
 class OperatorGrid(sv.ScaleVariationModeMixin):
     """Collection of evolution operators for several scales.
 
@@ -64,7 +73,7 @@ def __init__(
         use_fhmruvv: bool,
     ):
         # check
-        config = {}
+        config: Dict[str, Any] = {}
         config["order"] = order
         config["intrinsic_range"] = intrinsic_flavors
         config["xif2"] = xif**2
@@ -95,13 +104,13 @@ def __init__(
 
         self.config = config
         self.q2_grid = mu2grid
-        self.managers = dict(
-            thresholds_config=atlas,
+        self.managers = Managers(
+            atlas=atlas,
             couplings=couplings,
-            interpol_dispatcher=interpol_dispatcher,
+            interpolator=interpol_dispatcher,
         )
-        self._threshold_operators = {}
-        self._matching_operators = {}
+        self._threshold_operators: Dict[Segment, Operator] = {}
+        self._matching_operators: Dict[SquaredScale, OpMembers] = {}
 
     def get_threshold_operators(self, path: List[Segment]) -> List[Operator]:
         """Generate the threshold operators.
@@ -123,7 +132,6 @@ def get_threshold_operators(self, path: List[Segment]) -> List[Operator]:
         is_downward = is_downward_path(path)
         shift = flavor_shift(is_downward)
         for seg in path[:-1]:
-            new_op_key = astuple(seg)
             kthr = self.config["thresholds_ratios"][seg.nf - shift]
             ome = OperatorMatrixElement(
                 self.config,
@@ -134,13 +142,13 @@ def get_threshold_operators(self, path: List[Segment]) -> List[Operator]:
                 np.log(kthr),
                 self.config["HQ"] == "MSBAR",
             )
-            if new_op_key not in self._threshold_operators:
+            if seg not in self._threshold_operators:
                 # Compute the operator and store it
                 logger.info("Prepare threshold operator")
                 op_th = Operator(self.config, self.managers, seg, is_threshold=True)
                 op_th.compute()
-                self._threshold_operators[new_op_key] = op_th
-            thr_ops.append(self._threshold_operators[new_op_key])
+                self._threshold_operators[seg] = op_th
+            thr_ops.append(self._threshold_operators[seg])
 
             # Compute the matching conditions and store it
             if seg.target not in self._matching_operators:
@@ -159,7 +167,7 @@ def generate(self, q2: EPoint) -> OpDict:
 
         """
         # The lists of areas as produced by the thresholds
-        path = self.managers["thresholds_config"].path(q2)
+        path = self.managers.atlas.path(q2)
         # Prepare the path for the composition of the operator
         thr_ops = self.get_threshold_operators(path)
         # we start composing with the highest operator ...
diff --git a/src/eko/evolution_operator/operator_matrix_element.py b/src/eko/evolution_operator/operator_matrix_element.py
index 25e5df2dd..fc29aafc6 100644
--- a/src/eko/evolution_operator/operator_matrix_element.py
+++ b/src/eko/evolution_operator/operator_matrix_element.py
@@ -227,7 +227,7 @@ class OperatorMatrixElement(Operator):
 
     log_label = "Matching"
     # complete list of possible matching operators labels
-    full_labels = [
+    full_labels = (
         *br.singlet_labels,
         (br.matching_hplus_pid, 21),
         (br.matching_hplus_pid, 100),
@@ -238,7 +238,7 @@ class OperatorMatrixElement(Operator):
         (200, br.matching_hminus_pid),
         (br.matching_hminus_pid, 200),
         (br.matching_hminus_pid, br.matching_hminus_pid),
-    ]
+    )
     # still valid in QED since Sdelta and Vdelta matchings are diagonal
     full_labels_qed = copy.deepcopy(full_labels)
 
@@ -339,7 +339,7 @@ def a_s(self):
 
         Note that here you need to use :math:`a_s^{n_f+1}`
         """
-        sc = self.managers["couplings"]
+        sc = self.managers.couplings
         return sc.a_s(
             self.q2_from
             * (self.xif2 if self.sv_mode == sv.Modes.exponentiated else 1.0),
diff --git a/src/eko/io/inventory.py b/src/eko/io/inventory.py
index 5bb98842f..9b25ad75d 100644
--- a/src/eko/io/inventory.py
+++ b/src/eko/io/inventory.py
@@ -3,7 +3,7 @@
 import base64
 from dataclasses import asdict, dataclass, field
 from pathlib import Path
-from typing import Dict, Generic, Optional, Type, TypeVar
+from typing import Dict, Generic, Literal, Optional, Type, TypeVar
 
 import yaml
 
@@ -11,7 +11,7 @@
 from .items import Header, Operator
 
 NBYTES = 8
-ENDIANNESS = "little"
+ENDIANNESS: Literal["little", "big"] = "little"
 
 HEADER_EXT = ".yaml"
 ARRAY_EXT = [".npy", ".npz"]
diff --git a/src/eko/io/legacy.py b/src/eko/io/legacy.py
index c05d150ec..0458b6a78 100644
--- a/src/eko/io/legacy.py
+++ b/src/eko/io/legacy.py
@@ -13,15 +13,23 @@
 import numpy as np
 import yaml
 
-from eko.interpolation import XGrid
-from eko.io.runcards import flavored_mugrid
-from eko.quantities.heavy_quarks import HeavyInfo, HeavyQuarkMasses, MatchingRatios
-
+from ..interpolation import XGrid
+from ..io.runcards import flavored_mugrid
+from ..quantities.heavy_quarks import (
+    HeavyInfo,
+    HeavyQuarkMasses,
+    MatchingRatios,
+    QuarkMassScheme,
+)
 from . import raw
 from .dictlike import DictLike
 from .struct import EKO, Operator
 from .types import EvolutionPoint as EPoint
-from .types import RawCard
+from .types import RawCard, ReferenceRunning
+
+_MC = 1.51
+_MB = 4.92
+_MT = 172.5
 
 
 def load_tar(source: os.PathLike, dest: os.PathLike, errors: bool = False):
@@ -39,8 +47,8 @@ def load_tar(source: os.PathLike, dest: os.PathLike, errors: bool = False):
         whether to load also errors (default ``False``)
 
     """
-    with tempfile.TemporaryDirectory() as tmpdir:
-        tmpdir = pathlib.Path(tmpdir)
+    with tempfile.TemporaryDirectory() as tmpdirr:
+        tmpdir = pathlib.Path(tmpdirr)
 
         with tarfile.open(source, "r") as tar:
             raw.safe_extractall(tar, tmpdir)
@@ -60,7 +68,7 @@ def load_tar(source: os.PathLike, dest: os.PathLike, errors: bool = False):
     if op5 is None:
         op5 = metaold["mu2grid"]
     grid = op5to4(
-        flavored_mugrid(op5, theory.heavy.masses, theory.heavy.matching_ratios), arrays
+        flavored_mugrid(op5, [_MC, _MB, _MT], theory.heavy.matching_ratios), arrays
     )
 
     with EKO.create(dest) as builder:
@@ -93,10 +101,16 @@ def from_old(cls, old: RawCard):
         """Load from old metadata."""
         heavy = HeavyInfo(
             num_flavs_init=4,
-            num_flavs_max_pdf=None,
-            intrinsic_flavors=None,
-            masses=HeavyQuarkMasses([1.51, 4.92, 172.5]),
-            masses_scheme=None,
+            num_flavs_max_pdf=5,
+            intrinsic_flavors=[],
+            masses=HeavyQuarkMasses(
+                [
+                    ReferenceRunning([_MC, np.inf]),
+                    ReferenceRunning([_MB, np.inf]),
+                    ReferenceRunning([_MT, np.inf]),
+                ]
+            ),
+            masses_scheme=QuarkMassScheme.POLE,
             matching_ratios=MatchingRatios([1.0, 1.0, 1.0]),
         )
         return cls(heavy=heavy)
@@ -125,7 +139,7 @@ def from_old(cls, old: RawCard):
             mu2list = old["mu2grid"]
         mu2grid = np.array(mu2list)
         evolgrid = flavored_mugrid(
-            np.sqrt(mu2grid).tolist(), [1.51, 4.92, 172.5], [1, 1, 1]
+            np.sqrt(mu2grid).tolist(), [_MC, _MB, _MT], [1, 1, 1]
         )
 
         xgrid = XGrid(old["interpolation_xgrid"])
diff --git a/src/eko/io/struct.py b/src/eko/io/struct.py
index 57d2ae9df..076139d71 100644
--- a/src/eko/io/struct.py
+++ b/src/eko/io/struct.py
@@ -273,7 +273,7 @@ def approx(
         Raises
         ------
         ValueError
-            if multiple values are find in the neighbourhood
+            if multiple values are found in the neighbourhood
 
         """
         eps = np.array([ep_ for ep_ in self if ep_[1] == ep[1]])
@@ -281,7 +281,9 @@ def approx(
         close = eps[np.isclose(ep[0], mu2s, rtol=rtol, atol=atol)]
 
         if len(close) == 1:
-            return tuple(close[0])
+            found = close[0]
+            assert isinstance(found[0], float)
+            return (found[0], int(found[1]))
         if len(close) == 0:
             return None
         raise ValueError(f"Multiple values of Q2 have been found close to {ep}")
@@ -374,17 +376,17 @@ def open(cls, path: os.PathLike, mode="r"):
             raise ValueError(f"Unknown file mode: {mode}")
 
         tmpdir = pathlib.Path(tempfile.mkdtemp(prefix=TEMP_PREFIX))
-        if load:
-            cls.load(path, tmpdir)
-            metadata = Metadata.load(tmpdir)
-            opened = cls(
-                **inventories(tmpdir, access),
-                metadata=metadata,
-                access=access,
-            )
-            opened.operators.sync()
-        else:
-            opened = Builder(path=tmpdir, access=access)
+        if not load:
+            return Builder(path=tmpdir, access=access)
+        # load existing instead
+        cls.load(path, tmpdir)
+        metadata = Metadata.load(tmpdir)
+        opened: EKO = cls(
+            **inventories(tmpdir, access),
+            metadata=metadata,
+            access=access,
+        )
+        opened.operators.sync()
 
         return opened
 
diff --git a/src/eko/io/types.py b/src/eko/io/types.py
index d8377a4e6..20ac52db5 100644
--- a/src/eko/io/types.py
+++ b/src/eko/io/types.py
@@ -1,8 +1,7 @@
 """Common type definitions, only used for static analysis."""
 
 import enum
-import typing
-from typing import Any, Dict, Generic, Tuple, TypeVar
+from typing import Any, Dict, Generic, List, Tuple, TypeVar
 
 # Energy scales
 # -------------
@@ -23,8 +22,9 @@
 Order = Tuple[int, int]
 FlavorsNumber = int
 FlavorIndex = int
-IntrinsicFlavors = typing.List[FlavorIndex]
-N3LOAdVariation = typing.Tuple[int, int, int, int]
+IntrinsicFlavors = List[FlavorIndex]
+N3LOAdVariation = Tuple[int, int, int, int]
+OperatorLabel = Tuple[int, int]
 
 # Evolution coordinates
 # ---------------------
diff --git a/src/eko/msbar_masses.py b/src/eko/msbar_masses.py
index 26fb9874e..de605764b 100644
--- a/src/eko/msbar_masses.py
+++ b/src/eko/msbar_masses.py
@@ -396,7 +396,7 @@ def sc(thr_masses):
     heavy_quarks = quark_names[3:]
     hq_idxs = np.arange(0, 3)
     if nf_ref > 4:
-        heavy_quarks = reversed(heavy_quarks)
+        heavy_quarks = "".join(reversed(heavy_quarks))
         hq_idxs = reversed(hq_idxs)
 
     # loop on heavy quarks and compute the msbar masses
diff --git a/src/eko/quantities/couplings.py b/src/eko/quantities/couplings.py
index aec1d083f..b223ce2a7 100644
--- a/src/eko/quantities/couplings.py
+++ b/src/eko/quantities/couplings.py
@@ -2,7 +2,6 @@
 
 import dataclasses
 import enum
-from typing import Optional
 
 from ..io.dictlike import DictLike
 from ..io.types import FlavorsNumber, LinearScale, ReferenceRunning, Scalar
@@ -23,7 +22,7 @@ class CouplingsInfo(DictLike):
     alphaem: Coupling
     scale: LinearScale
     max_num_flavs: FlavorsNumber
-    num_flavs_ref: Optional[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
diff --git a/src/eko/runner/operators.py b/src/eko/runner/operators.py
index 07ac34093..1efa8bb07 100644
--- a/src/eko/runner/operators.py
+++ b/src/eko/runner/operators.py
@@ -17,7 +17,9 @@ def retrieve(
     elements = []
     for head in headers:
         inv = parts if isinstance(head, Evolution) else parts_matching
-        elements.append(inv[head])
+        op = inv[head]
+        assert op is not None
+        elements.append(op)
 
     return elements
 
diff --git a/src/eko/runner/parts.py b/src/eko/runner/parts.py
index f9379bf76..798353f71 100644
--- a/src/eko/runner/parts.py
+++ b/src/eko/runner/parts.py
@@ -16,13 +16,14 @@
 from ..evolution_operator import matching_condition
 from ..evolution_operator import operator_matrix_element as ome
 from ..evolution_operator import physical
+from ..evolution_operator.grid import Managers
 from ..io import EKO
 from ..io.items import Evolution, Matching, Operator
 from ..quantities.heavy_quarks import QuarkMassScheme
 from . import commons
 
 
-def managers(eko: EKO) -> dict:
+def managers(eko: EKO) -> Managers:
     """Collect managers for operator computation.
 
     .. todo::
@@ -31,10 +32,10 @@ def managers(eko: EKO) -> dict:
     """
     tcard = eko.theory_card
     ocard = eko.operator_card
-    return dict(
-        thresholds_config=commons.atlas(tcard, ocard),
+    return Managers(
+        atlas=commons.atlas(tcard, ocard),
         couplings=commons.couplings(tcard, ocard),
-        interpol_dispatcher=commons.interpolator(ocard),
+        interpolator=commons.interpolator(ocard),
     )
 
 
diff --git a/src/eko/scale_variations/__init__.py b/src/eko/scale_variations/__init__.py
index 21d56a503..8d3f01556 100644
--- a/src/eko/scale_variations/__init__.py
+++ b/src/eko/scale_variations/__init__.py
@@ -5,6 +5,7 @@
 """
 
 import enum
+from typing import Any, Dict
 
 from ..io.types import ScaleVariationsMethod
 from . import expanded, exponentiated
@@ -40,6 +41,8 @@ def sv_mode(s: ScaleVariationsMethod) -> Modes:
 class ScaleVariationModeMixin:
     """Mixin to cast scale variation mode."""
 
+    config: Dict[str, Any]
+
     @property
     def sv_mode(self) -> Modes:
         """Return the scale variation mode."""
diff --git a/src/ekobox/apply.py b/src/ekobox/apply.py
index 336900798..a260a4cec 100644
--- a/src/ekobox/apply.py
+++ b/src/ekobox/apply.py
@@ -1,10 +1,14 @@
 """Apply operator evolution to PDF set."""
 
+from dataclasses import dataclass
+from typing import Dict, Optional, Union
+
 import numpy as np
 
 from eko import basis_rotation as br
 from eko import interpolation
 from eko.io import EKO
+from eko.io.types import EvolutionPoint
 
 
 def apply_pdf(
@@ -49,6 +53,17 @@ def apply_pdf(
 
 CONTRACTION = "ajbk,bk"
 
+_PdfLabel = Union[int, str]
+"""PDF identifier: either PID or label."""
+
+
+@dataclass
+class PdfResult:
+    """Helper class to collect PDF results."""
+
+    pdfs: Dict[_PdfLabel, float]
+    errors: Optional[Dict[_PdfLabel, float]] = None
+
 
 def apply_pdf_flavor(
     eko: EKO, lhapdf_like, targetgrid=None, flavor_rotation=None, labels=None
@@ -85,35 +100,34 @@ def apply_pdf_flavor(
         )
 
     # build output
-    out_grid = {}
+    out_grid: Dict[EvolutionPoint, PdfResult] = {}
     for ep, elem in eko.items():
         pdf_final = np.einsum(CONTRACTION, elem.operator, pdfs, optimize="optimal")
         if elem.error is not None:
             error_final = np.einsum(CONTRACTION, elem.error, pdfs, optimize="optimal")
         else:
             error_final = None
-        out_grid[ep] = {
-            "pdfs": dict(zip(eko.bases.targetpids, pdf_final)),
-            "errors": None,
-        }
+        out_grid[ep] = PdfResult(dict(zip(eko.bases.targetpids, pdf_final)))
         if error_final is not None:
-            out_grid[ep]["errors"] = dict(zip(eko.bases.targetpids, error_final))
+            out_grid[ep].errors = dict(zip(eko.bases.targetpids, error_final))
 
     qed = eko.theory_card.order[1] > 0
     # rotate to evolution basis
     if flavor_rotation is not None:
         for q2, op in out_grid.items():
             pdf = flavor_rotation @ np.array(
-                [op["pdfs"][pid] for pid in br.flavor_basis_pids]
+                [op.pdfs[pid] for pid in br.flavor_basis_pids]
             )
-            if labels is None:
-                labels = list(range(flavor_rotation.shape[0]))
-            op["pdfs"] = dict(zip(labels, pdf))
-            if op["errors"] is not None:
+            if not qed:
+                evol_basis = br.evol_basis
+            else:
+                evol_basis = br.unified_evol_basis
+            op.pdfs = dict(zip(evol_basis, pdf))
+            if op.errors is not None:
                 errors = flavor_rotation @ np.array(
-                    [op["errors"][pid] for pid in br.flavor_basis_pids]
+                    [op.errors[pid] for pid in br.flavor_basis_pids]
                 )
-                op["errors"] = dict(zip(labels, errors))
+                op.errors = dict(zip(evol_basis, errors))
 
     # rotate/interpolate to target grid
     if targetgrid is not None:
@@ -125,11 +139,12 @@ def apply_pdf_flavor(
 
         rot = b.get_interpolation(targetgrid)
         for evpdf in out_grid.values():
-            for pdf_label in evpdf["pdfs"]:
-                evpdf["pdfs"][pdf_label] = np.matmul(rot, evpdf["pdfs"][pdf_label])
-                if evpdf["errors"] is not None:
-                    evpdf["errors"][pdf_label] = np.matmul(
-                        rot, evpdf["errors"][pdf_label]
-                    )
-
-    return out_grid
+            for pdf_label in evpdf.pdfs:
+                evpdf.pdfs[pdf_label] = np.matmul(rot, evpdf.pdfs[pdf_label])
+                if evpdf.errors is not None:
+                    evpdf.errors[pdf_label] = np.matmul(rot, evpdf.errors[pdf_label])
+    # cast back to be backward compatible
+    real_out_grid = {}
+    for ep, res in out_grid.items():
+        real_out_grid[ep] = {"pdfs": res.pdfs, "errors": res.errors}
+    return real_out_grid
diff --git a/src/ekobox/cli/inspect.py b/src/ekobox/cli/inspect.py
index 5787c5f1e..c355fce61 100644
--- a/src/ekobox/cli/inspect.py
+++ b/src/ekobox/cli/inspect.py
@@ -28,7 +28,7 @@ def subcommand(ctx, path: pathlib.Path):
 @click.pass_obj
 def sub_mu2(operator: EKO):
     """Check operator's mu2grid."""
-    rich.print_json(data=operator.mu2grid.tolist())
+    rich.print_json(data=operator.mu2grid)
 
 
 @subcommand.command("cards")
diff --git a/src/ekobox/cli/run.py b/src/ekobox/cli/run.py
index 9a91e36ab..44cfe216e 100644
--- a/src/ekobox/cli/run.py
+++ b/src/ekobox/cli/run.py
@@ -53,11 +53,11 @@ def subcommand(paths: Sequence[pathlib.Path]):
     else:
         output = operator.parent / OUTPUT
 
-    theory = yaml.safe_load(theory.read_text(encoding="utf-8"))
-    if "order" in theory:
-        theory = TheoryCard.from_dict(theory)
-    operator = yaml.safe_load(operator.read_text(encoding="utf-8"))
-    if "configs" in operator:
-        operator = OperatorCard.from_dict(operator)
-
-    eko.solve(theory, operator, path=output)
+    tc = yaml.safe_load(theory.read_text(encoding="utf-8"))
+    if "order" in tc:
+        tc = TheoryCard.from_dict(tc)
+    oc = yaml.safe_load(operator.read_text(encoding="utf-8"))
+    if "configs" in oc:
+        oc = OperatorCard.from_dict(oc)
+
+    eko.solve(tc, oc, path=output)
diff --git a/src/ekobox/cli/runcards.py b/src/ekobox/cli/runcards.py
index 020a58917..227e012f4 100644
--- a/src/ekobox/cli/runcards.py
+++ b/src/ekobox/cli/runcards.py
@@ -3,6 +3,8 @@
 import logging
 import pathlib
 
+import numpy as np
+
 from .. import cards
 from . import library as lib
 from .base import command
@@ -35,6 +37,6 @@ def sub_example(destination: pathlib.Path):
     cards.dump(theory.raw, path=destination / "theory.yaml")
     operator = cards.example.operator()
     operator.mu0 = 1.65
-    operator.mu2grid = [1e5]
+    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/evol_pdf.py b/src/ekobox/evol_pdf.py
index abdb2ec6a..9da6dda65 100644
--- a/src/ekobox/evol_pdf.py
+++ b/src/ekobox/evol_pdf.py
@@ -73,12 +73,16 @@ def evolve_pdfs(
 
     # apply PDF to eko
     evolved_PDF_list = []
+    q2block_per_nf = {}
     with EKO.read(eko_path) as eko_output:
         for initial_PDF in initial_PDF_list:
             evolved_PDF_list.append(
                 apply.apply_pdf(eko_output, initial_PDF, targetgrid)
             )
 
+        # separate by nf the evolgrid (and order per nf/q)
+        q2block_per_nf = regroup_evolgrid(eko_output.evolgrid)  # pylint: disable=E1101
+
     # update info file
     if targetgrid is None:
         targetgrid = operators_card.xgrid
@@ -93,8 +97,6 @@ def evolve_pdfs(
         info_update=info_update,
     )
 
-    # in the eko scales are squared
-    q2block_per_nf = {nf: np.power(q2s, 2) for nf, q2s in q2block_per_nf.items()}
     # write all replicas
     all_member_blocks = []
     for evolved_PDF in evolved_PDF_list:
@@ -134,7 +136,7 @@ def pdf_xq2(pid, x, Q2):
             pdf_xq2,
             xgrid=xgrid,
             sorted_q2grid=q2grid,
-            pids=br.flavor_basis_pids,
+            pids=np.array(br.flavor_basis_pids),
         )
         all_blocks.append(block)
     return all_blocks
diff --git a/src/ekobox/genpdf/__init__.py b/src/ekobox/genpdf/__init__.py
index 4666891bc..d6c427989 100644
--- a/src/ekobox/genpdf/__init__.py
+++ b/src/ekobox/genpdf/__init__.py
@@ -17,7 +17,7 @@
 
 
 def take_data(
-    parent_pdf_set: Optional[Union[str, dict]] = None,
+    parent_pdf_set=None,
     members: bool = False,
     xgrid: Optional[List[float]] = None,
     evolgrid: Optional[List[EPoint]] = None,
@@ -63,7 +63,7 @@ def take_data(
             all_blocks.append(
                 [
                     generate_block(
-                        toylh.xfxQ2, xgrid, sorted_q2grid, br.flavor_basis_pids
+                        toylh.xfxQ2, xgrid, sorted_q2grid, list(br.flavor_basis_pids)
                     )
                 ]
             )
@@ -87,7 +87,7 @@ def take_data(
                     ),
                     xgrid,
                     sorted_q2grid,
-                    br.flavor_basis_pids,
+                    list(br.flavor_basis_pids),
                 )
             ]
         )
diff --git a/src/ekobox/info_file.py b/src/ekobox/info_file.py
index fc5c17f1d..56439599d 100644
--- a/src/ekobox/info_file.py
+++ b/src/ekobox/info_file.py
@@ -2,6 +2,7 @@
 
 import copy
 import math
+from typing import Any, Dict
 
 import numpy as np
 
@@ -81,7 +82,7 @@ def build_alphas(
         info file section in lhapdf format
     """
     # start with meta stuff
-    template_info = {}
+    template_info: Dict[str, Any] = {}
     template_info["AlphaS_MZ"] = theory_card.couplings.alphas
     template_info["AlphaS_OrderQCD"] = theory_card.order[0] - 1
     # prepare
diff --git a/src/ekomark/plots.py b/src/ekomark/plots.py
index 0395bd03d..15fa19192 100644
--- a/src/ekomark/plots.py
+++ b/src/ekomark/plots.py
@@ -2,6 +2,7 @@
 
 import io
 import pprint
+from typing import Dict, List, Tuple
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -231,7 +232,7 @@ def save_operators_to_pdf(path, theory, ops, me: EKO, skip_pdfs, change_lab=Fals
     """
     ops_names = list(me.bases.targetpids)
     if np.allclose(ops_names, br.rotate_flavor_to_evolution):
-        ops_names = br.evol_basis_pids
+        ops_names = list(br.evol_basis_pids)
     else:
         raise ValueError("Can not reconstruct PDF names")
     ops_id = f"o{ops['hash'][:6]}_t{theory['hash'][:6]}"
@@ -245,16 +246,16 @@ def save_operators_to_pdf(path, theory, ops, me: EKO, skip_pdfs, change_lab=Fals
 
         # plot the operators
         # it's necessary to reshuffle the eko output
-        for mu2 in me.mu2grid:
-            results = me[mu2].operator
-            errors = me[mu2].error
+        for ep, op in me.items():
+            results = op.operator
+            errors = op.error
 
             # loop on pids
             for label_out, res, res_err in zip(ops_names, results, errors):
                 if label_out in skip_pdfs:
                     continue
-                new_op = {}
-                new_op_err = {}
+                new_op: Dict[Tuple[int, ...], List[float]] = {}
+                new_op_err: Dict[Tuple[int, ...], List[float]] = {}
                 # loop on xgrid point
                 for j in range(len(me.xgrid)):
                     # loop on pid in
@@ -280,7 +281,7 @@ def save_operators_to_pdf(path, theory, ops, me: EKO, skip_pdfs, change_lab=Fals
                     fig = None
                     try:
                         fig = plot_operator(
-                            f"Operator ({lab_in};{lab_out}) µ_F^2 = {mu2} GeV^2",
+                            f"Operator ({lab_in};{lab_out}) µ_F^2 = {ep[0]} GeV^2, nf = {ep[1]}",
                             new_op[label_in],
                             new_op_err[label_in],
                         )
diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py
index 23395f66b..ed003f925 100644
--- a/tests/eko/evolution_operator/test_init.py
+++ b/tests/eko/evolution_operator/test_init.py
@@ -1,4 +1,5 @@
 import os
+from dataclasses import dataclass
 
 import numpy as np
 import pytest
@@ -227,7 +228,12 @@ def compute(self, a_ref, nf, scale_from, scale_to):
         return a_ref
 
 
-fake_managers = {"couplings": FakeCoupling()}
+@dataclass(frozen=True)
+class FakeManagers:
+    couplings: FakeCoupling
+
+
+fake_managers = FakeManagers(couplings=FakeCoupling())
 
 
 class TestOperator:
diff --git a/tests/eko/quantities/test_couplings.py b/tests/eko/quantities/test_couplings.py
index 0a3b97415..00ccae9d1 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=None)
+    d = dict(alphas=0.1, alphaem=0.01, scale=scale, max_num_flavs=6, num_flavs_ref=5)
     couplings = CouplingsInfo.from_dict(d)
     assert couplings.scale == scale
     assert not couplings.em_running
diff --git a/tests/eko/test_couplings.py b/tests/eko/test_couplings.py
index ef32fb0d0..f8e051e4c 100644
--- a/tests/eko/test_couplings.py
+++ b/tests/eko/test_couplings.py
@@ -55,7 +55,7 @@ def test_init(self):
                 alphas=alpharef[0],
                 alphaem=alpharef[1],
                 scale=muref,
-                num_flavs_ref=None,
+                num_flavs_ref=5,
                 max_num_flavs=6,
             )
         )
@@ -153,18 +153,19 @@ def test_ref(self):
             (0, np.inf, np.inf),
             (2, 4, 175),
         ]
+        nfrefs = (3, 4, 5)
         alpharef = (0.118, 0.00781)
         muref = 91.0
-        couplings = CouplingsInfo.from_dict(
-            dict(
-                alphas=alpharef[0],
-                alphaem=alpharef[1],
-                scale=muref,
-                num_flavs_ref=None,
-                max_num_flavs=6,
+        for thresh_setup, nfref in zip(thresh_setups, nfrefs):
+            couplings = CouplingsInfo.from_dict(
+                dict(
+                    alphas=alpharef[0],
+                    alphaem=alpharef[1],
+                    scale=muref,
+                    num_flavs_ref=nfref,
+                    max_num_flavs=6,
+                )
             )
-        )
-        for thresh_setup in thresh_setups:
             for order_s in [1, 2, 3, 4]:
                 for order_em in [0, 1, 2]:
                     for evmod in CouplingEvolutionMethod:
@@ -221,9 +222,10 @@ def test_exact(self):
             (0, np.inf, np.inf),
             (2, 4, 175),
         ]
+        nfrefs = (3, 4, 5)
         alpharef = np.array([0.118, 0.00781])
         muref = 91.0
-        for thresh_setup in thresh_setups:
+        for thresh_setup, nfref in zip(thresh_setups, nfrefs):
             for qcd in range(1, 4 + 1):
                 for qed in range(2 + 1):
                     for em_running in [
@@ -236,7 +238,7 @@ def test_exact(self):
                                 alphas=alpharef[0],
                                 alphaem=alpharef[1],
                                 scale=muref,
-                                num_flavs_ref=None,
+                                num_flavs_ref=nfref,
                                 max_num_flavs=6,
                                 em_running=em_running,
                             )
@@ -299,7 +301,7 @@ def benchmark_expanded_n3lo(self):
             alphas=alpharef[0],
             alphaem=alpharef[1],
             scale=muref,
-            num_flavs_ref=None,
+            num_flavs_ref=5,
             max_num_flavs=6,
         )
         m2c = 2