some cleaning in scale variations

benchmarks/bench_cli.py

@@ -113,18 +113,15 @@ def benchmark_scaffold_cli(test_empty_proj):
     assert "Success: All the folders are correctly configured" in res.output
 
 
-def benchmark_gen_sv_cli(test_files, tmp_path, test_pdf, lhapdf_path):
+def benchmark_gen_sv_cli(test_files, tmp_path):
     runner = CliRunner()
-    pdf_name = "NNPDF40_nlo_as_01180"
     max_as = "2"
     nf = "5"
     name_grid = "ATLAS_TTB_8TEV_LJ_TRAP_norensv_fixed.pineappl.lz4"
     grid_path = test_files / "data" / "grids" / "400" / name_grid
     new_grid_path = tmp_path / name_grid
     target_path = tmp_path
     shutil.copy(grid_path, new_grid_path)
-    with lhapdf_path(test_pdf):
-        pdf = lhapdf.mkPDF(pdf_name)
     res = runner.invoke(
         command,
         ["ren_sv_grid", str(new_grid_path), str(target_path), max_as, nf, "False"],

src/pineko/kfactor.py

@@ -9,6 +9,7 @@
 from pineappl import import_only_subgrid
 
 from . import scale_variations
+from . scale_variations import orders_as_tuple
 
 DEFAULT_PDF_SET = "NNPDF40_nnlo_as_01180"
 
@@ -156,7 +157,7 @@ def construct_new_order(grid, order, order_to_update, central_kfactor, alphas):
         alphas
     """
     # extract the relevant order to rescale from the grid for each lumi and bin
-    grid_orders = [order.as_tuple() for order in grid.orders()]
+    grid_orders = orders_as_tuple(grid)
 
     new_grid = scale_variations.initialize_new_grid(grid, order_to_update)
     orginal_order_index = grid_orders.index(order)
@@ -215,14 +216,15 @@ def do_it(
     order_exists: bool
         True if the order to update is already present
     """
-    grid_orders = [order.as_tuple() for order in grid.orders()]
+    grid_orders = orders_as_tuple(grid)
 
     # remove not necessary orders
-    # NOTE: eventual QED corrections are not supported
+    # NOTE: eventual QED kfactors are not supported
     order_mask = pineappl.grid.Order.create_mask(grid.orders(), pto_to_update, 0, True)
     grid_orders_filtered = list(np.array(grid_orders)[order_mask])
     grid_orders_filtered.sort(key=scale_variations.qcd)
     min_as = grid_orders_filtered[0][0]
+    # TODO: this is always going to be 0, given the mask above ...
     min_al = grid_orders_filtered[0][1]
 
     # the actual alpha_s order to update
@@ -243,7 +245,7 @@ def do_it(
     max_as = grid_orders_filtered[-1][0]
     orders_list = [(de, min_al, 0, 0) for de in range(min_as, max_as + 1)]
     # create an empty grid and add the rescaled order
-    order_to_update = (order_to_update, grid_orders_filtered[0][1], 0, 0)
+    order_to_update = (order_to_update, min_al, 0, 0)
     new_order_grid = None
     for i, as_order in enumerate(orders_list):
         order_grid = construct_new_order(

src/pineko/scale_variations.py

@@ -30,6 +30,8 @@ def qcd(order: OrderTuple) -> int:
     """Extract the QCD order from an OrderTuple."""
     return order[0]
 
+def orders_as_tuple(grid: pineappl.grid.Grid) -> list[OrderTuple]:
+    return [order.as_tuple() for order in grid.orders()]
 
 def ren_sv_coeffs(m, max_as, logpart, which_part, nf):
     """Ren_sv coefficient for the requested part.
@@ -114,7 +116,7 @@ def create_svonly(grid, order, new_order, scalefactor):
     """Create a grid containing only the renormalization scale variations at a given order for a grid."""
     new_grid = initialize_new_grid(grid, new_order)
     # extract the relevant order to rescale from the grid for each lumi and bin
-    grid_orders = [order.as_tuple() for order in grid.orders()]
+    grid_orders = orders_as_tuple(grid)
     order_index = grid_orders.index(order)
     for lumi_index in range(len(new_grid.lumi())):
         for bin_index in range(grid.bins()):
@@ -137,7 +139,7 @@ def create_svonly(grid, order, new_order, scalefactor):
 def create_grids(gridpath, max_as, nf):
     """Create all the necessary scale variations grids for a certain starting grid."""
     grid = pineappl.grid.Grid.read(gridpath)
-    grid_orders = [orde.as_tuple() for orde in grid.orders()]
+    grid_orders =  orders_as_tuple(grid)
     order_mask = pineappl.grid.Order.create_mask(grid.orders(), max_as, 0, True)
     grid_orders_filtered = list(np.array(grid_orders)[order_mask])
     grid_orders_filtered.sort(key=qcd)
@@ -159,47 +161,7 @@ def create_grids(gridpath, max_as, nf):
             )
         grid_list[to_construct_order] = list_grid_order
 
-    return grid_list, nec_orders
-
-
-def write_grids(gridpath, grid_list):
-    """Write the single grids."""
-    base_name = gridpath.stem.split(".pineappl")[0]
-    final_part = ".pineappl.lz4"
-    grid_paths = []
-    for order in grid_list:
-        # For each scale variation order, if more than one grid contributes, merge them all together in a single one
-        if len(grid_list[order]) > 1:
-            for grid in grid_list[order][1:]:
-                tmp_path = gridpath.parent / ("tmp" + final_part)
-                grid.write_lz4(tmp_path)
-                grid_list[order][0].merge_from_file(tmp_path)
-                tmp_path.unlink()
-        new_grid_path = gridpath.parent / (
-            base_name + "_" + str(order[2]) + final_part
-        )  # order[2] is the ren_sv order
-        grid_paths.append(new_grid_path)
-        grid_list[order][0].write_lz4(new_grid_path)
-    return grid_paths
-
-
-def merge_grids(
-    gridpath, grid_list_path, target_path=None, nec_orders=None, order_exists=False
-):
-    """Merge the single grids in the original."""
-    grid = pineappl.grid.Grid.read(gridpath)
-    if target_path is None:
-        target_path = gridpath.parent / gridpath.name
-    else:
-        target_path = target_path / gridpath.name
-    if order_exists:
-        grid = construct_and_dump_order_exists_grid(
-            grid, list(nec_orders)[0] if nec_orders is not None else []
-        )
-    for grid_path in grid_list_path:
-        grid.merge_from_file(grid_path)
-        grid_path.unlink()
-    grid.write_lz4(target_path)
+    return grid_list
 
 
 def construct_and_dump_order_exists_grid(ori_grid, to_construct_order):
@@ -212,13 +174,11 @@ def construct_and_dump_order_exists_grid(ori_grid, to_construct_order):
     to_construct_order:
         order to delete
     """
-
-    # TODO: maybe we can make this function simpler ??
-
+    # TODO: can we make this function simpler ??
     bin_limits = [float(bin) for bin in range(ori_grid.bins() + 1)]
     lumi_grid = [pineappl.lumi.LumiEntry(mylum) for mylum in ori_grid.lumi()]
     subgrid_params = pineappl.subgrid.SubgridParams()
-    ori_grid_orders = [order.as_tuple() for order in ori_grid.orders()]
+    ori_grid_orders = orders_as_tuple(ori_grid)
     new_orders = [
         pineappl.grid.Order(*ord)
         for ord in ori_grid_orders
@@ -268,9 +228,8 @@ def compute_ren_sv_grid(
 ):
     """Generate renormalization scale variation terms for the given grid, according to the max_as."""
     # First let's check if the ren_sv are already there
-    checkres, max_as_effective = check.contains_sv(
-        pineappl.grid.Grid.read(grid_path), max_as, 0, check.Scale.REN
-    )
+    grid = pineappl.grid.Grid.read(grid_path)
+    checkres, max_as_effective = check.contains_sv(grid, max_as, 0, check.Scale.REN)
     # Usual different convention with max_as
     if max_as_effective == max_as and (checkres is not check.AvailableAtMax.CENTRAL):
         if not order_exists:
@@ -279,18 +238,21 @@ def compute_ren_sv_grid(
         return ReturnState.ORDER_EXISTS_FAILURE
     if max_as_effective < max_as and checkres is check.AvailableAtMax.SCVAR:
         return ReturnState.MISSING_CENTRAL
+    # Create all the necessary grids
     # With respect to the usual convention here max_as is max_as-1
-    max_as -= 1
-    # Creating all the necessary grids
-    grid_list, nec_orders = create_grids(grid_path, max_as, nf)
-    # Writing the sv grids
-    sv_grids_paths = write_grids(gridpath=grid_path, grid_list=grid_list)
-    # Merging all together
-    merge_grids(
-        gridpath=grid_path,
-        grid_list_path=sv_grids_paths,
-        target_path=target_path,
-        nec_orders=nec_orders,
-        order_exists=order_exists,
-    )
+    grid_list = create_grids(grid_path, max_as - 1, nf)
+
+    # merge them
+    for sv_order, sv_grids in grid_list.items():
+        # if the new order is there, clean the old one.
+        if sv_order in orders_as_tuple(grid):
+            grid = construct_and_dump_order_exists_grid(
+                grid, list(grid_list)[0]
+            )
+        for sv_grid in sv_grids:
+            grid.merge(sv_grid)
+    # save
+    if target_path is None:
+        target_path = grid_path.parent
+    grid.write_lz4(target_path / grid_path.name)
     return ReturnState.SUCCESS