add a default q0 grid for nf0=3 and associated tests

n3fit/src/evolven3fit_new/eko_utils.py

@@ -97,6 +97,7 @@ def construct_eko_cards(
             theory["mb"]: thresholds["b"],
             theory["mt"]: thresholds["t"],
         },
+        theory["nf0"],
     )
     op_card = default_op_card
     masses = np.array([theory["mc"], theory["mb"], theory["mt"]]) ** 2

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,
@@ -64,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:
     """
@@ -122,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:
@@ -147,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))

n3fit/src/n3fit/tests/test_evolven3fit.py

@@ -23,11 +23,6 @@ def assert_sorted(arr, title):
         raise ValueError(f"The values of {title} are not sorted!")
 
 
-def check_consecutive_members(grid, value):
-    """Check if the first occurrence of value in grid is followed by value again"""
-    return np.allclose(grid[list(grid).index(value) + 1], value)
-
-
 def check_lhapdf_info(info_path):
     """Check the LHAPDF info file is correct"""
     info = yaml.load(info_path.open("r", encoding="utf-8"))
@@ -77,24 +72,40 @@ def check_lhapdf_dat(dat_path, info):
     # Use allclose here to avoid failing because of a different in the 7th place
     np.testing.assert_allclose(q[-1], info["QMax"])
 
+
+def test_generate_q2grid():
+    """Tests the creation of the default grids for different values of nf
+    and whether the matched grid is generating points in the desired locations
+    """
+    # nf 3, q0 = 1.0
+    grid = utils.generate_q2grid(None, None, None, {}, 3)
+    assert grid[0] == 1.0**2
+    # nf 4, q0 = 1.65
+    grid = utils.generate_q2grid(None, None, None, {}, 4)
+    assert grid[0] == 1.65**2
+
+    for nf in [1, 2, 5, 6]:
+        with pytest.raises(NotImplementedError):
+            grid = utils.generate_q2grid(None, None, None, {}, nf)
+
+    with pytest.raises(ValueError):
+        grid = utils.generate_q2grid(None, None, None, {})
 
-def test_utils():
-    # Testing the default grid
-    grid = utils.generate_q2grid(1.65, None, None, {})
-    assert_allclose(1.65**2, grid[0])
-    assert len(grid) == 50
-    # We expect the bottom mass to be repeated twice because it is intended once in 4 flavors and once in 5 flavors.
-    assert check_consecutive_members(grid, 4.92**2)
-    # Testing if the points of the matching are correctly repeated twice
     matched_grid = utils.generate_q2grid(1.65, 1.0e5, 100, {4.92: 2.0, 100: 1.0})
-    assert len(matched_grid) == 100
+    t1 = 4.92 * 2.0
+    t2 = 100.0 * 1.0
+
+    assert_allclose((1.65) ** 2, matched_grid[0])
     assert_allclose((1.0e5) ** 2, matched_grid[-1])
-    assert check_consecutive_members(matched_grid, (4.92 * 2.0) ** 2)
-    assert check_consecutive_members(matched_grid, (100.0 * 1.0) ** 2)
+    assert t1**2 in matched_grid
+    assert t2**2 in matched_grid
+
+
+def test_utils():
     # Testing the fake LHAPDF class
     q20 = 1.65**2
     x_grid = np.geomspace(1.0e-7, 1.0, 30)
-    fake_grids = [[x * (1.0 - x) for x in x_grid] for pid in PIDS_DICT.keys()]
+    fake_grids = [[x * (1.0 - x) for x in x_grid] for _ in PIDS_DICT.keys()]
     pdf_grid = dict([(pid, v) for pid, v in zip(range(len(PIDS_DICT)), fake_grids)])
     my_PDF = utils.LhapdfLike(pdf_grid, q20, x_grid)
     assert my_PDF.hasFlavor(6)
@@ -133,10 +144,15 @@ def test_eko_utils(tmp_path):
         t_card_dict["order"][0] == pto + 1
     )  # This is due to a different convention in eko orders due to QED
     assert_allclose(op_card_dict["xgrid"], x_grid)
-    assert_allclose(op_card_dict["mugrid"][0], (1.65, 4))
+    # In theory 162 the charm threshold is at 1.51
+    # and we should find two entries, one for nf=3 and another one for nf=4
+    assert_allclose(op_card_dict["mugrid"][0], (1.51, 3))
+    assert_allclose(op_card_dict["mugrid"][1], (1.51, 4))
+    # Then (with the number of points we chosen it will happen in position 2,3
+    # we will find the bottom threshold at two different nf
+    assert_allclose(op_card_dict["mugrid"][2], (4.92, 4))
+    assert_allclose(op_card_dict["mugrid"][3], (4.92, 5))
     assert_allclose(op_card_dict["mugrid"][-1], (q_fin, 5))
-    # In this case there are not enough points to have twice the bottom matching scale
-    assert_allclose(op_card_dict["mugrid"][1], (4.92, 5))
     # Testing computation of eko
     save_path = tmp_path / "ekotest.tar"
     runner.solve(t_card, op_card, save_path)