correct sorting of level1 data

conda-recipe/meta.yaml

@@ -51,7 +51,7 @@ requirements:
         - sphinxcontrib-bibtex
         - curio >=1.0
         - pineappl >=0.5.8
-        - eko >=0.13.4,<0.14
+        - eko >=0.13.5,<0.14
         - banana-hep >=0.6.8
         - fiatlux
 

n3fit/runcards/examples/Basic_runcard_qed.yml

@@ -95,7 +95,7 @@ datacuts:
 
 ############################################################
 theory:
-  theoryid: 522     # database id
+  theoryid: 523     # database id
 
 ############################################################
 trvlseed: 1551864071
@@ -169,3 +169,4 @@ fiatlux:
   luxset: NNPDF40_nnlo_as_01180
   additional_errors: true # should be set to true only for the last iteration
   luxseed: 1234567890
+  eps_base: 1e-2

n3fit/src/evolven3fit_new/eko_utils.py

@@ -22,15 +22,14 @@
 }
 
 EVOLVEN3FIT_CONFIGS_DEFAULTS_EXA = {
-    "ev_op_iterations": 10,
+    "ev_op_iterations": 30,
     "ev_op_max_order": (1, 0),
     "evolution_method": "iterate-exact",
     "inversion_method": "exact",
     "n_integration_cores": 1,
 }
 
 NFREF_DEFAULT = 5
