Add possibility to produce DL3 with ctapipe

.codespell-ignores

@@ -3,3 +3,4 @@ usera
 nd
 studi
 referenc
+livetime

pyproject.toml

@@ -108,6 +108,7 @@ ctapipe-process = "ctapipe.tools.process:main"
 ctapipe-merge = "ctapipe.tools.merge:main"
 ctapipe-optimize-event-selection = "ctapipe.tools.optimize_event_selection:main"
 ctapipe-compute-irf = "ctapipe.tools.compute_irf:main"
+ctapipe-create-dl3 = "ctapipe.tools.create_dl3:main"
 ctapipe-fileinfo = "ctapipe.tools.fileinfo:main"
 ctapipe-quickstart = "ctapipe.tools.quickstart:main"
 ctapipe-calculate-pixel-statistics = "ctapipe.tools.calculate_pixel_stats:main"

src/ctapipe/conftest.py

@@ -783,7 +783,7 @@ def irf_event_loader_test_config():
                 "energy_reconstructor": "ExtraTreesRegressor",
                 "geometry_reconstructor": "HillasReconstructor",
                 "gammaness_classifier": "ExtraTreesClassifier",
-                "EventQualityQuery": {
+                "EventQualitySelection": {
                     "quality_criteria": [
                         (
                             "multiplicity 4",
@@ -816,9 +816,17 @@ def irf_events_table():
     N2 = 100
     N = N1 + N2
     epp = EventPreprocessor()
-    tab = epp.make_empty_table()
-
-    ids, bulk, unitless = tab.colnames[:2], tab.colnames[2:-2], tab.colnames[-2:]
+    keep_columns, _, _ = epp.get_columns_keep_rename_scheme(None, True)
+    tab = epp.make_empty_table(keep_columns)
+
+    ids, bulk, unitless = [], [], []
+    for c in tab.columns:
+        if "id" in c:
+            ids.append(c)
+        elif tab[c].unit == u.dimensionless_unscaled:
+            unitless.append(c)
+        else:
+            bulk.append(c)
 
     id_tab = QTable(
         data=np.zeros((N, len(ids)), dtype=np.uint64),

src/ctapipe/irf/__init__.py

@@ -18,6 +18,7 @@
     check_bins_in_range,
     make_bins_per_decade,
 )
+from .cuts import EventQualitySelection, EventSelection
 from .irfs import (
     BackgroundRate2dMaker,
     EffectiveArea2dMaker,
@@ -56,4 +57,6 @@
     "FLUX_UNIT",
     "check_bins_in_range",
     "make_bins_per_decade",
+    "EventSelection",
+    "EventQualitySelection",
 ]

src/ctapipe/irf/cuts/__init__.py

@@ -0,0 +1,4 @@
+from .quality_cuts import EventQualitySelection
+from .selection_cuts import EventSelection
+
+__all__ = ["EventSelection", "EventQualitySelection"]

src/ctapipe/irf/cuts/quality_cuts.py

@@ -0,0 +1,58 @@
+from astropy.table import Table
+
+from ...core import QualityQuery
+from ...core.traits import List, Tuple, Unicode
+
+
+class EventQualitySelection(QualityQuery):
+    """
+    Event pre-selection quality criteria for IRF and DL3 computation with different defaults.
+    """
+
+    quality_criteria = List(
+        Tuple(Unicode(), Unicode()),
+        default_value=[
+            (
+                "multiplicity 4",
+                "np.count_nonzero(HillasReconstructor_telescopes,axis=1) >= 4",
+            ),
+            ("valid classifier", "RandomForestClassifier_is_valid"),
+            ("valid geom reco", "HillasReconstructor_is_valid"),
+            ("valid energy reco", "RandomForestRegressor_is_valid"),
+        ],
+        help=QualityQuery.quality_criteria.help,
+    ).tag(config=True)
+
+    def calculate_selection(self, events: Table):
+        """
+        Add the selection columns to the events, will only compute quality selection
+
+        Parameters
+        ----------
+        events: Table
+            The table containing the events on which selection need to be applied
+
+        Returns
+        -------
+        Table
+            events with selection columns added.
+        """
+        return self.calculate_quality_selection(events)
+
+    def calculate_quality_selection(self, events: Table):
+        """
+        Add the selection columns to the events, will only compute quality selection
+
+        Parameters
+        ----------
+        events: Table
+            The table containing the events on which selection need to be applied
+
+        Returns
+        -------
+        Table
+            events with selection columns added.
+        """
+        events["selected_quality"] = self.get_table_mask(events)
+        events["selected"] = events["selected_quality"]
+        return events

src/ctapipe/irf/cuts/selection_cuts.py

@@ -0,0 +1,88 @@
+import operator
+
+from astropy.table import Table
+from pyirf.cuts import evaluate_binned_cut
+
+from ...core.traits import Path
+from ..optimize.results import OptimizationResult
+from .quality_cuts import EventQualitySelection
+
+__all__ = ["EventSelection"]
+
+
+class EventSelection(EventQualitySelection):
+    """
+    Event selection quality and gammaness criteria for IRF and DL3
+    """
+
+    cuts_file = Path(
+        allow_none=False,
+        directory_ok=False,
+        exists=True,
+        help="Path to the cuts file to apply to the observation.",
+    ).tag(config=True)
+
+    def __init__(self, config=None, parent=None, **kwargs):
+        super().__init__(config=config, parent=parent, **kwargs)
+        self.cuts = OptimizationResult.read(self.cuts_file)
+
+    def calculate_selection(self, events: Table, apply_spatial_selection: bool = False):
+        """
+        Add the selection columns to the events
+
+        Parameters
+        ----------
+        events: Table
+            The table containing the events on which selection need to be applied
+        apply_spatial_selection: bool
+            True if the theta cuts should be applied
+
+        Returns
+        -------
+        Table
+            events with selection columns added.
+        """
+        events = self.calculate_quality_selection(events)
+        events = self.calculate_gamma_selection(events, apply_spatial_selection)
+        events["selected"] = events["selected_quality"] & events["selected_gamma"]
+        return events
+
+    def calculate_gamma_selection(
+        self, events: Table, apply_spatial_selection: bool = False
+    ) -> Table:
+        """
+        Add the selection columns to the events, will compute only gamma criteria
+
+        Parameters
+        ----------
+        events: Table
+            The table containing the events on which selection need to be applied
+        apply_spatial_selection: bool
+            True if the theta cuts should be applied
+
+        Returns
+        -------
+        Table
+            events with selection columns added.
+        """
+
+        events["selected_gh"] = evaluate_binned_cut(
+            events["gh_score"],
+            events["reco_energy"],
+            self.cuts.gh_cuts,
+            operator.ge,
+        )
+
+        if apply_spatial_selection:
+            events["selected_theta"] = evaluate_binned_cut(
+                events["theta"],
+                events["reco_energy"],
+                self.cuts.spatial_selection_table,
+                operator.le,
+            )
+            events["selected_gamma"] = events["selected_theta"] & events["selected_gh"]
+        else:
+            events["selected_gamma"] = events["selected_gh"]
+
+        events["selected"] = events["selected_gamma"]
+        return events