Enable nf=3 with `evolven3fit_new`.

Created by: scarlehoff

This PR enables the possibility of running evolven3fit_new with a fit in nf=3.

In order to make it more general now the operator receives a card where every time there is a threshold the Q_threshold is passed for both nf and nf+1 in order to avoid the hack of simply repeating one value.

The problem at the bottom threshold still remains.

Testing a fit evolved with apfel and eko for theory 200/400, they agree with a relative precision of 1e-3 for all points in x and for Q != mb (for nf=5 the agreement is a bit worse when you are close to mb). I've done the same test for theory 211 (perturbative charm, fit done in nf=3), there we again have a problem at the bottom threshold but as far as I can see the same problem doesn't appear at the charm threshold not sure why.

-> The problem at Q=mb for nf=5 is being investigated. Might be connected to this issue: https://github.com/NNPDF/eko/issues/173

-> The problem at Q=mb for nf=4 is due to eko which at the moment returns an object after evolution which is pdf(q), and pdf(q, nf=4) get's overwritten by pdf(q, nf=5). (same thing happens at nf=3 in theory).

This is still a draft because I would like to:

• Organize better the Q-grids in a separate file.
• Decide whether to set to 0.0 from n3fit the exportgrid for nf > nf0.

n3fit/src/evolven3fit_new/eko_utils.py

@@ -30,7 +30,6 @@ EVOLVEN3FIT_CONFIGS_DEFAULTS_EXA = {
 }
 
 NFREF_DEFAULT = 5
-NF0_DEFAULT = 4
 
 
 def construct_eko_cards(
@@ -58,31 +57,46 @@ 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
@@ -90,11 +104,11 @@ def construct_eko_cards(
 
     atlas = Atlas(
         matching_scales=MatchingScales(masses * thresholds_ratios),
-        origin=(theory["Q0"] ** 2, theory["nf0"]),
+        origin=(mu0**2, theory["nf0"]),
     )
     op_card.update(
         {
-            "mu0": theory["Q0"],
+            "mu0": mu0,
             "mugrid": [(float(np.sqrt(q2)), int(nf_default(q2, atlas))) for q2 in q2_grid],
         }
     )