Add full cut optimization as introduced in pyirf 0.13

docs/changes/2789.feature.rst

@@ -0,0 +1,3 @@
+Add support for the full cut optimization introduced in pyirf 0.13.
+This can now be used via the ``PointSourceSensitvityOptimizer``,
+while the previous EventDisplay-like optimization can be used via the ``PointSourceSensitvityGhOptimizer``.

pyproject.toml

@@ -58,7 +58,7 @@ all = [
     "ctapipe[eventio]",
     "iminuit >=2",
     "matplotlib ~=3.0",
-    "pyirf ~=0.13.0"
+    "pyirf >=0.13.0",
 ]
 
 tests = [

src/ctapipe/core/tests/test_traits.py

@@ -333,7 +333,7 @@ class SomeBadComponentWithEnergyTrait(Component):
     with pytest.raises(
         TraitError,
         match=f"Given physical type {u.physical.energy} does not match"
-        + f" physical type of the default value, {u.get_physical_type(5 * u.m)}.",
+        + f" the default value's physical type {u.get_physical_type(5 * u.m)}.",
     ):
 
         class AnotherBadComponentWithEnergyTrait(Component):
@@ -347,7 +347,9 @@ def test_quantity_from_config(tmp_path):
     from ctapipe.core.traits import AstroQuantity
 
     class QuantityComponent(Component):
-        distance = AstroQuantity(u.cm, default_value=1 * u.m).tag(config=True)
+        distance = AstroQuantity(physical_type=u.cm, default_value=1 * u.m).tag(
+            config=True
+        )
 
     config = Config()
     config.QuantityComponent.distance = "5 cm"
@@ -362,7 +364,9 @@ class QuantityComponent(Component):
     assert q.distance.unit == u.m
 
     class QTool(Tool):
-        distance = AstroQuantity(u.cm, default_value=1 * u.m).tag(config=True)
+        distance = AstroQuantity(physical_type=u.cm, default_value=1 * u.m).tag(
+            config=True
+        )
 
     config_file = tmp_path / "config.json"
     config_file.write_text(
@@ -398,7 +402,9 @@ class MyTool(Tool):
 
 def test_quantity_none():
     class AllowNone(Component):
-        quantity = AstroQuantity(default_value=None, allow_none=True)
+        quantity = AstroQuantity(
+            default_value=None, allow_none=True, physical_type=u.TeV
+        )
 
     c = AllowNone()
     assert c.quantity is None

src/ctapipe/core/traits.py

@@ -97,27 +97,34 @@ class AstroQuantity(TraitType):
 
     def __init__(self, physical_type=None, **kwargs):
         super().__init__(**kwargs)
-        if physical_type is not None:
-            if isinstance(physical_type, u.PhysicalType):
-                self.physical_type = physical_type
-            elif isinstance(physical_type, u.UnitBase):
-                self.physical_type = u.get_physical_type(physical_type)
-            else:
-                raise TraitError(
-                    "Given physical type must be either of type"
-                    " astropy.units.PhysicalType or a subclass of"
-                    f" astropy.units.UnitBase, was {type(physical_type)}."
-                )
-        else:
+        if physical_type is None:
+            self.physical_type = None
+        elif isinstance(physical_type, u.PhysicalType):
             self.physical_type = physical_type
+        elif isinstance(physical_type, u.UnitBase):
+            self.physical_type = u.get_physical_type(physical_type)
+        else:
+            raise TraitError(
+                "Given physical type must be either of type"
+                " astropy.units.PhysicalType or a subclass of"
+                f" astropy.units.UnitBase, was {type(physical_type)}."
+            )
+
+        if self.default_value is not Undefined and self.default_value is not None:
+            self._validate_default_value()
 
-        if self.default_value is not Undefined and self.physical_type is not None:
+    def _validate_default_value(self):
+        if self.physical_type is not None:
             default_type = u.get_physical_type(self.default_value)
             if default_type != self.physical_type:
                 raise TraitError(
                     f"Given physical type {self.physical_type} does not match"
-                    f" physical type of the default value, {default_type}."
+                    f" the default value's physical type {default_type}."
                 )
+        else:
+            if not isinstance(self.default_value, u.Quantity):
+                self.default_value = u.Quantity(self.default_value)
+            self.physical_type = u.get_physical_type(self.default_value)
 
     @property
     def info_text(self):

src/ctapipe/io/dl2_tables_preprocessing.py

@@ -42,10 +42,6 @@ class DL2EventQualityQuery(QualityQuery):
     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"),
@@ -57,8 +53,6 @@ class DL2EventQualityQuery(QualityQuery):
 class DL2EventPreprocessor(Component):
     """Defines pre-selection cuts and the necessary renaming of columns."""
 
-    classes = [DL2EventQualityQuery]
-
     energy_reconstructor = Unicode(
         default_value="RandomForestRegressor",
         help="Prefix of the reco `_energy` column",
@@ -177,6 +171,13 @@ def normalise_column_names(self, events: QTable) -> QTable:
 
         fixed_columns = list(set(columns_to_keep) - set(rename_to))
         columns_to_read = fixed_columns + rename_from
+        if self.apply_derived_columns:
+            # read columns needed for "multiplicity"
+            columns_to_read += [
+                f"{self.energy_reconstructor}_telescopes",
+                f"{self.gammaness_classifier}_telescopes",
+                f"{self.geometry_reconstructor}_telescopes",
+            ]
         for col in columns_to_read:
             if col not in events.colnames:
                 raise ValueError(
@@ -230,6 +231,11 @@ def make_empty_table(self) -> QTable:
                         dtype=np.float64,
                         description="Event weight",
                     ),
+                    Column(
+                        name="multiplicity",
+                        dtype=np.int16,
+                        description="Event multiplicity",
+                    ),
                 ]
             )
 
@@ -250,8 +256,6 @@ class DL2EventLoader(Component):
     Component for loading events and simulation metadata, applying preselection and optional derived column logic.
     """
 
-    classes = [DL2EventPreprocessor]
-
     # User-selectable event reading function and kwargs
     event_reader_function = Unicode(
         default_value="read_subarray_events_chunked",
@@ -436,6 +440,29 @@ def make_derived_columns(self, events: QTable) -> QTable:
         events["reco_fov_lon"] = u.Quantity(-reco_nominal.fov_lon)  # minus for GADF
         events["reco_fov_lat"] = u.Quantity(reco_nominal.fov_lat)
         events["weight"] = 1.0  # defer calculation of proper weights to later
+        # define multiplicity as lowest multiplicity between the 3 reconstructions
+        events["multiplicity"] = np.min(
+            [
+                np.count_nonzero(
+                    events[f"{self.epp.energy_reconstructor}_telescopes"], axis=1
+                ),
+                np.count_nonzero(
+                    events[f"{self.epp.geometry_reconstructor}_telescopes"], axis=1
+                ),
+                np.count_nonzero(
+                    events[f"{self.epp.gammaness_classifier}_telescopes"], axis=1
+                ),
+            ],
+            axis=0,
+        )
+        # delete "_telescope" columns as they are not needed downstream
+        events.remove_columns(
+            [
+                f"{self.epp.energy_reconstructor}_telescopes",
+                f"{self.epp.geometry_reconstructor}_telescopes",
+                f"{self.epp.gammaness_classifier}_telescopes",
+            ]
+        )
         return events
 
     def make_event_weights(
@@ -469,6 +496,11 @@ def make_event_weights(
         ValueError
             If ``fov_offset_bins`` is required but not provided.
         """
+        # Re-initialize weights = 0 to exclude events outside fov bins
+        # in the cut optimization.
+        # This will become unnecessary once fov bins are used during cut optimization.
+        events["weight"] = 0.0
+
         if (
             kind == "gammas"
             and self.target_spectrum.normalization.unit.is_equivalent(

src/ctapipe/io/tests/test_preprocessing.py

@@ -14,11 +14,35 @@ def dummy_table():
             "event_id": [1, 2, 3, 1, 1, 2],
             "true_energy": [0.99, 10, 0.37, 2.1, 73.4, 1] * u.TeV,
             "dummy_energy": [1, 10, 0.4, 2.5, 73, 1] * u.TeV,
+            "dummy_telescopes": [
+                [True, True, True, False, False, True],
+                [True, True, True, True, False, True],
+                [False, True, True, True, True, False],
+                [False, False, False, False, True, True],
+                [True, True, True, True, True, True],
+                [True, False, True, False, True, True],
+            ],
             "classifier_prediction": [1, 0.3, 0.87, 0.93, 0, 0.1],
+            "classifier_telescopes": [
+                [True, True, True, False, False, True],
+                [True, True, True, True, False, True],
+                [False, True, True, True, True, False],
+                [False, False, False, False, True, True],
+                [True, True, True, True, True, True],
+                [True, False, True, False, True, True],
+            ],
             "true_alt": [60, 60, 60, 60, 60, 60] * u.deg,
             "geom_alt": [58.5, 61.2, 59, 71.6, 60, 62] * u.deg,
             "true_az": [13, 13, 13, 13, 13, 13] * u.deg,
             "geom_az": [12.5, 13, 11.8, 15.1, 14.7, 12.8] * u.deg,
+            "geom_telescopes": [
+                [True, True, True, False, False, True],
+                [True, True, True, True, False, True],
+                [False, True, True, True, True, False],
+                [False, False, False, False, True, True],
+                [True, True, True, True, True, True],
+                [True, False, True, False, True, True],
+            ],
             "subarray_pointing_frame": np.zeros(6),
             "subarray_pointing_lat": np.full(6, 20) * u.deg,
             "subarray_pointing_lon": np.full(6, 0) * u.deg,
@@ -119,6 +143,7 @@ def test_event_loader(gamma_diffuse_full_reco_file, irf_event_loader_test_config
         "true_source_fov_offset",
         "reco_source_fov_offset",
         "weight",
+        "multiplicity",
     ]
 
     assert sorted(columns) == sorted(events.colnames)
@@ -226,9 +251,12 @@ def test_name_overriding(dummy_table):
         "true_az",
         "true_alt",
         "dummy_energy",
+        "dummy_telescopes",
         "classifier_prediction",
+        "classifier_telescopes",
         "geom_alt",
         "geom_az",
+        "geom_telescopes",
         "subarray_pointing_frame",
         "subarray_pointing_lat",
         "subarray_pointing_lon",

src/ctapipe/irf/__init__.py

@@ -28,6 +28,7 @@
     GhPercentileCutCalculator,
     OptimizationResult,
     PercentileCuts,
+    PointSourceSensitivityGhOptimizer,
     PointSourceSensitivityOptimizer,
     ThetaPercentileCutCalculator,
 )
@@ -43,6 +44,7 @@
     "EffectiveArea2dMaker",
     "ResultValidRange",
     "OptimizationResult",
+    "PointSourceSensitivityGhOptimizer",
     "PointSourceSensitivityOptimizer",
     "PercentileCuts",
     "Spectra",