Fix FK computation whenmur != muf

benchmarks/bench_evolve.py

@@ -98,9 +98,11 @@ def benchmark_dglap(tmp_path, test_files, test_configs):
 
 
 def benchmark_evolve_grid(tmp_path, lhapdf_path, test_files, test_pdf):
-    pine_path = test_files / "data/grids/400/HERA_NC_225GEV_EP_SIGMARED.pineappl.lz4"
+    tid = 400
+    tcard = pineko.theory_card.load(tid)
+    pine_path = test_files / f"data/grids/{tid}/HERA_NC_225GEV_EP_SIGMARED.pineappl.lz4"
     pinegrid = pineappl.grid.Grid.read(pine_path)
-    eko_path = test_files / "data/ekos/400/HERA_NC_225GEV_EP_SIGMARED.tar"
+    eko_path = test_files / f"data/ekos/{tid}/HERA_NC_225GEV_EP_SIGMARED.tar"
     target_path = pathlib.Path(tmp_path / "test_fktable.pineappl.lz4")
     max_as = 3
     max_al = 0
@@ -119,6 +121,7 @@ def benchmark_evolve_grid(tmp_path, lhapdf_path, test_files, test_pdf):
                 xir=1.0,
                 xif=1.0,
                 xia=1.0,
+                theory_meta=tcard,
                 assumptions=assumptions,
                 comparison_pdfs=["NNPDF40_nnlo_as_01180"],
             )

src/pineko/cli/convolve.py

@@ -101,6 +101,7 @@ def subcommand(
             xir,
             xif,
             xia,
+            theory_meta={},  # TODO: fix or remove CLI completely
             assumptions=assumptions,
             comparison_pdfs=pdfs,
             min_as=min_as,

src/pineko/evolve.py

@@ -34,7 +34,6 @@ def sv_scheme(tcard):
     ----------
     tcard : dict
         theory card
-
     """
     modsv_list = {a.value for a in ScaleVariationsMethod}
     xif = tcard["XIF"]
@@ -49,6 +48,25 @@ def sv_scheme(tcard):
     return modsv
 
 
+def construct_atlas(tcard):
+    """Construct the atlas for heavy quarks matching.
+
+    Parameters
+    ----------
+    tcard : dict
+        theory card
+    """
+    masses = np.array([tcard["mc"], tcard["mb"], tcard["mt"]]) ** 2
+    thresholds_ratios = np.array([tcard["kcThr"], tcard["kbThr"], tcard["ktThr"]]) ** 2
+    for q in range(tcard["MaxNfPdf"] + 1, 6 + 1):
+        thresholds_ratios[q - 4] = np.inf
+    atlas = Atlas(
+        matching_scales=heavy_quarks.MatchingScales(masses * thresholds_ratios),
+        origin=(tcard["Q0"] ** 2, tcard["nf0"]),
+    )
+    return atlas
+
+
 def get_convolution_suffix(convolution: pineappl.convolutions.Conv) -> str:
     """Get the correct suffix for a given convolution.
 
@@ -200,19 +218,12 @@ def write_operator_card(
         raise ValueError("Template declares a value of Q0, nf0 different from theory")
 
     q2_grid = (xif * xif * muf2_grid).tolist()
-    masses = np.array([tcard["mc"], tcard["mb"], tcard["mt"]]) ** 2
-    thresholds_ratios = np.array([tcard["kcThr"], tcard["kbThr"], tcard["ktThr"]]) ** 2
-    for q in range(tcard["MaxNfPdf"] + 1, 6 + 1):
-        thresholds_ratios[q - 4] = np.inf
-    atlas = Atlas(
-        matching_scales=heavy_quarks.MatchingScales(masses * thresholds_ratios),
-        origin=(tcard["Q0"] ** 2, tcard["nf0"]),
-    )
     # If we are producing nFONLL FKs we need to look to NfFF...
     if check.is_num_fonll(tcard["FNS"]):
         nf = tcard["NfFF"]
         operators_card["mugrid"] = [(float(np.sqrt(q2)), int(nf)) for q2 in q2_grid]
     else:
+        atlas = construct_atlas(tcard)
         operators_card["mugrid"] = [
             (float(np.sqrt(q2)), nf_default(q2, atlas)) for q2 in q2_grid
         ]
@@ -287,9 +298,9 @@ def evolve_grid(
     xir: float,
     xif: float,
     xia: float,
+    theory_meta: dict,
     assumptions="Nf6Ind",
     comparison_pdfs: Optional[list[str]] = None,
-    meta_data=None,
     min_as=None,
 ):
     """Convolute grid with EKO from file paths.
@@ -312,12 +323,12 @@ def evolve_grid(
         factorization scale variation
     xia : float
         fragmentation scale variation
+    tcard: dict
+        card containing the theory parameters
     assumptions : str
         assumptions on the flavor dimension
     comparison_pdfs : list(str) or None
         if given, a comparison table (with / without evolution) will be printed
-    meta_data : None or dict
-        if given, additional meta data written to the FK table
     min_as: None or int
         minimum power of strong coupling
     """
@@ -335,6 +346,7 @@ def evolve_grid(
     if "integrability_version" in grid.metadata:
         x_grid = np.append(x_grid, 1.0)
 
+    muf2_grid = evol_info.fac1
     mur2_grid = evol_info.ren1
     xif = 1.0 if operators[0].operator_card.configs.scvar_method is not None else xif
     tcard = operators[0].theory_card
@@ -358,7 +370,17 @@ def evolve_grid(
     ren_grid2 = xir * xir * mur2_grid
     # NOTE: Currently, getting `nfgrid` from the first Operator is correct but this
     # might need to be addressed in the future
-    nfgrid = [x[1] for x in operators[0].operator_card.mugrid]
+    # TODO: Find a suitable way to avoid the following duplication. This should ideally
+    # be part of the operator card.
+    if not np.array_equal(muf2_grid, mur2_grid):
+        if check.is_num_fonll(theory_meta["FNS"]):
+            nfgrid = [int(theory_meta["NfFF"]) for _ in mur2_grid]
+        else:
+            q2mur_grid = (xir * xir * mur2_grid).tolist()
+            atlas = construct_atlas(theory_meta)
+            nfgrid = [nf_default(q2, atlas) for q2 in q2mur_grid]
+    else:
+        nfgrid = [x[1] for x in operators[0].operator_card.mugrid]
     alphas_values = [
         4.0 * np.pi * sc.a_s(mur2, nf_to=nf) for mur2, nf in zip(ren_grid2, nfgrid)
     ]
@@ -415,9 +437,7 @@ def prepare(operator, convolution_types):
         )
 
     fktable.set_metadata("pineko_version", version.__version__)
-    if meta_data is not None:
-        for k, v in meta_data.items():
-            fktable.set_metadata(k, v)
+    fktable.set_metadata("theory_card", json.dumps(theory_meta))
 
     # compare before/after
     comparison = None

src/pineko/theory.py

@@ -6,7 +6,6 @@
 commonly referred to as 'theory'.
 """
 
-import json
 import logging
 import time
 
@@ -563,9 +562,9 @@ def fk(self, name, grid_path, tcard, pdfs):
                 xir=xir,
                 xif=xif,
                 xia=xia,
+                theory_meta=tcard,
                 assumptions=assumptions,
                 comparison_pdfs=pdfs,
-                meta_data={"theory_card": json.dumps(tcard)},
             )
 
             if n_ekos > 1: