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``. 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/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 diff --git a/src/ctapipe/io/dl2_tables_preprocessing.py b/src/ctapipe/io/dl2_tables_preprocessing.py index c0776496d75..e12cac1f28b 100644 --- a/src/ctapipe/io/dl2_tables_preprocessing.py +++ b/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,23 @@ 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 + reconstructors = ( + self.epp.energy_reconstructor, + self.epp.geometry_reconstructor, + self.epp.gammaness_classifier, + ) + events["multiplicity"] = np.min( + [ + 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"{reconstructor}_telescopes" for reconstructor in reconstructors] + ) 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..65cfdd5f14f 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, @@ -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", 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..bf74b8b7030 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): + representation = f"" - ) - else: - return ( - f"" - ) + representation += f", and {len(self.spatial_selection_table)} theta bins" + + if self.multiplicity_cuts is not None: + representation += f", and {len(self.multiplicity_cuts)} multiplicity bins" + + representation += ( + 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 representation 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: @@ -120,16 +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]) - spatial_selection_table = QTable.read(hdul[5]) if len(hdul) > 5 else None + 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 "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, @@ -140,6 +152,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 +161,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' " @@ -213,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(), @@ -227,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( @@ -244,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) @@ -261,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, @@ -290,13 +293,11 @@ 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``. """ - 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) @@ -316,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( @@ -341,18 +342,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 +383,32 @@ 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): + """ + 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( + 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) @@ -399,8 +422,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), @@ -408,9 +431,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, ) @@ -427,8 +450,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, @@ -441,10 +464,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( @@ -466,14 +489,89 @@ 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): + """ + 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( + default_value=None, + allow_none=True, + physical_type=u.physical.angle, + help="Smallest angular cut value allowed (None means no cut)", + ).tag(config=True) + + theta_max_angle = AstroQuantity( + default_value=u.Quantity(0.5, u.deg), + allow_none=True, + physical_type=u.physical.angle, + help="Largest angular cut value allowed (None means no cut)", + ).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, 3, 4, 5, 6, 7, 8, 9, 10], + help="Event-multiplicity cuts used for scanning after optimal cut", + ).tag(config=True) + + def __call__( + self, + events: dict[str, QTable], + quality_query: DL2EventQualityQuery, + 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( + signal=events["signal"], + background=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/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"]) 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 e99ed0fe1df..26593abc832 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,10 +28,25 @@ 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: null + 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: -1 deg - theta_max_angle: 0.32 deg - theta_fill_value: 0.32 deg + theta_min_angle: null + theta_max_angle: 0.5 deg + theta_fill_value: 0.5 deg smoothing: target_percentile: 68 min_counts: 10 @@ -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 diff --git a/src/ctapipe/tools/compute_irf.py b/src/ctapipe/tools/compute_irf.py index dd47f8755c5..986a918061c 100644 --- a/src/ctapipe/tools/compute_irf.py +++ b/src/ctapipe/tools/compute_irf.py @@ -355,6 +355,10 @@ def calculate_selections(self, reduced_events: dict) -> dict: self.opt_result.gh_cuts, 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"], @@ -362,14 +366,23 @@ 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" ] + reduced_events["gammas"]["selected_gh_multiplicity"] &= reduced_events[ + "gammas" + ]["selected_multiplicity"] if self.do_background: backgrounds = ( @@ -383,9 +396,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" @@ -418,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 @@ -453,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: @@ -467,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, @@ -596,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, ) ) @@ -604,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, @@ -616,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,