-NF0_DEFAULT = 4
 
 
 def construct_eko_cards(
@@ -58,46 +57,77 @@ def construct_eko_cards(
     theory = Loader().check_theoryID(theoryID).get_description()
     theory.pop("FNS")
     theory.update(theory_card_dict)
-    if "nfref" not in theory:
-        theory["nfref"] = NFREF_DEFAULT
-    if "nf0" not in theory:
-        theory["nf0"] = NF0_DEFAULT
-
+
     # Prepare the thresholds according to MaxNfPdf
     thresholds = {"c": theory["kcThr"], "b": theory["kbThr"], "t": theory["ktThr"]}
     if theory["MaxNfPdf"] < 5:
         thresholds["b"] = np.inf
     if theory["MaxNfPdf"] < 6:
         thresholds["t"] = np.inf
-
+
+    if "nfref" not in theory:
+        theory["nfref"] = NFREF_DEFAULT
+
+    # Set nf_0 according to the fitting scale unless set explicitly
+    mu0 = theory["Q0"]
+    if "nf0" not in theory:
+        if mu0 < theory["mc"] * thresholds["c"]:
+            theory["nf0"] = 3
+        elif mu0 < theory["mb"] * thresholds["b"]:
+            theory["nf0"] = 4
+        elif mu0 < theory["mt"] * thresholds["t"]:
+            theory["nf0"] = 5
+        else:
+            theory["nf0"] = 6
+
     # Setting the thresholds in the theory card to inf if necessary
-    theory.update({"kbThr":thresholds["b"], "ktThr":thresholds["t"] })
-    
+    theory.update({"kbThr": thresholds["b"], "ktThr": thresholds["t"]})
+
     # The Legacy function is able to construct a theory card for eko starting from an NNPDF theory
     legacy_class = runcards.Legacy(theory, {})
     theory_card = legacy_class.new_theory
 
-    # construct operator card
+    # Generate the q2grid, if q_fin and q_points are None, use `nf0` to select a default
     q2_grid = utils.generate_q2grid(
-        theory["Q0"],
+        mu0,
         q_fin,
         q_points,
-        {theory["mb"]: thresholds["b"], theory["mt"]: thresholds["t"]},
+        {
+            theory["mc"]: thresholds["c"],
+            theory["mb"]: thresholds["b"],
+            theory["mt"]: thresholds["t"],
+        },
+        theory["nf0"],
     )
+
+    # construct operator card
     op_card = default_op_card
     masses = np.array([theory["mc"], theory["mb"], theory["mt"]]) ** 2
     thresholds_ratios = np.array([thresholds["c"], thresholds["b"], thresholds["t"]]) ** 2
 
     atlas = Atlas(
         matching_scales=MatchingScales(masses * thresholds_ratios),
-        origin=(theory["Q0"] ** 2, theory["nf0"]),
-    )
-    op_card.update(
-        {
-            "mu0": theory["Q0"],
-            "mugrid": [(float(np.sqrt(q2)), int(nf_default(q2, atlas))) for q2 in q2_grid],
-        }
+        origin=(mu0**2, theory["nf0"]),
     )
+
+    # Create the eko operator q2grid
+    # This is a grid which contains information on (q, nf)
+    # in VFNS values at the matching scales need to be doubled so that they are considered in both sides
+
+    ep = 1e-4
+    mugrid = []
+    for q2 in q2_grid:
+        q = float(np.sqrt(q2))
+        if nf_default(q2 + ep, atlas) != nf_default(q2 - ep, atlas):
+            nf_l = int(nf_default(q2 - ep, atlas))
+            nf_u = int(nf_default(q2 + ep, atlas))
+            mugrid.append((q, nf_l))
+            mugrid.append((q, nf_u))
+        else:
+            mugrid.append((q, int(nf_default(q2, atlas))))
+
+    op_card.update({"mu0": theory["Q0"], "mugrid": mugrid})
+
     op_card["xgrid"] = x_grid
     # Specific defaults for evolven3fit evolution
     if theory["ModEv"] == "TRN":
@@ -116,22 +146,3 @@ def construct_eko_cards(
 
     op_card = runcards.OperatorCard.from_dict(op_card)
     return theory_card, op_card
-
-
-def split_evolgrid(evolgrid):
-    """Split the evolgrid in blocks according to the number of flavors and repeating the last entry of one block in the first entry of the next."""
-    evolgrid_index_list = []
-    evolgrid.sort()
-    starting_nf = evolgrid[0][1]
-    for evo_point in evolgrid:
-        current_nf = evo_point[1]
-        if current_nf != starting_nf:
-            evolgrid_index_list.append(evolgrid.index(evo_point))
-            starting_nf = current_nf
-    start_index = 0
-    evolgrid_list = []
-    for index in evolgrid_index_list:
-        evolgrid_list.append(evolgrid[start_index : index + 1])
-        start_index = index
-    evolgrid_list.append(evolgrid[start_index:])
-    return evolgrid_list

n3fit/src/evolven3fit_new/evolve.py

@@ -1,3 +1,4 @@
+from collections import defaultdict
 import logging
 import pathlib
 import sys
@@ -145,7 +146,7 @@ def load_fit(usr_path):
     """
     nnfitpath = usr_path / "nnfit"
     pdf_dict = {}
-    for yaml_file in nnfitpath.glob("replica_*/*.exportgrid"):
+    for yaml_file in nnfitpath.glob(f"replica_*/{usr_path.name}.exportgrid"):
         data = yaml.safe_load(yaml_file.read_text(encoding="UTF-8"))
         pdf_dict[yaml_file.parent.stem] = data
     return pdf_dict
@@ -178,19 +179,27 @@ def evolve_exportgrid(exportgrid, eko, x_grid, qed):
     # generate block to dump
     targetgrid = eko.bases.targetgrid.tolist()
 
-    def ev_pdf(pid, x, Q2):
-        return x * evolved_pdf[Q2]["pdfs"][pid][targetgrid.index(x)]
+    # Finally separate by nf block (and order per nf/q)
+    by_nf = defaultdict(list)
+    for q, nf in sorted(eko.evolgrid, key=lambda ep: ep[1]):
+        by_nf[nf].append(q)
+    q2block_per_nf = {nf: sorted(qs) for nf, qs in by_nf.items()}
 
-    evolgrid_list = eko_utils.split_evolgrid(eko.evolgrid)
     blocks = []
-    for evgrid in evolgrid_list:
+    for nf, q2grid in q2block_per_nf.items():
+
+        def pdf_xq2(pid, x, Q2):
+            x_idx = targetgrid.index(x)
+            return x * evolved_pdf[(Q2, nf)]["pdfs"][pid][x_idx]
+
         block = genpdf.generate_block(
-            ev_pdf,
+            pdf_xq2,
             xgrid=targetgrid,
-            evolgrid=evgrid,
+            sorted_q2grid=q2grid,
             pids=basis_rotation.flavor_basis_pids,
         )
         blocks.append(block)
+
     return blocks
 
 

n3fit/src/evolven3fit_new/utils.py

@@ -7,7 +7,7 @@
 from reportengine.compat import yaml
 from validphys.pdfbases import PIDS_DICT
 
-DEFAULT_Q2GRID = (
+Q2GRID_Nf04 = (
     np.array(
         [
             1.6500000e00,
@@ -22,7 +22,6 @@
             3.8800751e00,
             4.3584516e00,
             4.9200000e00,
-            4.9200000e00,
             5.5493622e00,
             6.2897452e00,
             7.1650687e00,
@@ -65,6 +64,62 @@
     ** 2
 )
 
+Q2GRID_Nf03 = (
+    np.array(
+        [
+            1.0000000e00,
+            1.0768843e00,
+            1.1642787e00,
+            1.2640247e00,
+            1.3783565e00,
+            1.5100000e00,
+            1.6573843e00,
+            1.8279487e00,
+            2.0263188e00,
+            2.2582323e00,
+            2.5308507e00,
+            2.8531703e00,
+            3.2365690e00,
+            3.6955380e00,
+            4.2486693e00,
+            4.9200000e00,
+            5.6571821e00,
+            6.5475141e00,
+            7.6300446e00,
+            8.9555329e00,
+            1.0590474e01,
+            1.2622686e01,
+            1.5169120e01,
+            1.8386905e01,
+            2.2489085e01,
+            2.7767274e01,
+            3.4624624e01,
+            4.3624282e01,
+            5.5561424e01,
+            7.1571582e01,
+            9.3295496e01,
+            1.2313315e02,
+            1.6464038e02,
+            2.2315640e02,
+            3.0681103e02,
+            4.2816505e02,
+            6.0692308e02,
+            8.7449251e02,
+            1.2817733e03,
+            1.9127020e03,
+            2.9082314e03,
+            4.5095982e03,
+            7.1379509e03,
+            1.1543948e04,
+            1.9094934e04,
+            3.2338760e04,
+            5.6137084e04,
+            1.0000000e05,
+        ]
+    )
+    ** 2
+)
+
 
 class LhapdfLike:
     """
@@ -123,15 +178,22 @@ def get_theoryID_from_runcard(usr_path):
     return my_runcard["theory"]["theoryid"]
 
 
-def generate_q2grid(Q0, Qfin, Q_points, match_dict):
+def generate_q2grid(Q0, Qfin, Q_points, match_dict, nf0=None):
     """Generate the q2grid used in the final evolved pdfs or use the default grid if Qfin or Q_points is
     not provided.
 
     match_dict contains the couples (mass : factor) where factor is the number to be multiplied to mass
     in order to obtain the relative matching scale.
     """
     if Qfin is None and Q_points is None:
-        return DEFAULT_Q2GRID
+        if nf0 == 4:
+            return Q2GRID_Nf04
+        elif nf0 == 3:
+            return Q2GRID_Nf03
+        elif nf0 is None:
+            raise ValueError("In order to use a default grid, a value of nf0 must be provided")
+        else:
+            raise NotImplementedError(f"No default grid in Q available for {nf0=}")
     elif Qfin is None or Q_points is None:
         raise ValueError("q_fin and q_points must be specified either both or none of them")
     else:
@@ -148,7 +210,9 @@ def generate_q2grid(Q0, Qfin, Q_points, match_dict):
                 frac_of_point = np.log(match_scale / Q_ini) / np.log(Qfin / Q0)
                 num_points = int(Q_points * frac_of_point)
                 num_points_list.append(num_points)
-                grids.append(np.geomspace(Q_ini**2, match_scale**2, num=num_points))
+                grids.append(
+                    np.geomspace(Q_ini**2, match_scale**2, num=num_points, endpoint=False)
+                )
                 Q_ini = match_scale
         num_points = Q_points - sum(num_points_list)
         grids.append(np.geomspace(Q_ini**2, Qfin**2, num=num_points))