From 79448df794386b6d5f2ef1f78f919a0c32f4c792 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Wed, 18 Jun 2025 18:46:16 +0200 Subject: [PATCH 01/16] Start adding full cut opt from pyirf 0.13.0; remove unnecessay 'classes' attributes --- pyproject.toml | 2 +- src/ctapipe/io/dl2_tables_preprocessing.py | 47 +++++- src/ctapipe/io/tests/test_preprocessing.py | 24 +++ src/ctapipe/irf/__init__.py | 2 + src/ctapipe/irf/optimize.py | 185 ++++++++++++++++----- src/ctapipe/tools/compute_irf.py | 39 ++++- 6 files changed, 243 insertions(+), 56 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 33bee31c21c..ad75eee4f94 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -58,7 +58,7 @@ all = [ "ctapipe[eventio]", "iminuit >=2", "matplotlib ~=3.0", - "pyirf ~=0.14.0" + "pyirf ~=0.14.0", ] tests = [ diff --git a/src/ctapipe/io/dl2_tables_preprocessing.py b/src/ctapipe/io/dl2_tables_preprocessing.py index c0776496d75..2181a778aa8 100644 --- a/src/ctapipe/io/dl2_tables_preprocessing.py +++ b/src/ctapipe/io/dl2_tables_preprocessing.py @@ -27,7 +27,7 @@ from ..compat import COPY_IF_NEEDED from ..containers import CoordinateFrameType from ..coordinates import NominalFrame -from ..core import Component, QualityQuery +from ..core import Component, QualityQuery, ToolConfigurationError from ..core.traits import Bool, Dict, List, Tuple, Unicode from .tableloader import TableLoader @@ -57,8 +57,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 +175,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 +235,11 @@ def make_empty_table(self) -> QTable: dtype=np.float64, description="Event weight", ), + Column( + name="multiplicity", + dtype=np.int16, + description="Event multiplicity", + ), ] ) @@ -250,8 +260,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 +444,35 @@ 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 + events["multiplicity"] = np.count_nonzero( + events[f"{self.epp.energy_reconstructor}_telescopes"], axis=1 + ) + if not ( + np.array_equal( + events["multiplicity"], + np.count_nonzero( + events[f"{self.epp.geometry_reconstructor}_telescopes"], axis=1 + ), + ) + and np.array_equal( + events["multiplicity"], + np.count_nonzero( + events[f"{self.epp.gammaness_classifier}_telescopes"], axis=1 + ), + ) + ): + raise ToolConfigurationError( + "There are events for which not all the reconstructions were successful. " + "Please check your configuration of the EventQualityQuery." + ) + # 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( diff --git a/src/ctapipe/io/tests/test_preprocessing.py b/src/ctapipe/io/tests/test_preprocessing.py index 46879ac6350..9520634f1ee 100644 --- a/src/ctapipe/io/tests/test_preprocessing.py +++ b/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, diff --git a/src/ctapipe/irf/__init__.py b/src/ctapipe/irf/__init__.py index d4abdb43fac..42a40af2c19 100644 --- a/src/ctapipe/irf/__init__.py +++ b/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", diff --git a/src/ctapipe/irf/optimize.py b/src/ctapipe/irf/optimize.py index 7e937f106f8..a10f283bda4 100644 --- a/src/ctapipe/irf/optimize.py +++ b/src/ctapipe/irf/optimize.py @@ -8,11 +8,11 @@ import numpy as np from astropy.io import fits from astropy.table import QTable, Table -from pyirf.cut_optimization import optimize_gh_cut +from pyirf.cut_optimization import optimize_cuts, optimize_gh_cut from pyirf.cuts import calculate_percentile_cut, evaluate_binned_cut from ..core import Component, QualityQuery -from ..core.traits import AstroQuantity, Float, Integer, Path +from ..core.traits import AstroQuantity, Float, Integer, List, Path from ..io.dl2_tables_preprocessing import DL2EventQualityQuery from .binning import DefaultRecoEnergyBins, ResultValidRange @@ -21,7 +21,9 @@ "GhPercentileCutCalculator", "OptimizationResult", "PercentileCuts", + "PointSourceSensitivityGhOptimizer", "PointSourceSensitivityOptimizer", + "PointSourceSensitivityOptimizerBase", "ThetaPercentileCutCalculator", ] @@ -37,8 +39,9 @@ def __init__( valid_offset_max: u.Quantity, gh_cuts: QTable, clf_prefix: str, - spatial_selection_table: QTable | None = None, quality_query: QualityQuery | Sequence | None = None, + spatial_selection_table: QTable | None = None, + multiplicity_cuts: QTable | None = None, ) -> None: if quality_query: if isinstance(quality_query, QualityQuery): @@ -60,23 +63,23 @@ def __init__( self.gh_cuts = gh_cuts self.clf_prefix = clf_prefix self.spatial_selection_table = spatial_selection_table + self.multiplicity_cuts = multiplicity_cuts def __repr__(self): + repr = f"" - ) - else: - return ( - f"" - ) + repr += f", and {len(self.spatial_selection_table)} theta bins" + + if self.multiplicity_cuts is not None: + repr += f", and {len(self.multiplicity_cuts)} multiplicity bins" + + repr += ( + f" valid between {self.valid_offset.min} to {self.valid_offset.max}" + f" and {self.valid_energy.min} to {self.valid_energy.max}" + f" with {len(self.quality_query.quality_criteria)} quality criteria>" + ) + return repr def write(self, output_name: Path | str, overwrite: bool = False) -> None: """Write an ``OptimizationResult`` to a file in FITS format.""" @@ -109,6 +112,10 @@ def write(self, output_name: Path | str, overwrite: bool = False) -> None: self.spatial_selection_table.meta["EXTNAME"] = "RAD_MAX" results.append(self.spatial_selection_table) + if self.multiplicity_cuts is not None: + self.multiplicity_cuts.meta["EXTNAME"] = "MULTIPLICITY_CUTS" + results.append(self.multiplicity_cuts) + primary = fits.PrimaryHDU() hdus = [primary] for table in results: @@ -129,7 +136,16 @@ def read(cls, file_name): gh_cuts = QTable.read(hdul[2]) valid_energy = QTable.read(hdul[3]) valid_offset = QTable.read(hdul[4]) - spatial_selection_table = QTable.read(hdul[5]) if len(hdul) > 5 else None + spatial_selection_table = None + multiplicity_cuts = None + if len(hdul) == 6: + if hdul[5].header["EXTNAME"] == "RAD_MAX": + spatial_selection_table = QTable.read(hdul[5]) + else: + multiplicity_cuts = QTable.read(hdul[5]) + elif len(hdul) == 7: + spatial_selection_table = QTable.read(hdul[5]) + multiplicity_cuts = QTable.read(hdul[6]) return cls( quality_query=quality_query, @@ -140,6 +156,7 @@ def read(cls, file_name): gh_cuts=gh_cuts, clf_prefix=gh_cuts.meta["CLFNAME"], spatial_selection_table=spatial_selection_table, + multiplicity_cuts=multiplicity_cuts, ) @@ -148,10 +165,8 @@ class CutOptimizerBase(DefaultRecoEnergyBins): needs_background = False - def __init__(self, config=None, parent=None, **kwargs): - super().__init__(config=config, parent=parent, **kwargs) - def _check_events(self, events: dict[str, QTable]): + """To be called at the start `__call__` in all subclasses.""" if "signal" not in events.keys(): raise ValueError( "Calculating G/H and/or spatial selection cuts requires 'signal' " @@ -295,8 +310,6 @@ class PercentileCuts(CutOptimizerBase): events surviving this G/H cut. """ - classes = [GhPercentileCutCalculator, ThetaPercentileCutCalculator] - def __init__(self, config=None, parent=None, **kwargs): super().__init__(config=config, parent=parent, **kwargs) self.gh_cut_calculator = GhPercentileCutCalculator(parent=self) @@ -341,18 +354,14 @@ def __call__( return result -class PointSourceSensitivityOptimizer(CutOptimizerBase): +class PointSourceSensitivityOptimizerBase(CutOptimizerBase): """ - Optimizes a G/H cut for maximum point source sensitivity and - calculates a percentile cut on theta. + Base class for optimizers optimizing cuts for maximum point source efficiency. + This holds some configurable traits shared between different optimization + algorithms implemented in subclasses of this class. """ needs_background = True - classes = [ThetaPercentileCutCalculator] - - initial_gh_cut_efficency = Float( - default_value=0.4, help="Start value of gamma efficiency before optimization" - ).tag(config=True) max_gh_cut_efficiency = Float( default_value=0.8, help="Maximum gamma efficiency requested" @@ -386,6 +395,30 @@ class PointSourceSensitivityOptimizer(CutOptimizerBase): physical_type=u.physical.angle, ).tag(config=True) + def _get_valid_energy_range(self, opt_sensitivity): + keep_mask = np.isfinite(opt_sensitivity["significance"]) + + count = np.arange(start=0, stop=len(keep_mask), step=1) + if all(np.diff(count[keep_mask]) == 1): + return [ + opt_sensitivity["reco_energy_low"][keep_mask][0], + opt_sensitivity["reco_energy_high"][keep_mask][-1], + ] + else: + raise ValueError("Optimal significance curve has internal NaN bins") + + +class PointSourceSensitivityGhOptimizer(PointSourceSensitivityOptimizerBase): + """ + EventDisplay-like procedure of optimizing a G/H cut for maximum point source + sensitivity, calculating only a percentile cut on theta, and no cut on + multiplicity. + """ + + initial_gh_cut_efficency = Float( + default_value=0.4, help="Start value of gamma efficiency before optimization" + ).tag(config=True) + def __init__(self, config=None, parent=None, **kwargs): super().__init__(config=config, parent=parent, **kwargs) self.theta_cut_calculator = ThetaPercentileCutCalculator(parent=self) @@ -466,14 +499,84 @@ def __call__( ) return result - def _get_valid_energy_range(self, opt_sensitivity): - keep_mask = np.isfinite(opt_sensitivity["significance"]) - count = np.arange(start=0, stop=len(keep_mask), step=1) - if all(np.diff(count[keep_mask]) == 1): - return [ - opt_sensitivity["reco_energy_low"][keep_mask][0], - opt_sensitivity["reco_energy_high"][keep_mask][-1], - ] - else: - raise ValueError("Optimal significance curve has internal NaN bins") +class PointSourceSensitivityOptimizer(PointSourceSensitivityOptimizerBase): + """ + Finds the combination of G/H cut, theta cut, and cut on the event-multiplicity + for each energy bin with the maximum point source sensitivity in a grid search. + """ + + theta_min_angle = AstroQuantity( + default_value=u.Quantity(-1, u.deg), + physical_type=u.physical.angle, + help="Smallest angular cut value allowed (-1 means no cut)", + ).tag(config=True) + + theta_max_angle = AstroQuantity( + default_value=u.Quantity(0.32, u.deg), + physical_type=u.physical.angle, + help="Largest angular cut value allowed", + ).tag(config=True) + + max_theta_cut_efficiency = Float( + default_value=0.8, help="Maximum gamma efficiency of the theta cut requested" + ).tag(config=True) + + theta_cut_efficiency_step = Float( + default_value=0.1, + help="Efficiency-stepsize used for scanning after optimal theta cut", + ).tag(config=True) + + multiplicity_cuts = List( + Integer(), + default_value=[2, 4, 6, 8, 10], + help="Event-multiplicity cuts used for scanning after optimal cut", + ).tag(config=True) + + def __call__( + self, + events: dict[str, QTable], + quality_query: EventQualityQuery, + clf_prefix: str, + ) -> OptimizationResult: + self._check_events(events) + + gh_cut_efficiencies = np.arange( + self.gh_cut_efficiency_step, + self.max_gh_cut_efficiency + self.gh_cut_efficiency_step / 2, + self.gh_cut_efficiency_step, + ) + theta_cut_efficiencies = np.arange( + self.theta_cut_efficiency_step, + self.max_theta_cut_efficiency + self.theta_cut_efficiency_step / 2, + self.theta_cut_efficiency_step, + ) + + opt_sensitivity, multiplicity_cuts, theta_cuts, gh_cuts = optimize_cuts( + events["signal"], + events["background"], + reco_energy_bins=self.reco_energy_bins, + multiplicity_cuts=self.multiplicity_cuts, + gh_cut_efficiencies=gh_cut_efficiencies, + theta_cut_efficiencies=theta_cut_efficiencies, + theta_min_value=self.theta_min_angle, + theta_max_value=self.theta_max_angle, + alpha=self.alpha, + fov_offset_max=self.max_background_fov_offset, + fov_offset_min=self.min_background_fov_offset, + ) + valid_energy = self._get_valid_energy_range(opt_sensitivity) + + result = OptimizationResult( + quality_query=quality_query, + clf_prefix=clf_prefix, + gh_cuts=gh_cuts, + valid_energy_min=valid_energy[0], + valid_energy_max=valid_energy[1], + # A single set of cuts is calculated for the whole fov atm + valid_offset_min=self.min_background_fov_offset, + valid_offset_max=self.max_background_fov_offset, + spatial_selection_table=theta_cuts, + multiplicity_cuts=multiplicity_cuts, + ) + return result diff --git a/src/ctapipe/tools/compute_irf.py b/src/ctapipe/tools/compute_irf.py index dd47f8755c5..6b954537acd 100644 --- a/src/ctapipe/tools/compute_irf.py +++ b/src/ctapipe/tools/compute_irf.py @@ -355,6 +355,7 @@ def calculate_selections(self, reduced_events: dict) -> dict: self.opt_result.gh_cuts, operator.ge, ) + reduced_events["gammas"]["selected"] = reduced_events["gammas"]["selected_gh"] if self.spatial_selection_applied: reduced_events["gammas"]["selected_theta"] = evaluate_binned_cut( reduced_events["gammas"]["theta"], @@ -362,13 +363,19 @@ def calculate_selections(self, reduced_events: dict) -> dict: self.opt_result.spatial_selection_table, operator.le, ) - reduced_events["gammas"]["selected"] = ( - reduced_events["gammas"]["selected_theta"] - & reduced_events["gammas"]["selected_gh"] + reduced_events["gammas"]["selected"] &= reduced_events["gammas"][ + "selected_theta" + ] + + if self.opt_result.multiplicity_cuts is not None: + reduced_events["gammas"]["selected_multiplicity"] = evaluate_binned_cut( + reduced_events["gammas"]["multiplicity"], + reduced_events["gammas"]["reco_energy"], + self.opt_result.multiplicity_cuts, + operator.ge, ) - else: - reduced_events["gammas"]["selected"] = reduced_events["gammas"][ - "selected_gh" + reduced_events["gammas"]["selected"] &= reduced_events["gammas"][ + "selected_multiplicity" ] if self.do_background: @@ -383,9 +390,23 @@ def calculate_selections(self, reduced_events: dict) -> dict: self.opt_result.gh_cuts, operator.ge, ) - n_sel[bkg_type] = np.count_nonzero( - reduced_events[bkg_type]["selected_gh"] - ) + reduced_events[bkg_type]["selected"] = reduced_events[bkg_type][ + "selected_gh" + ] + if self.opt_result.multiplicity_cuts is not None: + reduced_events[bkg_type][ + "selected_multiplicity" + ] = evaluate_binned_cut( + reduced_events[bkg_type]["multiplicity"], + reduced_events[bkg_type]["reco_energy"], + self.opt_result.multiplicity_cuts, + operator.ge, + ) + reduced_events[bkg_type]["selected"] &= reduced_events[bkg_type][ + "selected_multiplicity" + ] + + n_sel[bkg_type] = np.count_nonzero(reduced_events[bkg_type]["selected"]) self.log.info( "Keeping %d signal, %d proton events, and %d electron events" From db5b249c7fba0618960562b9b7e188cfbb75117a Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Thu, 26 Jun 2025 18:17:29 +0200 Subject: [PATCH 02/16] Update tests --- src/ctapipe/io/tests/test_preprocessing.py | 1 + src/ctapipe/irf/tests/test_optimize.py | 6 +++++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/io/tests/test_preprocessing.py b/src/ctapipe/io/tests/test_preprocessing.py index 9520634f1ee..3aaf2bf577b 100644 --- a/src/ctapipe/io/tests/test_preprocessing.py +++ b/src/ctapipe/io/tests/test_preprocessing.py @@ -143,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) diff --git a/src/ctapipe/irf/tests/test_optimize.py b/src/ctapipe/irf/tests/test_optimize.py index b4e17818dd0..8170b97e92f 100644 --- a/src/ctapipe/irf/tests/test_optimize.py +++ b/src/ctapipe/irf/tests/test_optimize.py @@ -4,7 +4,7 @@ from astropy.table import QTable from ctapipe.core import QualityQuery, non_abstract_children -from ctapipe.irf.optimize import CutOptimizerBase +from ctapipe.irf.optimize import CutOptimizerBase, PointSourceSensitivityOptimizer @pytest.mark.parametrize("file_format", [".fits", ".fits.gz"]) @@ -121,3 +121,7 @@ def test_cut_optimizer( assert result.valid_energy.min >= result.gh_cuts["low"][0] assert result.valid_energy.max <= result.gh_cuts["high"][-1] assert result.spatial_selection_table["cut"].unit == u.deg + + if isinstance(optimizer, PointSourceSensitivityOptimizer): + assert result.multiplicity_cuts is not None + assert len(result.multiplicity_cuts["low"]) == len(result.gh_cuts["low"]) From fe23ded557e3eecedf20d581761fec67a8e7fdbd Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Thu, 26 Jun 2025 19:44:33 +0200 Subject: [PATCH 03/16] Fix multiplicity computation; remove multiplicity precut --- src/ctapipe/io/dl2_tables_preprocessing.py | 30 ++++++++-------------- src/ctapipe/irf/optimize.py | 2 +- 2 files changed, 11 insertions(+), 21 deletions(-) diff --git a/src/ctapipe/io/dl2_tables_preprocessing.py b/src/ctapipe/io/dl2_tables_preprocessing.py index 2181a778aa8..715ede92aed 100644 --- a/src/ctapipe/io/dl2_tables_preprocessing.py +++ b/src/ctapipe/io/dl2_tables_preprocessing.py @@ -27,7 +27,7 @@ from ..compat import COPY_IF_NEEDED from ..containers import CoordinateFrameType from ..coordinates import NominalFrame -from ..core import Component, QualityQuery, ToolConfigurationError +from ..core import Component, QualityQuery from ..core.traits import Bool, Dict, List, Tuple, Unicode from .tableloader import TableLoader @@ -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"), @@ -444,27 +440,21 @@ 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 - events["multiplicity"] = np.count_nonzero( - events[f"{self.epp.energy_reconstructor}_telescopes"], axis=1 - ) - if not ( - np.array_equal( - events["multiplicity"], + # 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 ), - ) - and np.array_equal( - events["multiplicity"], np.count_nonzero( events[f"{self.epp.gammaness_classifier}_telescopes"], axis=1 ), - ) - ): - raise ToolConfigurationError( - "There are events for which not all the reconstructions were successful. " - "Please check your configuration of the EventQualityQuery." - ) + ], + axis=0, + ) # delete "_telescope" columns as they are not needed downstream events.remove_columns( [ diff --git a/src/ctapipe/irf/optimize.py b/src/ctapipe/irf/optimize.py index a10f283bda4..13878e135f3 100644 --- a/src/ctapipe/irf/optimize.py +++ b/src/ctapipe/irf/optimize.py @@ -529,7 +529,7 @@ class PointSourceSensitivityOptimizer(PointSourceSensitivityOptimizerBase): multiplicity_cuts = List( Integer(), - default_value=[2, 4, 6, 8, 10], + default_value=[2, 3, 4, 5, 6, 7, 8, 9, 10], help="Event-multiplicity cuts used for scanning after optimal cut", ).tag(config=True) From 2311a7aa7a3eed448ba2bbe63086dd768b7cfba0 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Thu, 26 Jun 2025 20:09:12 +0200 Subject: [PATCH 04/16] Add changelog --- docs/changes/2789.feature.rst | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 docs/changes/2789.feature.rst diff --git a/docs/changes/2789.feature.rst b/docs/changes/2789.feature.rst new file mode 100644 index 00000000000..88366a1cf3e --- /dev/null +++ b/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``. From ce0e5adb25a33906b8e51f948cf5f6f353847028 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Fri, 27 Jun 2025 18:39:54 +0200 Subject: [PATCH 05/16] Index by extname when reading an OptimizationResult --- src/ctapipe/irf/optimize.py | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/src/ctapipe/irf/optimize.py b/src/ctapipe/irf/optimize.py index 13878e135f3..436976c3ee2 100644 --- a/src/ctapipe/irf/optimize.py +++ b/src/ctapipe/irf/optimize.py @@ -127,25 +127,21 @@ def read(cls, file_name): """Read an ``OptimizationResult`` from a file in FITS format.""" with fits.open(file_name) as hdul: - cut_expr_tab = Table.read(hdul[1]) + cut_expr_tab = Table.read(hdul["QUALITY_CUTS_EXPR"]) cut_expr_lst = [(name, expr) for name, expr in cut_expr_tab.iterrows()] if (" ", " ") in cut_expr_lst: cut_expr_lst.remove((" ", " ")) quality_query = QualityQuery(quality_criteria=cut_expr_lst) - gh_cuts = QTable.read(hdul[2]) - valid_energy = QTable.read(hdul[3]) - valid_offset = QTable.read(hdul[4]) + gh_cuts = QTable.read(hdul["GH_CUTS"]) + valid_energy = QTable.read(hdul["VALID_ENERGY"]) + valid_offset = QTable.read(hdul["VALID_OFFSET"]) spatial_selection_table = None multiplicity_cuts = None - if len(hdul) == 6: - if hdul[5].header["EXTNAME"] == "RAD_MAX": - spatial_selection_table = QTable.read(hdul[5]) - else: - multiplicity_cuts = QTable.read(hdul[5]) - elif len(hdul) == 7: - spatial_selection_table = QTable.read(hdul[5]) - multiplicity_cuts = QTable.read(hdul[6]) + if "RAD_MAX" in hdul: + spatial_selection_table = QTable.read(hdul["RAD_MAX"]) + if "MULTIPLICITY_CUTS" in hdul: + multiplicity_cuts = QTable.read(hdul["MULTIPLICITY_CUTS"]) return cls( quality_query=quality_query, From 6e5ee958253450adc2d7f648be385edfa77bf429 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Fri, 27 Jun 2025 18:41:10 +0200 Subject: [PATCH 06/16] Fix application of multiplicity cut in irf tool --- src/ctapipe/tools/compute_irf.py | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/src/ctapipe/tools/compute_irf.py b/src/ctapipe/tools/compute_irf.py index 6b954537acd..c1991cb7f8c 100644 --- a/src/ctapipe/tools/compute_irf.py +++ b/src/ctapipe/tools/compute_irf.py @@ -356,6 +356,9 @@ def calculate_selections(self, reduced_events: dict) -> dict: operator.ge, ) reduced_events["gammas"]["selected"] = reduced_events["gammas"]["selected_gh"] + reduced_events["gammas"]["selected_gh_multiplicity"] = reduced_events["gammas"][ + "selected_gh" + ] if self.spatial_selection_applied: reduced_events["gammas"]["selected_theta"] = evaluate_binned_cut( reduced_events["gammas"]["theta"], @@ -377,6 +380,9 @@ def calculate_selections(self, reduced_events: dict) -> dict: reduced_events["gammas"]["selected"] &= reduced_events["gammas"][ "selected_multiplicity" ] + reduced_events["gammas"]["selected_gh_multiplicity"] &= reduced_events[ + "gammas" + ]["selected_multiplicity"] if self.do_background: backgrounds = ( @@ -439,7 +445,11 @@ def _make_signal_irf_hdus(self, hdus, sim_info): ) ) hdus.append( - self.psf_maker(events=self.signal_events[self.signal_events["selected_gh"]]) + self.psf_maker( + events=self.signal_events[ + self.signal_events["selected_gh_multiplicity"] + ] + ) ) if self.spatial_selection_applied: # TODO: Support fov binning @@ -474,7 +484,9 @@ def _make_benchmark_hdus(self, hdus): ) hdus.append( self.angular_resolution_maker( - events=self.signal_events[self.signal_events["selected_gh"]], + events=self.signal_events[ + self.signal_events["selected_gh_multiplicity"] + ], ) ) if self.do_background: @@ -488,7 +500,7 @@ def _make_benchmark_hdus(self, hdus): self.sensitivity_maker( signal_events=self.signal_events[self.signal_events["selected"]], background_events=self.background_events[ - self.background_events["selected_gh"] + self.background_events["selected"] ], spatial_selection_table=self.opt_result.spatial_selection_table, gamma_spectrum=self.gamma_target_spectrum, @@ -617,7 +629,7 @@ def start(self): if self.do_background: hdus.append( self.background_maker( - self.background_events[self.background_events["selected_gh"]], + self.background_events[self.background_events["selected"]], self.obs_time, ) ) @@ -625,7 +637,7 @@ def start(self): hdus.append( self.effective_area_maker( events=reduced_events["protons"][ - reduced_events["protons"]["selected_gh"] + reduced_events["protons"]["selected"] ], spatial_selection_applied=self.spatial_selection_applied, signal_is_point_like=False, @@ -637,7 +649,7 @@ def start(self): hdus.append( self.effective_area_maker( events=reduced_events["electrons"][ - reduced_events["electrons"]["selected_gh"] + reduced_events["electrons"]["selected"] ], spatial_selection_applied=self.spatial_selection_applied, signal_is_point_like=False, From 503e5bd74eccdb795971a477f9183180bf8ebef8 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Wed, 3 Dec 2025 14:55:49 +0100 Subject: [PATCH 07/16] Remove rebase artifacts --- src/ctapipe/irf/optimize.py | 2 +- src/ctapipe/tools/compute_irf.py | 14 +++++++------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/ctapipe/irf/optimize.py b/src/ctapipe/irf/optimize.py index 436976c3ee2..2455723a7b3 100644 --- a/src/ctapipe/irf/optimize.py +++ b/src/ctapipe/irf/optimize.py @@ -532,7 +532,7 @@ class PointSourceSensitivityOptimizer(PointSourceSensitivityOptimizerBase): def __call__( self, events: dict[str, QTable], - quality_query: EventQualityQuery, + quality_query: DL2EventQualityQuery, clf_prefix: str, ) -> OptimizationResult: self._check_events(events) diff --git a/src/ctapipe/tools/compute_irf.py b/src/ctapipe/tools/compute_irf.py index c1991cb7f8c..986a918061c 100644 --- a/src/ctapipe/tools/compute_irf.py +++ b/src/ctapipe/tools/compute_irf.py @@ -400,13 +400,13 @@ def calculate_selections(self, reduced_events: dict) -> dict: "selected_gh" ] if self.opt_result.multiplicity_cuts is not None: - reduced_events[bkg_type][ - "selected_multiplicity" - ] = evaluate_binned_cut( - reduced_events[bkg_type]["multiplicity"], - reduced_events[bkg_type]["reco_energy"], - self.opt_result.multiplicity_cuts, - operator.ge, + reduced_events[bkg_type]["selected_multiplicity"] = ( + evaluate_binned_cut( + reduced_events[bkg_type]["multiplicity"], + reduced_events[bkg_type]["reco_energy"], + self.opt_result.multiplicity_cuts, + operator.ge, + ) ) reduced_events[bkg_type]["selected"] &= reduced_events[bkg_type][ "selected_multiplicity" From 41d01f453dbf3f767becb30de9e56916f2625c49 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Wed, 3 Dec 2025 15:14:32 +0100 Subject: [PATCH 08/16] Fix test after rebase --- src/ctapipe/io/tests/test_preprocessing.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/ctapipe/io/tests/test_preprocessing.py b/src/ctapipe/io/tests/test_preprocessing.py index 3aaf2bf577b..65cfdd5f14f 100644 --- a/src/ctapipe/io/tests/test_preprocessing.py +++ b/src/ctapipe/io/tests/test_preprocessing.py @@ -251,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", From 8d19b620fbf9b78d07275f42009542e463e56caf Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Wed, 3 Dec 2025 16:15:25 +0100 Subject: [PATCH 09/16] Do not shadow builtin --- src/ctapipe/irf/optimize.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/irf/optimize.py b/src/ctapipe/irf/optimize.py index 2455723a7b3..578321a1371 100644 --- a/src/ctapipe/irf/optimize.py +++ b/src/ctapipe/irf/optimize.py @@ -66,20 +66,20 @@ def __init__( self.multiplicity_cuts = multiplicity_cuts def __repr__(self): - repr = f"" ) - return repr + return representation def write(self, output_name: Path | str, overwrite: bool = False) -> None: """Write an ``OptimizationResult`` to a file in FITS format.""" From e660fab5e75ef287ecf62a40c30cf1a147c87e53 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Mon, 5 Jan 2026 17:15:33 +0100 Subject: [PATCH 10/16] Enable allow_none for AstroQuantity with explicit physical_type and allow implicit definition of physical_type via default_value --- src/ctapipe/core/tests/test_traits.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/ctapipe/core/tests/test_traits.py b/src/ctapipe/core/tests/test_traits.py index 64ef4cc207a..428a29730d2 100644 --- a/src/ctapipe/core/tests/test_traits.py +++ b/src/ctapipe/core/tests/test_traits.py @@ -372,7 +372,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" @@ -387,7 +389,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( @@ -423,7 +427,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 From 754be2f29c6d23d1aa0e4b12d890443d9a58ebae Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Mon, 5 Jan 2026 17:18:21 +0100 Subject: [PATCH 11/16] Address comments --- src/ctapipe/irf/optimize.py | 82 +++++++++++------------- src/ctapipe/resources/optimize_cuts.yaml | 6 +- 2 files changed, 41 insertions(+), 47 deletions(-) diff --git a/src/ctapipe/irf/optimize.py b/src/ctapipe/irf/optimize.py index 578321a1371..282eccf7e6a 100644 --- a/src/ctapipe/irf/optimize.py +++ b/src/ctapipe/irf/optimize.py @@ -224,9 +224,9 @@ def __call__(self, gammaness, reco_energy, reco_energy_bins): self.smoothing = None return calculate_percentile_cut( - gammaness, - reco_energy, - reco_energy_bins, + values=gammaness, + bin_values=reco_energy, + bins=reco_energy_bins, smoothing=self.smoothing, percentile=100 - self.target_percentile, fill_value=gammaness.max(), @@ -238,15 +238,17 @@ class ThetaPercentileCutCalculator(Component): """Computes a percentile cut on theta.""" theta_min_angle = AstroQuantity( - default_value=u.Quantity(-1, u.deg), + default_value=None, + allow_none=True, physical_type=u.physical.angle, - help="Smallest angular cut value allowed (-1 means no cut)", + help="Smallest angular cut value allowed (None means no cut)", ).tag(config=True) theta_max_angle = AstroQuantity( - default_value=u.Quantity(0.32, u.deg), + default_value=u.Quantity(0.5, u.deg), + allow_none=True, physical_type=u.physical.angle, - help="Largest angular cut value allowed", + help="Largest angular cut value allowed (None means no cut)", ).tag(config=True) min_counts = Integer( @@ -255,7 +257,7 @@ class ThetaPercentileCutCalculator(Component): ).tag(config=True) theta_fill_value = AstroQuantity( - default_value=u.Quantity(0.32, u.deg), + default_value=u.Quantity(0.5, u.deg), physical_type=u.physical.angle, help="Angular cut value used for bins with too few events", ).tag(config=True) @@ -272,25 +274,15 @@ class ThetaPercentileCutCalculator(Component): ).tag(config=True) def __call__(self, theta, reco_energy, reco_energy_bins): - if self.theta_min_angle < 0 * u.deg: - theta_min_angle = None - else: - theta_min_angle = self.theta_min_angle - - if self.theta_max_angle < 0 * u.deg: - theta_max_angle = None - else: - theta_max_angle = self.theta_max_angle - if self.smoothing and self.smoothing < 0: self.smoothing = None return calculate_percentile_cut( - theta, - reco_energy, - reco_energy_bins, - min_value=theta_min_angle, - max_value=theta_max_angle, + values=theta, + bin_values=reco_energy, + bins=reco_energy_bins, + min_value=self.theta_min_angle, + max_value=self.theta_max_angle, smoothing=self.smoothing, percentile=self.target_percentile, fill_value=self.theta_fill_value, @@ -325,9 +317,9 @@ def __call__( self.reco_energy_bins, ) gh_mask = evaluate_binned_cut( - events["signal"]["gh_score"], - events["signal"]["reco_energy"], - gh_cuts, + values=events["signal"]["gh_score"], + bin_values=events["signal"]["reco_energy"], + cut_table=gh_cuts, op=operator.ge, ) spatial_selection_table = self.theta_cut_calculator( @@ -428,8 +420,8 @@ def __call__( self._check_events(events) initial_gh_cuts = calculate_percentile_cut( - events["signal"]["gh_score"], - events["signal"]["reco_energy"], + values=events["signal"]["gh_score"], + bin_values=events["signal"]["reco_energy"], bins=self.reco_energy_bins, fill_value=0.0, percentile=100 * (1 - self.initial_gh_cut_efficency), @@ -437,9 +429,9 @@ def __call__( smoothing=1, ) initial_gh_mask = evaluate_binned_cut( - events["signal"]["gh_score"], - events["signal"]["reco_energy"], - initial_gh_cuts, + values=events["signal"]["gh_score"], + bin_values=events["signal"]["reco_energy"], + cut_table=initial_gh_cuts, op=operator.gt, ) @@ -456,8 +448,8 @@ def __call__( self.gh_cut_efficiency_step, ) opt_sensitivity, gh_cuts = optimize_gh_cut( - events["signal"], - events["background"], + signal=events["signal"], + background=events["background"], reco_energy_bins=self.reco_energy_bins, gh_cut_efficiencies=gh_cut_efficiencies, op=operator.ge, @@ -470,10 +462,10 @@ def __call__( # Re-calculate theta cut with optimized g/h cut gh_mask = evaluate_binned_cut( - events["signal"]["gh_score"], - events["signal"]["reco_energy"], - gh_cuts, - operator.ge, + values=events["signal"]["gh_score"], + bin_values=events["signal"]["reco_energy"], + cut_table=gh_cuts, + op=operator.ge, ) events["signal"]["selected_gh"] = gh_mask spatial_selection_table_opt = self.theta_cut_calculator( @@ -499,19 +491,21 @@ def __call__( class PointSourceSensitivityOptimizer(PointSourceSensitivityOptimizerBase): """ Finds the combination of G/H cut, theta cut, and cut on the event-multiplicity - for each energy bin with the maximum point source sensitivity in a grid search. + for each energy bin with the best point source sensitivity in a grid search. """ theta_min_angle = AstroQuantity( - default_value=u.Quantity(-1, u.deg), + default_value=None, + allow_none=True, physical_type=u.physical.angle, - help="Smallest angular cut value allowed (-1 means no cut)", + help="Smallest angular cut value allowed (None means no cut)", ).tag(config=True) theta_max_angle = AstroQuantity( - default_value=u.Quantity(0.32, u.deg), + default_value=u.Quantity(0.5, u.deg), + allow_none=True, physical_type=u.physical.angle, - help="Largest angular cut value allowed", + help="Largest angular cut value allowed (None means no cut)", ).tag(config=True) max_theta_cut_efficiency = Float( @@ -549,8 +543,8 @@ def __call__( ) opt_sensitivity, multiplicity_cuts, theta_cuts, gh_cuts = optimize_cuts( - events["signal"], - events["background"], + signal=events["signal"], + background=events["background"], reco_energy_bins=self.reco_energy_bins, multiplicity_cuts=self.multiplicity_cuts, gh_cut_efficiencies=gh_cut_efficiencies, diff --git a/src/ctapipe/resources/optimize_cuts.yaml b/src/ctapipe/resources/optimize_cuts.yaml index e99ed0fe1df..91f28834ba9 100644 --- a/src/ctapipe/resources/optimize_cuts.yaml +++ b/src/ctapipe/resources/optimize_cuts.yaml @@ -30,9 +30,9 @@ DefaultRecoEnergyBins: reco_energy_n_bins_per_decade: 5 ThetaPercentileCutCalculator: - theta_min_angle: -1 deg - theta_max_angle: 0.32 deg - theta_fill_value: 0.32 deg + theta_min_angle: None + theta_max_angle: 0.5 deg + theta_fill_value: 0.5 deg smoothing: target_percentile: 68 min_counts: 10 From 3b8e7e566988b437cd436574dd4d0d3cd3835d4a Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Thu, 15 Jan 2026 15:51:16 +0100 Subject: [PATCH 12/16] Update resource configs and docstring --- src/ctapipe/irf/optimize.py | 4 ++-- src/ctapipe/resources/compute_irf.yaml | 1 - src/ctapipe/resources/optimize_cuts.yaml | 26 +++++++++++++++--------- 3 files changed, 18 insertions(+), 13 deletions(-) diff --git a/src/ctapipe/irf/optimize.py b/src/ctapipe/irf/optimize.py index 282eccf7e6a..414891a4da3 100644 --- a/src/ctapipe/irf/optimize.py +++ b/src/ctapipe/irf/optimize.py @@ -293,9 +293,9 @@ def __call__(self, theta, reco_energy, reco_energy_bins): class PercentileCuts(CutOptimizerBase): """ Calculates G/H separation cut based on the percentile of signal events - to keep in each bin. + to keep in each bin via ``GhPercentileCutCalculator``. Optionally also calculates a percentile cut on theta based on the signal - events surviving this G/H cut. + events surviving this G/H cut via ``ThetaPercentileCutCalculator``. """ def __init__(self, config=None, parent=None, **kwargs): diff --git a/src/ctapipe/resources/compute_irf.yaml b/src/ctapipe/resources/compute_irf.yaml index 411e1782440..9fe401d5414 100644 --- a/src/ctapipe/resources/compute_irf.yaml +++ b/src/ctapipe/resources/compute_irf.yaml @@ -25,7 +25,6 @@ DL2EventPreprocessor: DL2EventQualityQuery: quality_criteria: - - ["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"] diff --git a/src/ctapipe/resources/optimize_cuts.yaml b/src/ctapipe/resources/optimize_cuts.yaml index 91f28834ba9..091d7b8184f 100644 --- a/src/ctapipe/resources/optimize_cuts.yaml +++ b/src/ctapipe/resources/optimize_cuts.yaml @@ -10,7 +10,7 @@ EventSelectionOptimizer: proton_target_spectrum: IRFDOC_PROTON_SPECTRUM electron_target_spectrum: IRFDOC_ELECTRON_SPECTRUM obs_time: 50 hour - optimization_algorithm: "PointSourceSensitivityOptimizer" # Alternative: "PercentileCuts" + optimization_algorithm: "PointSourceSensitivityOptimizer" # Alternatives: "PointSourceSensitivityGhOptimizer", "PercentileCuts" DL2EventPreprocessor: energy_reconstructor: "RandomForestRegressor" @@ -19,7 +19,6 @@ DL2EventPreprocessor: DL2EventQualityQuery: quality_criteria: - - ["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"] @@ -29,6 +28,21 @@ DefaultRecoEnergyBins: reco_energy_max: 150 TeV reco_energy_n_bins_per_decade: 5 +PointSourceSensitivityOptimizer: + min_background_fov_offset: 0.0 deg + max_background_fov_offset: 5.0 deg + max_gh_cut_efficiency: 0.8 + gh_cut_efficiency_step: 0.1 + alpha: 0.2 + theta_min_angle: None + theta_max_angle: 0.5 deg + max_theta_cut_efficiency: 0.8 + theta_cut_efficiency_step: 0.1 + multiplicity_cuts: [2, 3, 4, 5, 6, 7, 8, 9, 10] + +PointSourceSensitivityGhOptimizer: + initial_gh_cut_efficency: 0.4 + ThetaPercentileCutCalculator: theta_min_angle: None theta_max_angle: 0.5 deg @@ -41,11 +55,3 @@ GhPercentileCutCalculator: target_percentile: 68 min_counts: 10 smoothing: - -PointSourceSensitivityOptimizer: - min_background_fov_offset: 0.0 deg - max_background_fov_offset: 5.0 deg - initial_gh_cut_efficency: 0.4 - max_gh_cut_efficiency: 0.8 - gh_cut_efficiency_step: 0.1 - alpha: 0.2 From 60782bbbb4d2fd9f0c1549f9e3a57d2f546ec28f Mon Sep 17 00:00:00 2001 From: Lukas Beiske <43672561+LukasBeiske@users.noreply.github.com> Date: Mon, 1 Jun 2026 14:41:53 +0200 Subject: [PATCH 13/16] None -> null in optimize_cuts.yaml Co-authored-by: Jonas Hackfeld <53918415+Hckjs@users.noreply.github.com> --- src/ctapipe/resources/optimize_cuts.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe/resources/optimize_cuts.yaml b/src/ctapipe/resources/optimize_cuts.yaml index 091d7b8184f..5839f0b7167 100644 --- a/src/ctapipe/resources/optimize_cuts.yaml +++ b/src/ctapipe/resources/optimize_cuts.yaml @@ -34,7 +34,7 @@ PointSourceSensitivityOptimizer: max_gh_cut_efficiency: 0.8 gh_cut_efficiency_step: 0.1 alpha: 0.2 - theta_min_angle: None + theta_min_angle: null theta_max_angle: 0.5 deg max_theta_cut_efficiency: 0.8 theta_cut_efficiency_step: 0.1 From 6829da715dc251dd4b58a543725d0b9397d0461c Mon Sep 17 00:00:00 2001 From: Lukas Beiske <43672561+LukasBeiske@users.noreply.github.com> Date: Mon, 1 Jun 2026 14:42:13 +0200 Subject: [PATCH 14/16] None -> null in optimize_cuts.yaml Co-authored-by: Jonas Hackfeld <53918415+Hckjs@users.noreply.github.com> --- src/ctapipe/resources/optimize_cuts.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe/resources/optimize_cuts.yaml b/src/ctapipe/resources/optimize_cuts.yaml index 5839f0b7167..26593abc832 100644 --- a/src/ctapipe/resources/optimize_cuts.yaml +++ b/src/ctapipe/resources/optimize_cuts.yaml @@ -44,7 +44,7 @@ PointSourceSensitivityGhOptimizer: initial_gh_cut_efficency: 0.4 ThetaPercentileCutCalculator: - theta_min_angle: None + theta_min_angle: null theta_max_angle: 0.5 deg theta_fill_value: 0.5 deg smoothing: From 2fe1e5f2b8b96ec7307bf73504c1c0f46b11dad6 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Mon, 1 Jun 2026 17:44:42 +0200 Subject: [PATCH 15/16] Address comments --- src/ctapipe/io/dl2_tables_preprocessing.py | 22 ++++++++-------------- 1 file changed, 8 insertions(+), 14 deletions(-) diff --git a/src/ctapipe/io/dl2_tables_preprocessing.py b/src/ctapipe/io/dl2_tables_preprocessing.py index 715ede92aed..e12cac1f28b 100644 --- a/src/ctapipe/io/dl2_tables_preprocessing.py +++ b/src/ctapipe/io/dl2_tables_preprocessing.py @@ -441,27 +441,21 @@ def make_derived_columns(self, events: QTable) -> QTable: 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 + reconstructors = ( + self.epp.energy_reconstructor, + self.epp.geometry_reconstructor, + self.epp.gammaness_classifier, + ) 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 - ), + np.count_nonzero(events[f"{reconstructor}_telescopes"], axis=1) + for reconstructor in reconstructors ], 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", - ] + [f"{reconstructor}_telescopes" for reconstructor in reconstructors] ) return events From aacf57827b164d4cc86374f7034790c70d9be429 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Fri, 12 Jun 2026 13:34:36 +0200 Subject: [PATCH 16/16] Improve docstrings --- src/ctapipe/irf/optimize.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/irf/optimize.py b/src/ctapipe/irf/optimize.py index 414891a4da3..bf74b8b7030 100644 --- a/src/ctapipe/irf/optimize.py +++ b/src/ctapipe/irf/optimize.py @@ -398,9 +398,11 @@ def _get_valid_energy_range(self, opt_sensitivity): class PointSourceSensitivityGhOptimizer(PointSourceSensitivityOptimizerBase): """ - EventDisplay-like procedure of optimizing a G/H cut for maximum point source - sensitivity, calculating only a percentile cut on theta, and no cut on - multiplicity. + Optimize gamma-hadron-separation cuts for best point-source sensitivity. + + This approach computes a simple, percentile based cut for theta and + then performs a grid search for best sensitivity only for the G/H cut + using `pyirf.cut_optimization.optimize_gh_cut`. """ initial_gh_cut_efficency = Float( @@ -490,8 +492,11 @@ def __call__( class PointSourceSensitivityOptimizer(PointSourceSensitivityOptimizerBase): """ - Finds the combination of G/H cut, theta cut, and cut on the event-multiplicity - for each energy bin with the best point source sensitivity in a grid search. + Optimize cuts for best point-source sensitivity. + + This optimizer employes a brute-force grid search over gammaness, + theta and event-multiplicity cuts in each bin of reconstructed energy + using `pyirf.cut_optimization.optimize_gh_cut`. """ theta_min_angle = AstroQuantity(