Allow using custom energy dependent cuts (e.g. based on sensitivity optimization) for IRFs and DL3

docs/examples/irf_dl3_tool_config_edepcuts.json

@@ -0,0 +1,56 @@
+{
+  "EventSelector": {
+    "filters": {
+      "intensity": [50, Infinity],
+      "width": [0, Infinity],
+      "length": [0, Infinity],
+      "r": [0, 1],
+      "wl": [0.01, 1],
+      "leakage_intensity_width_2": [0, 1],
+      "event_type": [32, 32]
+    }
+  },
+  "IRFFITSWriter": {
+    "cut_strategy_gh": "edep_custom",
+    "cut_strategy_theta": "edep_custom",
+    "cut_strategy_alpha": "edep_custom",
+    "interpolate_custom": "nearest"
+  },
+  "DL3Cuts": {
+    "min_event_p_en_bin": 100,
+    "energy_dependent_gh_cut": {
+      "energy_bins": [0.01, 0.1, 1, 10, 100],
+      "cut": [0.4, 0.45, 0.50, 0.55]
+    },
+    "energy_dependent_theta_cut": {
+      "energy_bins": [0.01, 0.1, 1, 10, 100],
+      "cut": [0.37, 0.35, 0.3, 0.25]
+    },
+    "energy_dependent_alpha_cut": {
+      "energy_bins": [0.01, 0.1, 1, 10, 100],
+      "cut": [4, 4.5, 5, 5.5]
+    },
+    "allowed_tels": [1]
+  },
+  "DataBinning": {
+    "true_energy_min": 0.005,
+    "true_energy_max": 500,
+    "true_energy_n_bins": 25,
+    "scale_true_energy": 1.0,
+    "reco_energy_min": 0.005,
+    "reco_energy_max": 500,
+    "reco_energy_n_bins": 25,
+    "energy_migration_min": 0.2,
+    "energy_migration_max": 5,
+    "energy_migration_n_bins": 30,
+    "fov_offset_min": 0.1,
+    "fov_offset_max": 1.1,
+    "fov_offset_n_edges": 9,
+    "bkg_fov_offset_min": 0,
+    "bkg_fov_offset_max": 10,
+    "bkg_fov_offset_n_edges": 21,
+    "source_offset_min": 0,
+    "source_offset_max": 1,
+    "source_offset_n_edges": 101
+  }
+}

lstchain/high_level/tests/test_interpolate.py

@@ -61,8 +61,8 @@ def test_interp_irf(
         simulated_dl2_file,
         "--output-irf-file",
         irf_file_en_1,
-        "--energy-dependent-gh",
-        "--energy-dependent-theta",
+        "--cut-strategy-gh=edep_efficiency",
+        "--cut-strategy-theta=edep_efficiency",
         "--point-like",
         "--gh-efficiency=0.8",
         "--theta-containment=0.8",
@@ -92,8 +92,8 @@ def test_interp_irf(
         simulated_srcdep_dl2_file,
         "--output-irf-file",
         irf_file_en_srcdep_1,
-        "--energy-dependent-gh",
-        "--energy-dependent-alpha",
+        "--cut-strategy-gh=edep_efficiency",
+        "--cut-strategy-alpha=edep_efficiency",
         "--point-like",
         "--source-dep",
         "--gh-efficiency=0.8",
@@ -393,8 +393,8 @@ def test_compare_irfs(
         simulated_dl2_file,
         "--output-irf-file",
         irf_file_en_1_cut2,
-        "--energy-dependent-gh",
-        "--energy-dependent-theta",
+        "--cut-strategy-gh=edep_efficiency",
+        "--cut-strategy-theta=edep_efficiency",
         "--point-like",
         "--gh-efficiency=0.7",
         "--theta-containment=0.7",
@@ -432,8 +432,8 @@ def test_compare_irfs(
         simulated_srcdep_dl2_file,
         "--output-irf-file",
         srcdep_irf_file_en_1_cut2,
-        "--energy-dependent-gh",
-        "--energy-dependent-alpha",
+        "--cut-strategy-gh=edep_efficiency",
+        "--cut-strategy-alpha=edep_efficiency",
         "--point-like",
         "--source-dep",
         "--gh-efficiency=0.7",

lstchain/io/event_selection.py

@@ -1,5 +1,6 @@
 import numpy as np
 import operator
+from astropy.table import QTable
 import astropy.units as u
 
 from ctapipe.core import Component
@@ -10,6 +11,8 @@
 
 from pyirf.cuts import calculate_percentile_cut, evaluate_binned_cut
 
+from scipy.interpolate import interp1d
+
 
 __all__ = ["EventSelector", "DL3Cuts", "DataBinning"]
 
@@ -67,6 +70,11 @@ class DL3Cuts(Component):
         default_value=0.95,
     ).tag(config=True)
 
+    energy_dependent_gh_cut = Dict(
+        help="Energy dependent selection cuts for gh_score (gammaness)",
+        default_value={},
+    ).tag(config=True)
+
     min_gh_cut = Float(
         help="Minimum gh_score (gammaness) cut in an energy bin",
         default_value=0.1,
@@ -105,6 +113,11 @@ class DL3Cuts(Component):
         default_value=0.2,
     ).tag(config=True)
 
+    energy_dependent_theta_cut = Dict(
+        help="Energy dependent selection cuts (deg) for theta",
+        default_value={},
+    ).tag(config=True)
+
     min_alpha_cut = Float(
         help="Minimum alpha cut (deg) in an energy bin",
         default_value=1,
@@ -133,6 +146,11 @@ class DL3Cuts(Component):
         default_value=20,
     ).tag(config=True)
 
+    energy_dependent_alpha_cut = Dict(
+        help="Energy dependent selection cuts (deg) for alpha",
+        default_value={},
+    ).tag(config=True)
+
     allowed_tels = List(
         help="List of allowed LST telescope ids",
         trait=Int(),
@@ -333,6 +351,47 @@ def apply_energy_dependent_alpha_cuts(self, data, alpha_cuts):
         )
         return data[data["selected_alpha"]]
 
+
+    def from_dict(self, parameter, energy_bins, interpolate_kind='nearest'):
+        """
+        Convert the cut dictionary to a QTable for the requested parameter.
+
+        Parameters
+        ----------
+        parameter: string
+            'gh', 'theta' or 'alpha'
+        energy_bins: `astropy.units.quantity.Quantity`
+            energy bins used for the IRFs
+        interpolate_kind: string
+            Interpolation strategy for `scipy.interpolate.interp1d`
+
+        Returns
+        -------
+        cut_table: QTable
+            Energy dependent cuts for the requested parameter.
+
+        """
+        cut_table = QTable()
+
+        cuts_dict = {"gh": self.energy_dependent_gh_cut,
+                     "theta": self.energy_dependent_theta_cut,
+                     "alpha":self.energy_dependent_alpha_cut}[parameter]
+
+        unit = {"gh": 1,
+                "theta": u.deg,
+                "alpha": u.deg}[parameter]
+
+        input_ebins_center = (cuts_dict["energy_bins"][:-1]* u.TeV + cuts_dict["energy_bins"][1:] * u.TeV) * 0.5
+
+        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(input_ebins_center, cuts_dict["cut"], kind=interpolate_kind, bounds_error=False,
+                       fill_value=(cuts_dict["cut"][0], cuts_dict["cut"][-1]), assume_sorted=True)
+        cut_table["cut"] = f_cut(cut_table["center"]) * unit
+
+        return cut_table
+
     def allowed_tels_filter(self, data):
         """
         Applying a filter on telescopes used for observation.

lstchain/tools/lstchain_create_dl3_file.py

@@ -361,10 +361,16 @@ def apply_srcindep_gh_cut(self):
             self.data = self.cuts.apply_energy_dependent_gh_cuts(
                 self.data, self.energy_dependent_gh_cuts
             )
-            self.log.info(
-                "Using gamma efficiency of "
-                f"{self.energy_dependent_gh_cuts.meta['GH_EFF']}"
-            )
+            if 'GH_EFF' in self.energy_dependent_gh_cuts.meta.keys():
+                self.log.info(
+                    "Using gamma efficiency of "
+                    f"{self.energy_dependent_gh_cuts.meta['GH_EFF']}"
+                )
+            else:
+                self.log.info(
+                    "Using energy dependent gammaness cuts : "
+                    f"{self.energy_dependent_gh_cuts.meta['GH_CUTS']}"
+                )
         else:
             self.cuts.global_gh_cut = self.irf_final_hdu[1].header["GH_CUT"]
             self.data = self.cuts.apply_global_gh_cut(self.data)
@@ -391,10 +397,16 @@ def apply_srcdep_gh_alpha_cut(self):
                 data_temp = self.cuts.apply_energy_dependent_gh_cuts(
                     data_temp, self.energy_dependent_gh_cuts
                 )
-                self.log.info(
-                    "Using gamma efficiency of "
-                    f"{self.energy_dependent_gh_cuts.meta['GH_EFF']}"
-                )
+                if 'GH_EFF' in self.energy_dependent_gh_cuts.meta.keys():
+                    self.log.info(
+                        "Using gamma efficiency of "
+                        f"{self.energy_dependent_gh_cuts.meta['GH_EFF']}"
+                    )
+                else:
+                    self.log.info(
+                        "Using energy dependent gammaness cuts : "
+                        f"{self.energy_dependent_gh_cuts.meta['GH_CUTS']}"
+                    )
             else:
                 self.cuts.global_gh_cut = self.irf_final_hdu[1].header["GH_CUT"]
                 data_temp = self.cuts.apply_global_gh_cut(data_temp)
@@ -406,10 +418,16 @@ def apply_srcdep_gh_alpha_cut(self):
                 data_temp = self.cuts.apply_energy_dependent_alpha_cuts(
                     data_temp, self.energy_dependent_alpha_cuts
                 )
-                self.log.info(
-                    "Using alpha containment region of "
-                    f'{self.energy_dependent_alpha_cuts.meta["AL_CONT"]}'
-                )
+                if 'AL_CONT' in self.energy_dependent_alpha_cuts.meta.keys():
+                    self.log.info(
+                        "Using alpha containment region of "
+                        f'{self.energy_dependent_alpha_cuts.meta["AL_CONT"]}'
+                    )
+                else:
+                    self.log.info(
+                        "Using energy dependent gammaness cuts : "
+                        f"{self.energy_dependent_gh_cuts.meta['AL_CUTS']}"
+                    )
             else:
                 self.cuts.global_alpha_cut = self.irf_final_hdu[1].header["AL_CUT"]
                 data_temp = self.cuts.apply_global_alpha_cut(data_temp)