From 29b3f196d9a4d1a3069322963657fe4411c1b4df Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Thu, 7 Sep 2023 14:54:48 +0200 Subject: [PATCH] Added specific configuration of the optimisation grid to the cut optimiser component --- ctapipe/irf/__init__.py | 9 ++- ctapipe/irf/irf_classes.py | 109 ++++++++++++++++++++++++++----------- ctapipe/tools/info.py | 1 + ctapipe/tools/make_irf.py | 31 ++--------- 4 files changed, 90 insertions(+), 60 deletions(-) diff --git a/ctapipe/irf/__init__.py b/ctapipe/irf/__init__.py index b1056dbb807..6cd70eb6a8c 100644 --- a/ctapipe/irf/__init__.py +++ b/ctapipe/irf/__init__.py @@ -1,3 +1,8 @@ -from .irf_classes import CutOptimising, DataBinning, EnergyBinning, EventPreProcessor +from .irf_classes import ( + CutOptimising, + DataBinning, + EventPreProcessor, + OutputEnergyBinning, +) -__all__ = ["CutOptimising", "DataBinning", "EnergyBinning", "EventPreProcessor"] +__all__ = ["CutOptimising", "DataBinning", "OutputEnergyBinning", "EventPreProcessor"] diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index 6c4658d0658..88a1e470e61 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -8,29 +8,72 @@ from astropy.table import QTable from pyirf.binning import create_bins_per_decade from pyirf.cut_optimization import optimize_gh_cut -from pyirf.cuts import calculate_percentile_cut +from pyirf.cuts import calculate_percentile_cut, evaluate_binned_cut from ..core import Component, QualityQuery from ..core.traits import Float, Integer, List, Unicode class CutOptimising(Component): - """Collects settings related to the cut configuration""" + """Performs cut optimisation""" max_gh_cut_efficiency = Float( default_value=0.8, help="Maximum gamma efficiency requested" ).tag(config=True) + gh_cut_efficiency_step = Float( default_value=0.1, help="Stepsize used for scanning after optimal gammaness cut", ).tag(config=True) + initial_gh_cut_efficency = Float( default_value=0.4, help="Start value of gamma efficiency before optimisation" ).tag(config=True) - def optimise_gh_cut( - self, signal, background, Et_bins, Er_bins, bins, alpha, max_bg_radius - ): + reco_energy_min = Float( + help="Minimum value for Reco Energy bins in TeV units", + default_value=0.005, + ).tag(config=True) + + reco_energy_max = Float( + help="Maximum value for Reco Energy bins in TeV units", + default_value=200, + ).tag(config=True) + + reco_energy_n_bins_per_decade = Float( + help="Number of edges per decade for Reco Energy bins", + default_value=5, + ).tag(config=True) + + theta_min_angle = Float( + default_value=0.05, help="Smallest angular cut value allowed" + ).tag(config=True) + + theta_max_angle = Float( + default_value=0.32, help="Largest angular cut value allowed" + ).tag(config=True) + + theta_min_counts = Integer( + default_value=10, + help="Minimum number of events in a bin to attempt to find a cut value", + ).tag(config=True) + + theta_fill_value = Float( + default_value=0.32, help="Angular cut value used for bins with too few events" + ).tag(config=True) + + def reco_energy_bins(self): + """ + Creates bins per decade for reconstructed MC energy using pyirf function. + """ + reco_energy = create_bins_per_decade( + self.reco_energy_min * u.TeV, + self.reco_energy_max * u.TeV, + self.reco_energy_n_bins_per_decade, + ) + return reco_energy + + def optimise_gh_cut(self, signal, background, alpha, max_bg_radius): INITIAL_GH_CUT = np.quantile( signal["gh_score"], (1 - self.initial_gh_cut_efficency) ) @@ -43,11 +86,11 @@ def optimise_gh_cut( theta_cuts = calculate_percentile_cut( signal["theta"][mask_theta_cuts], signal["reco_energy"][mask_theta_cuts], - bins=Et_bins, - min_value=bins.theta_min_angle * u.deg, - max_value=bins.theta_max_angle * u.deg, - fill_value=bins.theta_fill_value * u.deg, - min_events=bins.theta_min_counts, + bins=self.reco_energy_bins(), + min_value=self.theta_min_angle * u.deg, + max_value=self.theta_max_angle * u.deg, + fill_value=self.theta_fill_value * u.deg, + min_events=self.theta_min_counts, percentile=68, ) @@ -61,14 +104,33 @@ def optimise_gh_cut( sens2, gh_cuts = optimize_gh_cut( signal, background, - reco_energy_bins=Er_bins, + reco_energy_bins=self.reco_energy_bins(), gh_cut_efficiencies=gh_cut_efficiencies, op=operator.ge, theta_cuts=theta_cuts, alpha=alpha, fov_offset_max=max_bg_radius * u.deg, ) - return gh_cuts, sens2 + + # now that we have the optimized gh cuts, we recalculate the theta + # cut as 68 percent containment on the events surviving these cuts. + self.log.info("Recalculating theta cut for optimized GH Cuts") + for tab in (signal, background): + tab["selected_gh"] = evaluate_binned_cut( + tab["gh_score"], tab["reco_energy"], gh_cuts, operator.ge + ) + + theta_cuts = calculate_percentile_cut( + signal[signal["selected_gh"]]["theta"], + signal[signal["selected_gh"]]["reco_energy"], + self.reco_energy_bins(), + percentile=68, + min_value=self.theta_min_angle * u.deg, + max_value=self.theta_max_angle * u.deg, + fill_value=self.theta_fill_value * u.deg, + min_events=self.theta_min_counts, + ) + return gh_cuts, theta_cuts, sens2 class EventPreProcessor(QualityQuery): @@ -166,7 +228,7 @@ def make_empty_table(self): return QTable(names=columns, units=units) -class EnergyBinning(Component): +class OutputEnergyBinning(Component): """Collects energy binning settings""" true_energy_min = Float( @@ -186,12 +248,12 @@ class EnergyBinning(Component): reco_energy_min = Float( help="Minimum value for Reco Energy bins in TeV units", - default_value=0.005, + default_value=0.006, ).tag(config=True) reco_energy_max = Float( help="Maximum value for Reco Energy bins in TeV units", - default_value=200, + default_value=190, ).tag(config=True) reco_energy_n_bins_per_decade = Float( @@ -256,23 +318,6 @@ class DataBinning(Component): Stolen from LSTChain """ - theta_min_angle = Float( - default_value=0.05, help="Smallest angular cut value allowed" - ).tag(config=True) - - theta_max_angle = Float( - default_value=0.32, help="Largest angular cut value allowed" - ).tag(config=True) - - theta_min_counts = Integer( - default_value=10, - help="Minimum number of events in a bin to attempt to find a cut value", - ).tag(config=True) - - theta_fill_value = Float( - default_value=0.32, help="Angular cut value used for bins with too few events" - ).tag(config=True) - fov_offset_min = Float( help="Minimum value for FoV Offset bins in degrees", default_value=0.0, diff --git a/ctapipe/tools/info.py b/ctapipe/tools/info.py index 66b28ba97c8..4365dc4481c 100644 --- a/ctapipe/tools/info.py +++ b/ctapipe/tools/info.py @@ -30,6 +30,7 @@ "iminuit", "tables", "eventio", + "pyirf", ] ) diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index cbf4ae84ce4..f9ad3871ca6 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -8,7 +8,7 @@ from astropy.table import vstack from pyirf.benchmarks import angular_resolution, energy_bias_resolution from pyirf.binning import create_histogram_table -from pyirf.cuts import calculate_percentile_cut, evaluate_binned_cut +from pyirf.cuts import evaluate_binned_cut from pyirf.io import ( create_aeff2d_hdu, create_background_2d_hdu, @@ -36,7 +36,7 @@ from ..core import Provenance, Tool, traits from ..core.traits import Bool, Float, Integer, Unicode from ..io import TableLoader -from ..irf import CutOptimising, DataBinning, EnergyBinning, EventPreProcessor +from ..irf import CutOptimising, DataBinning, EventPreProcessor, OutputEnergyBinning class Spectra(Enum): @@ -115,7 +115,7 @@ class IrfTool(Tool): default_value=3.0, help="Radius used to calculate background rate in degrees" ).tag(config=True) - classes = [CutOptimising, DataBinning, EnergyBinning, EventPreProcessor] + classes = [CutOptimising, DataBinning, OutputEnergyBinning, EventPreProcessor] def make_derived_columns(self, kind, events, spectrum, target_spectrum, obs_conf): @@ -213,7 +213,7 @@ def load_preselected_events(self): def setup(self): self.co = CutOptimising(parent=self) - self.e_bins = EnergyBinning(parent=self) + self.e_bins = OutputEnergyBinning(parent=self) self.bins = DataBinning(parent=self) self.epp = EventPreProcessor(parent=self) @@ -227,34 +227,13 @@ def setup(self): def start(self): self.load_preselected_events() - self.gh_cuts, sens2 = self.co.optimise_gh_cut( + self.gh_cuts, self.theta_cuts_opt, sens2 = self.co.optimise_gh_cut( self.signal, self.background, - self.true_energy_bins, - self.reco_energy_bins, - self.bins, self.alpha, self.max_bg_radius, ) - # now that we have the optimized gh cuts, we recalculate the theta - # cut as 68 percent containment on the events surviving these cuts. - self.log.info("Recalculating theta cut for optimized GH Cuts") - for tab in (self.signal, self.background): - tab["selected_gh"] = evaluate_binned_cut( - tab["gh_score"], tab["reco_energy"], self.gh_cuts, operator.ge - ) - - self.theta_cuts_opt = calculate_percentile_cut( - self.signal[self.signal["selected_gh"]]["theta"], - self.signal[self.signal["selected_gh"]]["reco_energy"], - self.true_energy_bins, - percentile=68, - min_value=self.bins.theta_min_angle * u.deg, - max_value=self.bins.theta_max_angle * u.deg, - fill_value=self.bins.theta_fill_value * u.deg, - min_events=self.bins.theta_min_counts, - ) self.signal["selected_theta"] = evaluate_binned_cut( self.signal["theta"], self.signal["reco_energy"],