Merge pull request #239 from NNPDF/fix-mur-neq-muf

Fix FK computation when`mur != muf`

Author: scarlehoff

benchmarks/bench_cli.py

@@ -86,16 +86,27 @@ def benchmark_compare_cli(lhapdf_path, test_files, test_pdf):
     assert "right" in result.output
 
 
-def benchmark_convolve_cli(test_files, tmp_path):
+def benchmark_convolve_cli(test_files, tmp_path, test_configs):
+    theoryid = 400
     grid_path = pathlib.Path(
-        test_files / "data/grids/400/HERA_NC_225GEV_EP_SIGMARED.pineappl.lz4"
+        test_files / f"data/grids/{theoryid}/HERA_NC_225GEV_EP_SIGMARED.pineappl.lz4"
+    )
+    eko_path = pathlib.Path(
+        test_files / f"data/ekos/{theoryid}/HERA_NC_225GEV_EP_SIGMARED.tar"
     )
-    eko_path = pathlib.Path(test_files / "data/ekos/400/HERA_NC_225GEV_EP_SIGMARED.tar")
     fk_path = tmp_path / "testfk.pineappl.lz4"
     runner = CliRunner()
     result = runner.invoke(
         command,
-        ["convolve", str(fk_path), str(grid_path), "2", "0", str(eko_path)],
+        [
+            "convolve",
+            str(fk_path),
+            str(grid_path),
+            str(theoryid),
+            "2",
+            "0",
+            str(eko_path),
+        ],
     )
     assert "Optimizing for Nf6Ind" in result.output
 

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"],
             )

poetry.lock

@@ -28,7 +28,7 @@ packages = [{ include = "pineko", from = "src" }]
 [tool.poetry.dependencies]
 python = ">=3.9,<3.13"
 eko = "^0.15.0"
-pineappl = "^1.0.0"
+pineappl = "^1.3.0"
 PyYAML = "^6.0"
 numpy = "^1.21.0"
 pandas = "^2.1"
@@ -59,7 +59,7 @@ pytest-cov = "^4.0.0"
 pytest-env = "^0.6.2"
 pylint = "^3.1.0"
 banana-hep = "^0.6.13"
-pineappl-cli = "^1.0.0"
+pineappl-cli = "^1.3.0"
 
 [tool.poetry.group.dev.dependencies]
 pdbpp = "^0.11.6"

pyproject.toml

@@ -7,13 +7,14 @@
 import rich
 import rich_click as click
 
-from .. import evolve
+from .. import evolve, theory_card
 from ._base import command
 
 
 @command.command("convolve")
 @click.argument("fktable", type=click.Path())
 @click.argument("grid_path", type=click.Path(exists=True))
+@click.argument("theoryid", type=int)
 @click.argument("max_as", type=int)
 @click.argument("max_al", type=int)
 @click.argument("op_paths", type=click.Path(exists=True), nargs=-1)
@@ -37,6 +38,7 @@
 def subcommand(
     fktable,
     grid_path,
+    theoryid,
     max_as,
     max_al,
     op_paths,
@@ -73,6 +75,7 @@ def subcommand(
     grid = pineappl.grid.Grid.read(grid_path)
     grid.optimize()
     n_ekos = len(op_paths)
+    tcard = theory_card.load(theoryid)
     with eko.EKO.edit(pathlib.Path(op_paths[0])) as first_operator:
         operators = [first_operator]
         path_operators = f"[+] {op_paths[0]}\n"
@@ -101,6 +104,7 @@ def subcommand(
             xir,
             xif,
             xia,
+            theory_meta=tcard,
             assumptions=assumptions,
             comparison_pdfs=pdfs,
             min_as=min_as,

src/pineko/cli/convolve.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
@@ -356,9 +368,12 @@ def evolve_grid(
     # To compute the alphas values we are first reverting the factorization scale shift
     # and then obtaining the renormalization scale using xir.
     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]
+    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]
     alphas_values = [
         4.0 * np.pi * sc.a_s(mur2, nf_to=nf) for mur2, nf in zip(ren_grid2, nfgrid)
     ]
@@ -415,9 +430,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/evolve.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: