Merge pull request #280 from cta-observatory/edep_cuts_irfs

Implement custom energy-dependent cuts for IRFs and DL3 creations

Author: jsitarek

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 strategy.
+        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,41 @@ 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)
+        for cut_table in (cut_table_gh, cut_table_theta):
+            assert isinstance(cut_table, astropy.table.QTable)
+            assert cut_table["low"].unit == "TeV"
+            assert cut_table["high"].unit == "TeV"
+            assert cut_table["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"][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"

magicctapipe/scripts/lst1_magic/lst1_magic_create_irf.py

@@ -23,14 +23,17 @@
 The "software(_only_3tel)" types are supposed for the software event
 coincidence with LST-mono and MAGIC-stereo observations, allowing for
 only the events triggering both M1 and M2. The "software" type allows
-for the events of the any 2-tel combinations except the MAGIC-stereo
+for the events of any 2-tel combinations except the MAGIC-stereo
 combination at the moment. The "software_only_3tel" type allows for only
 the events of the 3-tel combination. The "magic_only" type allows for
 only the events of the MAGIC-stereo combination.
 
-There are two types of gammaness and theta cuts, "global" and "dynamic".
+There are three types of gammaness and theta cuts, "global", "dynamic" and
+"custom".
 In case of the dynamic cuts, the optimal cut satisfying a given
 efficiency will be calculated for every energy bin.
+In case of custom cuts, the provided cuts with be applied to the
+target energy binning using the selected scipy.interpolate.interpd1d method.
 
 Usage:
 $ python lst1_magic_create_irf.py
@@ -74,7 +77,12 @@
     calculate_event_weights,
 )
 
-from magicctapipe.io import create_gh_cuts_hdu, format_object, load_mc_dl2_data_file
+from magicctapipe.io import (
+    create_gh_cuts_hdu,
+    format_object,
+    get_custom_cuts,
+    load_mc_dl2_data_file,
+)
 
 __all__ = ["create_irf"]
 
@@ -110,6 +118,9 @@ def create_irf(
         If the input type of gammaness or theta cut is not known
     """
 
+    dynamic_keys = ["efficiency", "min_cut", "max_cut"]
+    custom_keys = ["custom_cuts", "interpolate_kind"]
+
     config_irf = config["create_irf"]
 
     quality_cuts = config_irf["quality_cuts"]
@@ -336,10 +347,12 @@ def create_irf(
             event_table_bkg = event_table_bkg[mask_gh]
 
     elif cut_type_gh == "dynamic":
-        config_gh_cuts.pop("global_cut_value", None)
-
         logger.info("\nDynamic gammaness cuts:")
-        logger.info(format_object(config_gh_cuts))
+        logger.info(
+            format_object(
+                {k: config_gh_cuts[k] for k in dynamic_keys if k in config_gh_cuts}
+            )
+        )
 
         gh_efficiency = config_gh_cuts["efficiency"]
         gh_cut_min = config_gh_cuts["min_cut"]
@@ -364,6 +377,20 @@ def create_irf(
             max_value=gh_cut_max,
         )
 
+    elif cut_type_gh == "custom":
+        cut_table_gh = get_custom_cuts(config_gh_cuts, energy_bins)
+        extra_header["GH_CUTS"] = "Custom energy-dependent"
+        output_suffix = "gh_custom"
+        logger.info("\nCustom gammaness cuts.")
+        logger.info(
+            format_object(
+                {k: config_gh_cuts[k] for k in custom_keys if k in config_gh_cuts}
+            )
+        )
+    else:
+        raise ValueError(f"Unknown gammaness-cut type '{cut_type_gh}'.")
+
+    if cut_type_gh == "dynamic" or "custom":
         logger.info(f"\nGammaness-cut table:\n\n{cut_table_gh}")
 
         # Apply the dynamic gammaness cuts
@@ -401,9 +428,6 @@ def create_irf(
 
         irf_hdus.append(hdu_gh_cuts)
 
-    else:
-        raise ValueError(f"Unknown gammaness-cut type '{cut_type_gh}'.")
-
     if is_point_like:
         # Apply the theta cut
         config_theta_cuts = config_irf["theta"]
@@ -423,10 +447,17 @@ def create_irf(
             event_table_gamma = event_table_gamma[mask_theta]
 
         elif cut_type_theta == "dynamic":
-            config_theta_cuts.pop("global_cut_value", None)
 
             logger.info("\nDynamic theta cuts:")
-            logger.info(format_object(config_theta_cuts))
+            logger.info(
+                format_object(
+                    {
+                        k: config_theta_cuts[k]
+                        for k in dynamic_keys
+                        if k in config_theta_cuts
+                    }
+                )
+            )
 
             theta_efficiency = config_theta_cuts["efficiency"]
             theta_cut_min = u.Quantity(config_theta_cuts["min_cut"])
@@ -450,7 +481,24 @@ def create_irf(
                 min_value=theta_cut_min,
                 max_value=theta_cut_max,
             )
+        elif cut_type_theta == "custom":
+            cut_table_theta = get_custom_cuts(config_theta_cuts, energy_bins)
+            extra_header["THETA_CUTS"] = "Custom energy-dependent"
+            output_suffix += "_theta_custom"
+            logger.info("\nCustom theta cuts.")
+            logger.info(
+                format_object(
+                    {
+                        k: config_theta_cuts[k]
+                        for k in custom_keys
+                        if k in config_theta_cuts
+                    }
+                )
+            )
+        else:
+            raise ValueError(f"Unknown theta-cut type '{cut_type_theta}'.")
 
+        if cut_type_theta == "dynamic" or "custom":
             logger.info(f"\nTheta-cut table:\n\n{cut_table_theta}")
 
             # Apply the dynamic theta cuts
@@ -479,9 +527,6 @@ def create_irf(
 
             irf_hdus.append(hdu_rad_max)
 
-        else:
-            raise ValueError(f"Unknown theta-cut type '{cut_type_theta}'.")
-
     # Create an effective-area HDU
     logger.info("\nCreating an effective-area HDU...")