Numerical FONLL

src/pineko/check.py

@@ -116,6 +116,28 @@ def check_grid_and_eko_compatible(pineappl_grid, operators, xif, max_as, max_al)
         raise ValueError("x grid in pineappl grid and eko operator are NOT compatible!")
 
 
+def is_dis(lumi):
+    """Check if the fktable we are computing is a DIS fktable.
+
+    Parameters
+    ----------
+    lumi : list(list(tuple))
+           luminosity info
+
+    Returns
+    -------
+    bool
+        true if the fktable is a DIS fktable
+    """
+    for lists in lumi:
+        for el in lists:
+            if (not islepton(el[0])) and (not islepton(el[1])):
+                # If neither of the incoming particles is a lepton we are sure
+                # it is not DIS
+                return False
+    return True
+
+
 def is_fonll_b(fns, lumi):
     """Check if the fktable we are computing is a DIS FONLL-B fktable.
 
@@ -131,12 +153,7 @@ def is_fonll_b(fns, lumi):
     bool
         true if the fktable is a FONLL-B DIS fktable
     """
-    for lists in lumi:
-        for el in lists:
-            if (not islepton(el[0])) and (not islepton(el[1])):
-                # in this case we are sure it is not DIS so for sure it is not FONLL-B
-                return False
-    if fns == "FONLL-B":
+    if is_dis(lumi) and fns == "FONLL-B":
         return True
     return False
 

src/pineko/cli/__init__.py

@@ -1,3 +1,13 @@
 """CLI entry point."""
-from . import check, compare, convolute, gen_sv, kfactor, opcard, scaffold, theory_
+from . import (
+    check,
+    compare,
+    convolute,
+    fonll,
+    gen_sv,
+    kfactor,
+    opcard,
+    scaffold,
+    theory_,
+)
 from ._base import command

src/pineko/cli/fonll.py

@@ -0,0 +1,52 @@
+"""CLI entry point to FONLL."""
+import click
+
+from .. import configs, fonll, theory_card
+from ._base import command
+
+
+class TheoryCardError(Exception):
+    """Raised when asked for FONLL theory cards with an original tcard as input that is not asking for FONLL."""
+
+    pass
+
+
+class InconsistentInputsError(Exception):
+    """Raised if the inputs are not consistent with FONLL accounting for only charm mass, or both charm and bottom masses."""
+
+    pass
+
+
+@command.command("combine_fonll")
+@click.argument("FFNS3", type=click.Path(exists=True))
+@click.argument("FFN03", type=click.Path(exists=True))
+@click.argument("FFNS4til", type=click.Path(exists=True))
+@click.argument("theoryID", type=int)
+@click.option("--FFNS4bar", type=click.Path(exists=True))
+@click.option("--FFN04", type=click.Path(exists=True))
+@click.option("--FFNS5til", type=click.Path(exists=True))
+@click.option("--FFNS5bar", type=click.Path(exists=True))
+def subcommand(ffns3, ffn03, ffns4til, theoryid, ffns4bar, ffn04, ffns5til, ffns5bar):
+    """Combine the different FKs needed to produce the FONLL prescription."""
+    if (ffn04 is None and ffns5til is not None) or (
+        ffn04 is not None and ffns5til is None
+    ):
+        print(
+            "One between ffn04 and ffns5til has been provided without the other. Since they are both needed to construct FONLL, this does not make sense."
+        )
+        raise InconsistentInputsError
+    fonll.produce_combined_fk(
+        ffns3, ffn03, ffns4til, ffns4bar, ffn04, ffns5til, ffns5bar, theoryid
+    )
+
+
+@command.command("fonll_tcards")
+@click.argument("theoryID", type=int)
+def fonll_tcards(theoryid):
+    """Produce the FONLL tcards starting from the original tcard given by the theoryID."""
+    configs.configs = configs.defaults(configs.load())
+    tcard = theory_card.load(theoryid)
+    tcard_parent_path = theory_card.path(theoryid).parent
+    if "FONLL" not in tcard["FNS"]:
+        raise TheoryCardError
+    fonll.produce_fonll_tcards(tcard, tcard_parent_path)

src/pineko/evolve.py

