Implement custom energy-dependent cuts for IRFs and DL3 creations

magicctapipe/io/__init__.py

@@ -15,6 +15,7 @@
 from .io import (
     check_input_list,
     format_object,
+    get_custom_cuts,
     get_dl2_mean,
     get_stereo_events,
     get_stereo_events_old,
@@ -40,6 +41,7 @@
     "create_gti_hdu",
     "create_pointing_hdu",
     "format_object",
+    "get_custom_cuts",
     "get_dl2_mean",
     "get_stereo_events",
     "get_stereo_events_old",

magicctapipe/io/io.py

@@ -22,6 +22,7 @@
 from pyirf.binning import join_bin_lo_hi
 from pyirf.simulations import SimulatedEventsInfo
 from pyirf.utils import calculate_source_fov_offset, calculate_theta
+from scipy.interpolate import interp1d
 
 try:
     from importlib.resources import files
@@ -33,6 +34,7 @@
 __all__ = [
     "check_input_list",
     "format_object",
+    "get_custom_cuts",
     "get_dl2_mean",
     "get_stereo_events",
     "get_stereo_events_old",
@@ -249,6 +251,53 @@ def format_object(input_object):
     return string
 
 
+def get_custom_cuts(config, energy_bins):
+    """
+    Extract and interpolate custom, energy-dependent cuts for the IRFs creation.
+
+    Parameters
+    ----------
+    config : dict
+        Dictionary containing the custom cuts and the interpolation tragedy.
+        Required entries are a sub-dictionary 'custom_cuts' and a field
+        'interpolate_kind' used in scipy.interpolate.interp1d.
+    energy_bins : astropy.units.Quantity
+        Energy bin edges targeted in the IRFs
+
+    Returns
+    -------
+    astropy.table.QTable
+        Custom cut table in target 'energy_bins'
+    """
+    cut_table = QTable()
+    interpolate_kind = config["interpolate_kind"]
+    energy_bins_cuts = config["custom_cuts"]["energy_bins"] * u.Unit(
+        config["custom_cuts"]["energy_unit"]
+    )
+    energy_bins_cuts_center = 0.5 * (energy_bins_cuts[:-1] + energy_bins_cuts[1:]).to(
+        energy_bins.unit
+    )
+    cuts = config["custom_cuts"]["cut"] * u.Unit(config["custom_cuts"]["cut_unit"])
+
+    cut_table["low"] = energy_bins[:-1]
+    cut_table["high"] = energy_bins[1:]
+    cut_table["center"] = 0.5 * (cut_table["low"] + cut_table["high"])
+    f_cut = interp1d(
+        energy_bins_cuts_center,
+        cuts,
+        kind=interpolate_kind,
+        bounds_error=False,
+        fill_value=(cuts[0], cuts[-1]),
+        assume_sorted=True,
+    )
+
+    cut_table["cut"] = f_cut(cut_table["center"]) * u.Unit(
+        config["custom_cuts"]["cut_unit"]
+    )
+
+    return cut_table
+
+
 def get_stereo_events_old(
     event_data, quality_cuts=None, group_index=["obs_id", "event_id"]
 ):

magicctapipe/io/tests/test_io.py

@@ -1,13 +1,17 @@
 import glob
 
+import astropy.table
+import astropy.units as u
 import numpy as np
 import pandas as pd
 import pytest
+import yaml
 from ctapipe_io_lst import REFERENCE_LOCATION
 
 from magicctapipe.io.io import (
     check_input_list,
     format_object,
+    get_custom_cuts,
     get_dl2_mean,
     get_stereo_events,
     load_dl2_data_file,
@@ -152,6 +156,42 @@ def test_check_input_list(self):
                 }
             )
 
+    def test_get_custom_cuts(self, config):
+        with open(config, "rb") as f:
+            config = yaml.safe_load(f)
+        # Create energy bins
+        config_eng_bins = config["create_irf"]["energy_bins"]
+        eng_bins_start = u.Quantity(config_eng_bins["start"])
+        eng_bins_stop = u.Quantity(config_eng_bins["stop"])
+
+        energy_bins = u.TeV * np.geomspace(
+            start=eng_bins_start.to_value("TeV").round(3),
+            stop=eng_bins_stop.to_value("TeV").round(3),
+            num=config_eng_bins["n_edges"],
+        )
+
+        # Load the cuts configs. 'cut_type' will be ignored
+        config_gh_irf = config["create_irf"]["gammaness"]
+        config_theta_irf = config["create_irf"]["theta"]
+
+        # Get the custom cuts without interpolating on energy_bins
+        cut_table_gh = get_custom_cuts(config_gh_irf, energy_bins)
+        cut_table_theta = get_custom_cuts(config_theta_irf, energy_bins)
+        assert isinstance(cut_table_gh, astropy.table.QTable)
+        assert isinstance(cut_table_theta, astropy.table.QTable)
+        assert cut_table_gh["low"].unit == "TeV"
+        assert cut_table_gh["high"].unit == "TeV"
+        assert cut_table_gh["center"].unit == "TeV"
+        assert cut_table_gh["cut"].unit == u.dimensionless_unscaled
+        assert cut_table_gh["cut"][15] == 0.5
+        assert cut_table_theta["cut"].unit == u.deg
+        assert cut_table_theta["cut"][15] == 0.3 * u.deg
+
+        # Get the custom cuts with linear interpolating on energy_bins
+        config_gh_irf["interpolate_kind"] = "linear"
+        cut_table_gh = get_custom_cuts(config_gh_irf, energy_bins)
+        assert np.isclose(cut_table_gh["cut"][15].value, 0.5074, atol=0.0001)
+
     def test_telescope_combinations(self, config_gen, config_gen_4lst):
         """
         Simple check on telescope combinations

magicctapipe/resources/config.yaml

@@ -256,18 +256,34 @@ create_irf:
         n_edges: 21
 
     gammaness:
-        cut_type: "dynamic"  # select "global" or "dynamic"
+        cut_type: "dynamic"  # select "global" or "dynamic" or "custom"
         global_cut_value: 0.8  # used for the global cut
         efficiency: 0.9  # used for the dynamic cuts
         min_cut: 0.05  # used for the dynamic cuts
         max_cut: 0.85  # used for the dynamic cuts
+        custom_cuts: {  # used for custom cuts
+          energy_bins: [0.01, 0.1, 1, 10, 100],
+          energy_unit: "TeV",
+          cut: [0.4, 0.45, 0.50, 0.55],
+          cut_unit: "",
+        }
+        interpolate_kind: "nearest"  # used for custom cuts, option for scipy.interpolate.interp1d
+                                     # Standard values are "nearest", "linear" and "cubic"
 
     theta:
-        cut_type: "dynamic"  # select "global" or "dynamic"
+        cut_type: "dynamic"  # select "global" or "dynamic" or "custom"
         global_cut_value: "0.2 deg"  # used for the global cut
         efficiency: 0.75  # used for the dynamic cuts
         min_cut: "0.1 deg"  # used for the dynamic cuts
         max_cut: "0.25 deg"  # used for the dynamic cuts
+        custom_cuts: {  # used for custom cuts
+          energy_bins: [0.01, 0.1, 1, 10, 100],
+          energy_unit: "TeV",
+          cut: [0.4, 0.35, 0.30, 0.25],
+          cut_unit: "deg",
+        }
+        interpolate_kind: "nearest"  # used for custom cuts, option for scipy.interpolate.interp1d
+                                     # Standard values are "nearest", "linear" and "cubic"
 
 
 dl2_to_dl3:

magicctapipe/resources/test_config.yaml

@@ -226,19 +226,34 @@ create_irf:
         n_edges: 21
 
     gammaness:
-        cut_type: "dynamic"  # select "global" or "dynamic"
+        cut_type: "dynamic"  # select "global" or "dynamic" or "custom"
         global_cut_value: 0.8  # used for the global cut
         efficiency: 0.9  # used for the dynamic cuts
         min_cut: 0.05  # used for the dynamic cuts
         max_cut: 0.85  # used for the dynamic cuts
+        custom_cuts: {  # used for custom cuts
+          energy_bins: [0.01, 0.1, 1, 10, 100],
+          energy_unit: "TeV",
+          cut: [0.4, 0.45, 0.50, 0.55],
+          cut_unit: "",
+        }
+        interpolate_kind: "nearest"  # used for custom cuts, option for scipy.interpolate.interp1d
+                                     # Standard values are "nearest", "linear" and "cubic"
 
     theta:
-        cut_type: "dynamic"  # select "global" or "dynamic"
+        cut_type: "dynamic"  # select "global" or "dynamic" or "custom"
         global_cut_value: "0.2 deg"  # used for the global cut
         efficiency: 0.75  # used for the dynamic cuts
         min_cut: "0.1 deg"  # used for the dynamic cuts
         max_cut: "0.3 deg"  # used for the dynamic cuts
-
+        custom_cuts: {  # used for custom cuts
+          energy_bins: [0.01, 0.1, 1, 10, 100],
+          energy_unit: "TeV",
+          cut: [0.4, 0.35, 0.30, 0.25],
+          cut_unit: "deg",
+        }
+        interpolate_kind: "nearest"  # used for custom cuts, option for scipy.interpolate.interp1d
+                                     # Standard values are "nearest", "linear" and "cubic"
 
 dl2_to_dl3:
     interpolation_method: "nearest"  # select "nearest", "linear" or "cubic"