Enable nf=3 with evolven3fit_new.

n3fit/src/evolven3fit_new/eko_utils.py

@@ -30,7 +30,6 @@
 }
 
 NFREF_DEFAULT = 5
-NF0_DEFAULT = 4
 
 
 def construct_eko_cards(
@@ -58,46 +57,74 @@ 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
     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"],
+        },
     )
     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":
@@ -119,7 +146,7 @@ def construct_eko_cards(
 
 
 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."""
+    """Split the evolgrid in blocks according to the number of flavors."""
     evolgrid_index_list = []
     evolgrid.sort()
     starting_nf = evolgrid[0][1]
@@ -131,7 +158,7 @@ def split_evolgrid(evolgrid):
     start_index = 0
     evolgrid_list = []
     for index in evolgrid_index_list:
-        evolgrid_list.append(evolgrid[start_index : index + 1])
+        evolgrid_list.append(evolgrid[start_index:index])
         start_index = index
     evolgrid_list.append(evolgrid[start_index:])
     return evolgrid_list

n3fit/src/evolven3fit_new/evolve.py

@@ -145,7 +145,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

n3fit/src/evolven3fit_new/utils.py

@@ -22,7 +22,6 @@
             3.8800751e00,
             4.3584516e00,
             4.9200000e00,
-            4.9200000e00,
             5.5493622e00,
             6.2897452e00,
             7.1650687e00,