@@ -106,6 +106,7 @@ def write_operator_card(pineappl_grid, default_card, card_path, tcard):
     """
     # Add a +1 to the orders for the difference in convention between nnpdf and pineappl
     is_fns = int(check.is_fonll_b(tcard["FNS"], pineappl_grid.lumi()))
+    # NB: This would not happen for simFONLL
     max_as = 1 + tcard["PTO"] + is_fns
     max_al = 1 + tcard["QED"]
     # ... in order to create a mask ...
@@ -128,9 +129,14 @@ def write_operator_card(pineappl_grid, default_card, card_path, tcard):
         matching_scales=heavy_quarks.MatchingScales(masses * thresholds_ratios),
         origin=(tcard["Q0"] ** 2, tcard["nf0"]),
     )
-    operators_card["mugrid"] = [
-        (float(np.sqrt(q2)), int(nf_default(q2, atlas))) for q2 in q2_grid
-    ]
+    # If we are producing simFONLL FKs we need to look to NfFF...
+    if issimFONLL(tcard["FNS"]):
+        nfff = tcard["NfFF"]
+        operators_card["mugrid"] = [(float(np.sqrt(q2)), int(nfff)) for q2 in q2_grid]
+    else:
+        operators_card["mugrid"] = [
+            (float(np.sqrt(q2)), int(nf_default(q2, atlas))) for q2 in q2_grid
+        ]
     if "integrability_version" in pineappl_grid.key_values():
         x_grid = evol_info.x1
         x_grid = np.append(x_grid, 1.0)

src/pineko/fonll.py

@@ -0,0 +1,190 @@
+"""Module to manage FONLL predictions."""
+
+import copy
+import json
+import tempfile
+from pathlib import Path
+
+import numpy as np
+import pineappl
+import yaml
+
+from . import configs
+
+
+class FONLLInfo:
+    """Class containing all the information for FONLL predictions."""
+
+    def __init__(
+        self, ffns3, ffn03, ffns4til, ffns4bar, ffn04, ffns5til, ffns5bar
+    ) -> None:
+        self.ffns3 = ffns3
+        self.ffn03 = ffn03
+        self.ffns4til = ffns4til
+        self.ffns4bar = ffns4bar
+        self.ffn04 = ffn04
+        self.ffns5til = ffns5til
+        self.ffns5bar = ffns5bar
+
+    @property
+    def fk_paths(self):
+        """Returns the list of the fk paths needed to produce FONLL predictions."""
+        always_needed = [self.ffns3, self.ffn03, self.ffns4til]
+        if self.ffns4bar:
+            paths = always_needed + [self.ffns4bar]
+        if self.ffn04 and self.ffns5til:
+            paths = paths + [self.ffn04, self.ffns5til]
+            # It does not make sense to have ffns5bar without ffns5til
+            if self.ffns5bar:
+                paths = paths + [self.self.ffns5bar]
+        return [Path(p) for p in paths]
+
+    @property
+    def fks(self):
+        """Returns the pineappl.Grid objects reading the paths provided by self.fk_paths."""
+        return [pineappl.grid.Grid.read(fk_path) for fk_path in self.fk_paths]
+
+    @property
+    def dataset_name(self):
+        """Return the name of the dataset given by the fktables name, if all the fktables have the same name."""
+        if len({p.name for p in self.fk_paths}) == 1:
+            return self.fk_paths[0].name
+        else:
+            raise ValueError("not all fktables have the same name")
+
+    @property
+    def theorycard_no_fns_pto(self):
+        """Return the common theory info between the different FONLL FKs."""
+        theorycards = [json.loads(fk.key_values()["theory"]) for fk in self.fks]
+        # Only these should differ
+        for card in theorycards:
+            del card["FNS"]
+            del card["PTO"]
+            del card["NfFF"]
+            del card["ID"]
+            del card["fonll-parts"]
+        if not all([theorycards[0] == card in theorycards[1:]]):
+            raise ValueError("theorycards not the same")
+        return theorycards[0]
+
+    @property
+    def Q2grid(self):
+        """Return the Q2grid of the fktables given by self.fks ."""
+        aa = json.loads(self.fks[0].raw.key_values()["runcard"])["observables"]
+        bb = list(aa.values())[
+            0
+        ]  # there is only a single obseravble because it's a dis fktable
+        cc = np.array([i["Q2"] for i in bb])
+        return cc
+
+
+# TODO: check if I want to damp the top mass corrections in 'mc' (i.e. 6)
+FK_TO_DAMP = {"mc": [1, 2, 4, 5], "mb": [4, 5]}
+FK_WITH_MINUS = [1, 4]  # asy terms should be subtracted, therefore the sign
+
+
+def produce_combined_fk(
+    ffns3, ffn03, ffns4til, ffns4bar, ffn04, ffns5til, ffns5bar, theoryid
+):
+    """Combine the FONLL FKs to produce one single final FONLL FK."""
+    fonll_info = FONLLInfo(ffns3, ffn03, ffns4til, ffns4bar, ffn04, ffns5til, ffns5bar)
+
+    theorycard_constituent_fks = fonll_info.theorycard_no_fns_pto
+    if theorycard_constituent_fks["DAMP"] == "False":
+        # then there is no damping, not even Heaviside only
+        combined_fk = fonll_info.fks[0]
+        for fk_path in fonll_info.fk_paths[1:]:
+            combined_fk.merge_from_file(fk_path)
+    else:
+        mc = theorycard_constituent_fks["mc"]
+        mb = theorycard_constituent_fks["mb"]
+        q2grid = fonll_info.Q2grid
+        step_function_charm = mc**2 < q2grid
+        step_function_bottom = mb**2 < q2grid
+        damping_factor_charm = (1 - q2grid / mc) ** theorycard_constituent_fks["DAMP"]
+        damping_factor_bottom = (1 - q2grid / mb) ** theorycard_constituent_fks["DAMP"]
+        damping_factor_charm *= step_function_charm
+        damping_factor_bottom *= step_function_bottom
+        dampings = {"mc": damping_factor_charm, "mb": damping_factor_bottom}
+        for i, fk in enumerate(fonll_info.fks):
+            sign = 1
+            if i in FK_WITH_MINUS:
+                sign = -1
+            for mass in FK_TO_DAMP:
+                if i in FK_TO_DAMP[mass]:
+                    fk.scale_by_bin(sign * dampings[mass])
+        # pineappl does not support operating with two grids in memory:
+        # https://github.com/NNPDF/pineappl/blob/8a672bef6d91b07a4edfdefbe4e30e4b1dd1f976/pineappl_py/src/grid.rs#L614-L617
+        with tempfile.TemporaryDirectory() as tmpdirname:
+            combined_fk = fonll_info.fks[0]
+            for i, fk in enumerate(fonll_info.fks[1:]):
+                tmpfile_path = Path(tmpdirname) / f"{i}.pineappl.lz4"
+                fk.write_lz4(tmpfile_path)
+                combined_fk.merge_from_file(tmpfile_path)
+
+    # update theorycard entries for the combined fktable by reading the yamldb of the original theory
+    input_theorycard_path = (
+        Path(configs.load(configs.detect())["paths"]["theory_cards"])
+        / f"{theoryid}.yaml"
+    )
+    with open(input_theorycard_path) as f:
+        final_theorycard = yaml.safe_load(f)
+    theorycard = json.loads(combined_fk.key_values()["theory"])
+    theorycard["FNS"] = final_theorycard["FNS"]
+    theorycard["PTO"] = final_theorycard["PTO"]
+    theorycard["NfFF"] = final_theorycard["NfFF"]
+    theorycard["ID"] = final_theorycard["ID"]
+    # Update the theorycard with the entries set above
+    combined_fk.set_key_value("theory", str(theorycard))
+
+    # save final "fonll" fktable
+    fk_folder = Path(configs.load(configs.detect())["paths"]["fktables"]) / str(
+        theoryid
+    )
+    fk_folder.mkdir(exist_ok=True)
+    output_path_fk = fk_folder / fonll_info.dataset_name
+    combined_fk.write_lz4(output_path_fk)
+
+
+def produce_fonll_recipe(fns):
+    """Produce the different theory cards according to which FONLL is asked for."""
+    fonll_recipe = []
+    if fns == "FONLL-A":
+        pto = [1 for _ in range(7)]
+    elif fns == "FONLL-C":
+        pto = [2 for _ in range(7)]
+    elif fns == "FONLL-B":
+        pto = [2, 2, 1, 2, 2, 1, 2]
+    FNS_list = [
+        "FONLL-FFNS",
+        "FONLL-FFN0",
+        "FONLL-FFNS",
+        "FONLL-FFNS",
+        "FONLL-FFN0",
+        "FONLL-FFNS",
+        "FONLL-FFNS",
+    ]
+    NfFF_list = [3, 3, 4, 4, 4, 5, 5, 5]
+    parts_list = ["full", "full", "massless", "massive", "full", "massless", "massive"]
+    for fns, nfff, po, part in zip(FNS_list, NfFF_list, pto, parts_list):
+        fonll_recipe.append(
+            {
+                "FNS": fns,
+                "NfFF": nfff,
+                "PTO": po,
+                "fonll-parts": part,
+            }
+        )
+    return fonll_recipe
+
+
+def produce_fonll_tcards(tcard, tcard_parent_path):
+    """Produce the seven fonll tcards from an original tcard and dump them in tcard_parent_path with names from '1000.yaml' to '1006.yaml'."""
+    theorycards = [copy.deepcopy(tcard) for _ in range(7)]
+    fonll_recipe = produce_fonll_recipe(tcard["FNS"])
+    for theorycard, recipe in zip(theorycards, fonll_recipe):
+        theorycard.update(recipe)
+    paths_list = [tcard_parent_path / f"100{num}.yaml" for num in range(7)]
+    for newtcard, path in zip(theorycards, paths_list):
+        with open(path, "w", encoding="UTF-8") as f:
+            yaml.safe_dump(newtcard, f)

src/pineko/parser.py

@@ -57,7 +57,6 @@ def get_yaml_information(yaml_file, grids_folder):
         List (of lists) with all the grids that will need to be loaded
     """
     yaml_content = _load_yaml(yaml_file)
-
     # Turn the operands and the members into paths (and check all of them exist)
     ret = []
     for operand in yaml_content["operands"]:

src/pineko/theory.py

@@ -388,6 +388,7 @@ def fk(self, name, grid_path, tcard, pdf):
             grid.lumi(),
         ):
             max_as += 1
+        # NB: This would not happen for simFONLL
         max_al = 0
         # check for sv
         if not np.isclose(xir, 1.0):