From 26c0f63f0786da57e4fb8fd4608575daff8f6036 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Tue, 18 Apr 2023 17:28:39 +0200 Subject: [PATCH 01/37] First kludge of protopipe irf code as a ctapipe tool --- ctapipe/irf/__init__.py | 3 + ctapipe/irf/irf_classes.py | 404 +++++++++++++++++++++++++++++++++++++ ctapipe/tools/make-irf.py | 342 +++++++++++++++++++++++++++++++ setup.cfg | 2 + 4 files changed, 751 insertions(+) create mode 100644 ctapipe/irf/__init__.py create mode 100644 ctapipe/irf/irf_classes.py create mode 100644 ctapipe/tools/make-irf.py diff --git a/ctapipe/irf/__init__.py b/ctapipe/irf/__init__.py new file mode 100644 index 00000000000..362baef0e80 --- /dev/null +++ b/ctapipe/irf/__init__.py @@ -0,0 +1,3 @@ +from .irf_classes import DataBinning, IrfToolBase + +__all__ = ["IrfToolBase", "DataBinning"] diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py new file mode 100644 index 00000000000..ce572f7c6f1 --- /dev/null +++ b/ctapipe/irf/irf_classes.py @@ -0,0 +1,404 @@ +""" +Define a parent IrfTool class to hold all the options +""" +import astropy.units as u +import numpy as np +from astropy.table import QTable +from pyirf.binning import create_bins_per_decade + +from ..core import Component, QualityQuery, Tool, traits +from ..core.traits import Bool, Float, Integer, List, Unicode + + +class IrfToolBase(Tool): + + gamma_file = traits.Path( + default_value=None, directory_ok=False, help="Gamma input filename and path" + ).tag(config=True) + gamma_sim_spectrum = traits.Unicode( + default_value="CRAB_HEGRA", + help="Name of the pyrif spectra used for the simulated gamma spectrum", + ).tag(config=True) + proton_file = traits.Path( + default_value=None, directory_ok=False, help="Gamma input filename and path" + ).tag(config=True) + proton_sim_spectrum = traits.Unicode( + default_value="IRFDOC_PROTON_SPECTRUM", + help="Name of the pyrif spectra used for the simulated proton spectrum", + ).tag(config=True) + electron_file = traits.Path( + default_value=None, directory_ok=False, help="Gamma input filename and path" + ).tag(config=True) + electron_sim_spectrum = traits.Unicode( + default_value="IRFDOC_ELECTRON_SPECTRUM", + help="Name of the pyrif spectra used for the simulated electron spectrum", + ).tag(config=True) + + chunk_size = Integer( + default_value=100000, + allow_none=True, + help="How many subarray events to load at once for making predictions.", + ).tag(config=True) + + output_path = traits.Path( + default_value=None, + allow_none=False, + directory_ok=False, + help="Output file", + ).tag(config=True) + output_file = Unicode( + default_value=None, allow_none=False, help="Name for the output file" + ).tag(config=True) + overwrite = Bool( + False, + help="Overwrite the output file if it exists", + ).tag(config=True) + + obs_time = Float(default_value=50.0, help="Observation time").tag(config=True) + obs_time_unit = Unicode( + default_value="hour", + help="Unit used to specify observation time as an astropy unit string.", + ).tag(config=True) + + alpha = Float( + default_value=5.0, help="Ratio between size of on and off regions" + ).tag(config=True) + + max_bg_radius = Float( + default_value=5.0, help="Radius used to calculate background rate in degrees" + ).tag(config=True) + + max_gh_cut_efficiency = Float( + default_value=0.8, help="Maximum gamma purity requested" + ).tag(config=True) + gh_cut_efficiency_step = Float( + default_value=0.01, + 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 purity before optimisatoin" + ).tag(config=True) + + energy_reconstructor = Unicode( + default_value="RandomForestRegressor", + help="Prefix of the reco `_energy` column", + ).tag(config=True) + geometry_reconstructor = Unicode( + default_value="HillasReconstructor", + help="Prefix of the `_alt` and `_az` reco geometry columns", + ).tag(config=True) + gammaness_classifier = Unicode( + default_value="RandomForestClassifier", + help="Prefix of the classifier `_prediction` column", + ).tag(config=True) + + preselect_criteria = List( + default_value=[ + ("multiplicity 4", "np.count_nonzero(tels,axis=1) >= 4"), + ("valid classifier", "valid_classer"), + ("valid geom reco", "valid_geom"), + ("valid energy reco", "valid_energy"), + ], + help=QualityQuery.quality_criteria.help, + ).tag(config=True) + + rename_columns = List( + help="List containing translation pairs of quality columns" + "used for quality filters and their names as given in the input file used." + "Ex: [('valid_geom','HillasReconstructor_is_valid')]", + default_value=[ + ("valid_geom", "HillasReconstructor_is_valid"), + ("valid_energy", "RandomForestRegressor_is_valid"), + ("valid_classer", "RandomForestClassifier_is_valid"), + ], + ) + + def _preselect_events(self, events): + keep_columns = [ + "obs_id", + "event_id", + "true_energy", + "true_az", + "true_alt", + ] + rename_from = [ + f"{self.energy_reconstructor}_energy", + f"{self.geometry_reconstructor}_az", + f"{self.geometry_reconstructor}_alt", + f"{self.gammaness_classifier}_prediction", + ] + rename_to = ["reco_energy", "reco_az", "reco_alt", "gh_score"] + + for new, old in self.rename_columns: + rename_from.append(old) + rename_to.append(new) + + keep_columns.append(rename_from) + events = QTable(events[keep_columns], copy=False) + events.rename_columns(rename_from, rename_to) + keep = QualityQuery(quality_criteria=self.preselect_criteria).get_table_mask( + events + ) + + return events[keep] + + def _make_empty_table(self): + columns = [ + "obs_id", + "event_id", + "true_energy", + "true_az", + "true_alt", + "reco_energy", + "reco_az", + "reco_alt", + "gh_score", + "pointing_az", + "pointing_alt", + "true_source_fov_offset", + "reco_source_fov_offset", + "weights", + ] + units = [ + None, + None, + u.TeV, + u.deg, + u.deg, + u.TeV, + u.deg, + u.deg, + None, + u.deg, + u.deg, + u.deg, + u.deg, + None, + ] + + return QTable(names=columns, units=units) + + +class ThetaSettings(Component): + + min_angle = Float( + default_value=0.05, help="Smallest angular cut value allowed" + ).tag(config=True) + max_angle = Float(default_value=0.32, help="Largest angular cut value allowed").tag( + config=True + ) + min_counts = Integer( + default_value=10, + help="Minimum number of events in a bin to attempt to find a cut value", + ).tag(config=True) + fill_value = Float( + default_value=0.32, help="Angular cut value used for bins with too few events" + ).tag(config=True) + + +class DataBinning(Component): + """ + Collects information on generating energy and angular bins for + generating IRFs as per pyIRF requirements. + + Stolen from LSTChain + """ + + true_energy_min = Float( + help="Minimum value for True Energy bins in TeV units", + default_value=0.005, + ).tag(config=True) + + true_energy_max = Float( + help="Maximum value for True Energy bins in TeV units", + default_value=200, + ).tag(config=True) + + true_energy_n_bins_per_decade = Float( + help="Number of edges per decade for True Energy bins", + default_value=5, + ).tag(config=True) + + 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) + + energy_migration_min = Float( + help="Minimum value of Energy Migration matrix", + default_value=0.2, + ).tag(config=True) + + energy_migration_max = Float( + help="Maximum value of Energy Migration matrix", + default_value=5, + ).tag(config=True) + + energy_migration_n_bins = Integer( + help="Number of bins in log scale for Energy Migration matrix", + default_value=31, + ).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) + + fov_offset_min = Float( + help="Minimum value for FoV Offset bins", + default_value=0.1, + ).tag(config=True) + + fov_offset_max = Float( + help="Maximum value for FoV offset bins", + default_value=1.1, + ).tag(config=True) + + fov_offset_n_edges = Integer( + help="Number of edges for FoV offset bins", + default_value=9, + ).tag(config=True) + + bkg_fov_offset_min = Float( + help="Minimum value for FoV offset bins for Background IRF", + default_value=0, + ).tag(config=True) + + bkg_fov_offset_max = Float( + help="Maximum value for FoV offset bins for Background IRF", + default_value=10, + ).tag(config=True) + + bkg_fov_offset_n_edges = Integer( + help="Number of edges for FoV offset bins for Background IRF", + default_value=21, + ).tag(config=True) + + source_offset_min = Float( + help="Minimum value for Source offset for PSF IRF", + default_value=0, + ).tag(config=True) + + source_offset_max = Float( + help="Maximum value for Source offset for PSF IRF", + default_value=1, + ).tag(config=True) + + source_offset_n_edges = Integer( + help="Number of edges for Source offset for PSF IRF", + default_value=101, + ).tag(config=True) + + def true_energy_bins(self): + """ + Creates bins per decade for true MC energy using pyirf function. + The overflow binning added is not needed at the current stage. + + Examples + -------- + It can be used as: + + >>> add_overflow_bins(***)[1:-1] + """ + true_energy = create_bins_per_decade( + self.true_energy_min * u.TeV, + self.true_energy_max * u.TeV, + self.true_energy_n_bins_per_decade, + ) + return true_energy + + def reco_energy_bins(self): + """ + Creates bins per decade for reconstructed MC energy using pyirf function. + The overflow binning added is not needed at the current stage. + + Examples + -------- + It can be used as: + + >>> add_overflow_bins(***)[1:-1] + """ + 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 energy_migration_bins(self): + """ + Creates bins for energy migration. + """ + energy_migration = np.geomspace( + self.energy_migration_min, + self.energy_migration_max, + self.energy_migration_n_bins, + ) + return energy_migration + + def fov_offset_bins(self): + """ + Creates bins for single/multiple FoV offset. + """ + fov_offset = ( + np.linspace( + self.fov_offset_min, + self.fov_offset_max, + self.fov_offset_n_edges, + ) + * u.deg + ) + return fov_offset + + def bkg_fov_offset_bins(self): + """ + Creates bins for FoV offset for Background IRF, + Using the same binning as in pyirf example. + """ + background_offset = ( + np.linspace( + self.bkg_fov_offset_min, + self.bkg_fov_offset_max, + self.bkg_fov_offset_n_edges, + ) + * u.deg + ) + return background_offset + + def source_offset_bins(self): + """ + Creates bins for source offset for generating PSF IRF. + Using the same binning as in pyirf example. + """ + + source_offset = ( + np.linspace( + self.source_offset_min, + self.source_offset_max, + self.source_offset_n_edges, + ) + * u.deg + ) + return source_offset diff --git a/ctapipe/tools/make-irf.py b/ctapipe/tools/make-irf.py new file mode 100644 index 00000000000..f5608f09c24 --- /dev/null +++ b/ctapipe/tools/make-irf.py @@ -0,0 +1,342 @@ +"""Tool to generate IRFs""" +import operator +from pathlib import Path + +import astropy.units as u +import numpy as np +from astropy.io import fits +from astropy.table import vstack +from pyirf.benchmarks import angular_resolution, energy_bias_resolution +from pyirf.binning import ( + add_overflow_bins, + create_bins_per_decade, + create_histogram_table, +) +from pyirf.cut_optimization import optimize_gh_cut +from pyirf.cuts import calculate_percentile_cut, evaluate_binned_cut +from pyirf.io import ( + create_aeff2d_hdu, + create_background_2d_hdu, + create_energy_dispersion_hdu, + create_psf_table_hdu, + create_rad_max_hdu, +) +from pyirf.irf import ( + background_2d, + effective_area_per_energy, + energy_dispersion, + psf_table, +) +from pyirf.sensitivity import calculate_sensitivity, estimate_background +from pyirf.simulations import SimulatedEventsInfo +from pyirf.spectral import ( + CRAB_HEGRA, + IRFDOC_ELECTRON_SPECTRUM, + IRFDOC_PROTON_SPECTRUM, + PowerLaw, + calculate_event_weights, +) +from pyirf.utils import calculate_source_fov_offset, calculate_theta + +from ..core import Provenance +from ..io import TableLoader +from ..irf import DataBinning, IrfToolBase + +PYIRF_SPECTRA = { + "CRAB_HEGRA": CRAB_HEGRA, + "IRFDOC_ELECTRON_SPECTRUM": IRFDOC_ELECTRON_SPECTRUM, + "IRFDOC_PROTON_SPECTRUM": IRFDOC_PROTON_SPECTRUM, +} + + +class IrfTool(IrfToolBase, DataBinning): + name = "ctapipe-make-irfs" + description = "Tool to create IRF files in GAD format" + + def make_derived_columns(self, events, spectrum, target_spectrum): + events["pointing_az"] = 0 * u.deg + events["pointing_alt"] = 70 * u.deg + + events["theta"] = calculate_theta( + events, + assumed_source_az=events["true_az"], + assumed_source_alt=events["true_alt"], + ) + + events["true_source_fov_offset"] = calculate_source_fov_offset( + events, prefix="true" + ) + events["reco_source_fov_offset"] = calculate_source_fov_offset( + events, prefix="reco" + ) + events["weights"] = calculate_event_weights( + events["true_energy"], + target_spectrum=target_spectrum, + simulated_spectrum=spectrum, + ) + + return events + + def get_sim_info_and_spectrum(self, loader): + sim = loader.read_simulation_configuration() + + sim_info = SimulatedEventsInfo( + n_showers=sum(sim["n_showers"] * sim["shower_reuse"]), + energy_min=sim["energy_range_min"].quantity[0], + energy_max=sim["energy_range_max"].quantity[0], + max_impact=sim["max_scatter_range"].quantity[0], + spectral_index=sim["spectral_index"][0], + viewcone=sim["max_viewcone_radius"].quantity[0], + ) + + return sim_info, PowerLaw.from_simulation( + sim_info, obstime=self.obs_time * u.Unit(self.obs_time_unit) + ) + + def setup(self): + opts = dict(load_dl2=True, load_simulated=True, load_dl1_parameters=False) + reduced_events = dict() + for kind, file, target_spectrum in [ + ("gamma", self.gamma_file, PYIRF_SPECTRA[self.gamma_sim_spectrum]), + ("proton", self.proton_file, PYIRF_SPECTRA[self.proton_sim_spectrum]), + ("electron", self.electron_file, PYIRF_SPECTRA[self.electron_sim_spectrum]), + ]: + with TableLoader(file, **opts) as load: + Provenance().add_input_file(file) + table = self._make_empty_table() + sim_info, spectrum = self.get_sim_info_and_spectrum(load) + if kind == "gamma": + self.sim_info = sim_info + self.spectrum = spectrum + for start, stop, events in load.read_subarray_events_chunked( + self.chunk_size + ): + selected = self._preselect_events(events) + selected = self.make_derived_columns( + selected, spectrum, target_spectrum + ) + table = vstack(table, selected) + + reduced_events[kind] = table + + self.signal = reduced_events["gamma"] + self.background = vstack(reduced_events["proton"], reduced_events["electron"]) + + self.theta_bins = add_overflow_bins( + create_bins_per_decade( + self.sim_info.energy_min, self.sim_info.energy_max, 50 + ) + ) + + self.energy_reco_bins = self.reco_energy_bins() + self.energy_true_bins = self.true_energy_bins() + self.source_offset_bins = self.source_offset_bins() + self.fov_offset_bins = self.fov_offset_bins() + self.energy_migration_bins = self.energy_migration_bins() + + def start(self): + + INITIAL_GH_CUT = np.quantile( + self.signal["gh_score"], (1 - self.initial_gh_cut_efficency) + ) + self.log.info( + f"Using fixed G/H cut of {INITIAL_GH_CUT} to calculate theta cuts" + ) + + mask_theta_cuts = self.signal["gh_score"] >= INITIAL_GH_CUT + + theta_cuts = calculate_percentile_cut( + self.signal["theta"][mask_theta_cuts], + self.signal["reco_energy"][mask_theta_cuts], + bins=self.theta_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, + ) + + self.log.info("Optimizing G/H separation cut for best sensitivity") + 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, + ) + + sens2, self.gh_cuts = optimize_gh_cut( + self.signal, + self.background, + reco_energy_bins=self.energy_reco_bins, + gh_cut_efficiencies=gh_cut_efficiencies, + op=operator.ge, + theta_cuts=theta_cuts, + alpha=self.alpha, + background_radius=self.max_bg_radius * u.deg, + ) + + # 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.theta_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, + ) + self.signal["selected_theta"] = evaluate_binned_cut( + self.signal["theta"], + self.signal["reco_energy"], + self.theta_cuts_opt, + operator.le, + ) + self.signal["selected"] = ( + self.signal["selected_theta"] & self.signal["selected_gh"] + ) + + # calculate sensitivity + signal_hist = create_histogram_table( + self.signal[self.signal["selected"]], bins=self.energy_reco_bins + ) + background_hist = estimate_background( + self.background[self.background["selected_gh"]], + reco_energy_bins=self.energy_reco_bins, + theta_cuts=self.theta_cuts_opt, + alpha=self.alpha, + fov_offset_min=self.fov_offset_min, + fov_offset_max=self.fov_offset_max, + ) + self.sensitivity = calculate_sensitivity( + signal_hist, background_hist, alpha=self.alpha + ) + + # scale relative sensitivity by Crab flux to get the flux sensitivity + for s in (sens2, self.sensitivity): + s["flux_sensitivity"] = s["relative_sensitivity"] * self.spectrum( + s["reco_energy_center"] + ) + + def finalise(self): + + masks = { + "": self.signal["selected"], + "_NO_CUTS": slice(None), + "_ONLY_GH": self.signal["selected_gh"], + "_ONLY_THETA": self.signal["selected_theta"], + } + hdus = [ + fits.PrimaryHDU(), + fits.BinTableHDU(self.sensitivity, name="SENSITIVITY"), + # fits.BinTableHDU(sensitivity_step_2, name="SENSITIVITY_STEP_2"), + # fits.BinTableHDU(self.theta_cuts, name="THETA_CUTS"), + fits.BinTableHDU(self.theta_cuts_opt, name="THETA_CUTS_OPT"), + fits.BinTableHDU(self.gh_cuts, name="GH_CUTS"), + ] + + for label, mask in masks.items(): + effective_area = effective_area_per_energy( + self.signal[mask], + self.sim_info, + true_energy_bins=self.true_energy_bins, + ) + hdus.append( + create_aeff2d_hdu( + effective_area[..., np.newaxis], # +1 dimension for FOV offset + self.true_energy_bins, + self.fov_offset_bins, + extname="EFFECTIVE AREA" + label, + ) + ) + edisp = energy_dispersion( + self.signal[mask], + true_energy_bins=self.true_energy_bins, + fov_offset_bins=self.fov_offset_bins, + migration_bins=self.energy_migration_bins, + ) + hdus.append( + create_energy_dispersion_hdu( + edisp, + true_energy_bins=self.true_energy_bins, + migration_bins=self.energy_migration_bins, + fov_offset_bins=self.fov_offset_bins, + extname="ENERGY_DISPERSION" + label, + ) + ) + # Here we use reconstructed energy instead of true energy for the sake of + # current pipelines comparisons + bias_resolution = energy_bias_resolution( + self.signal[self.signal["selected"]], + self.reco_energy_bins, + energy_type="reco", + ) + + # Here we use reconstructed energy instead of true energy for the sake of + # current pipelines comparisons + ang_res = angular_resolution( + self.signal[self.signal["selected_gh"]], + self.reco_energy_bins, + energy_type="reco", + ) + + psf = psf_table( + self.signal[self.signal["selected_gh"]], + self.true_energy_bins, + fov_offset_bins=self.fov_offset_bins, + source_offset_bins=self.source_offset_bins, + ) + + background_rate = background_2d( + self.background[self.background["selected_gh"]], + self.reco_energy_bins, + fov_offset_bins=np.arange(0, 11) * u.deg, + t_obs=self.obs_time * u.Unit(self.obs_time_unit), + ) + + hdus.append( + create_background_2d_hdu( + background_rate, + self.reco_energy_bins, + fov_offset_bins=np.arange(0, 11) * u.deg, + ) + ) + + hdus.append( + create_psf_table_hdu( + psf, + self.true_energy_bins, + self.source_offset_bins, + self.fov_offset_bins, + ) + ) + hdus.append( + create_rad_max_hdu( + self.theta_cuts_opt["cut"][:, np.newaxis], + self.theta_bins, + self.fov_offset_bins, + ) + ) + hdus.append(fits.BinTableHDU(ang_res, name="ANGULAR_RESOLUTION")) + hdus.append(fits.BinTableHDU(bias_resolution, name="ENERGY_BIAS_RESOLUTION")) + + self.log.info("Writing outputfile") + fits.HDUList(hdus).writeto( + self.output_path / Path(self.output_file + ".fits.gz"), + overwrite=self.overwrite, + ) + + +def main(): + tool = IrfTool() + tool.run() + + +if __name__ == "main": + main() diff --git a/setup.cfg b/setup.cfg index 8085914b2d9..a00303c3d42 100644 --- a/setup.cfg +++ b/setup.cfg @@ -36,6 +36,7 @@ install_requires= numba ~=0.56.0 numpy ~=1.16 psutil + pyirf pyyaml >=5.1 requests scikit-learn @@ -93,6 +94,7 @@ console_scripts = ctapipe-display-dl1 = ctapipe.tools.display_dl1:main ctapipe-process = ctapipe.tools.process:main ctapipe-merge = ctapipe.tools.merge:main + ctapipe-make-irfs = ctapipe.tools.make-irf:main ctapipe-fileinfo = ctapipe.tools.fileinfo:main ctapipe-quickstart = ctapipe.tools.quickstart:main ctapipe-train-energy-regressor = ctapipe.tools.train_energy_regressor:main From 7970c748e3bade4f05008432f65b86539c5be407 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Wed, 19 Apr 2023 15:27:04 +0200 Subject: [PATCH 02/37] Fixed names so the tool can install properly --- ctapipe/tools/{make-irf.py => make_irf.py} | 0 setup.cfg | 2 +- 2 files changed, 1 insertion(+), 1 deletion(-) rename ctapipe/tools/{make-irf.py => make_irf.py} (100%) diff --git a/ctapipe/tools/make-irf.py b/ctapipe/tools/make_irf.py similarity index 100% rename from ctapipe/tools/make-irf.py rename to ctapipe/tools/make_irf.py diff --git a/setup.cfg b/setup.cfg index a00303c3d42..54d5b4f703a 100644 --- a/setup.cfg +++ b/setup.cfg @@ -94,7 +94,7 @@ console_scripts = ctapipe-display-dl1 = ctapipe.tools.display_dl1:main ctapipe-process = ctapipe.tools.process:main ctapipe-merge = ctapipe.tools.merge:main - ctapipe-make-irfs = ctapipe.tools.make-irf:main + ctapipe-make-irfs = ctapipe.tools.make_irf:main ctapipe-fileinfo = ctapipe.tools.fileinfo:main ctapipe-quickstart = ctapipe.tools.quickstart:main ctapipe-train-energy-regressor = ctapipe.tools.train_energy_regressor:main From dd967bdd5cceda6573206dd8fb6da37b55f7f516 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Tue, 25 Apr 2023 12:41:07 +0200 Subject: [PATCH 03/37] Fixing two small comments from Karl --- ctapipe/irf/irf_classes.py | 2 +- ctapipe/tools/make_irf.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index ce572f7c6f1..3dcb527d646 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -76,7 +76,7 @@ class IrfToolBase(Tool): 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 purity before optimisatoin" + default_value=0.4, help="Start value of gamma purity before optimisation" ).tag(config=True) energy_reconstructor = Unicode( diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index f5608f09c24..46f3a202e8f 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -224,7 +224,7 @@ def start(self): s["reco_energy_center"] ) - def finalise(self): + def finish(self): masks = { "": self.signal["selected"], From 28c66f54b902d1dba647285bc8d6d0913b1c934c Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Wed, 3 May 2023 11:03:06 +0200 Subject: [PATCH 04/37] Various small fixes --- ctapipe/irf/irf_classes.py | 8 +++----- ctapipe/tools/make_irf.py | 5 ++--- environment.yml | 1 + 3 files changed, 6 insertions(+), 8 deletions(-) diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index 3dcb527d646..82785a5b329 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -41,14 +41,12 @@ class IrfToolBase(Tool): ).tag(config=True) output_path = traits.Path( - default_value=None, + default_value="./IRF.fits.gz", allow_none=False, directory_ok=False, help="Output file", ).tag(config=True) - output_file = Unicode( - default_value=None, allow_none=False, help="Name for the output file" - ).tag(config=True) + overwrite = Bool( False, help="Overwrite the output file if it exists", @@ -133,7 +131,7 @@ def _preselect_events(self, events): rename_from.append(old) rename_to.append(new) - keep_columns.append(rename_from) + keep_columns.extend(rename_from) events = QTable(events[keep_columns], copy=False) events.rename_columns(rename_from, rename_to) keep = QualityQuery(quality_criteria=self.preselect_criteria).get_table_mask( diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index 46f3a202e8f..593466622cc 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -1,6 +1,5 @@ """Tool to generate IRFs""" import operator -from pathlib import Path import astropy.units as u import numpy as np @@ -326,9 +325,9 @@ def finish(self): hdus.append(fits.BinTableHDU(ang_res, name="ANGULAR_RESOLUTION")) hdus.append(fits.BinTableHDU(bias_resolution, name="ENERGY_BIAS_RESOLUTION")) - self.log.info("Writing outputfile") + self.log.info("Writing outputfile '%s'" % self.output_path) fits.HDUList(hdus).writeto( - self.output_path / Path(self.output_file + ".fits.gz"), + self.output_path, overwrite=self.overwrite, ) diff --git a/environment.yml b/environment.yml index e2d90f14896..2e851f10fa5 100644 --- a/environment.yml +++ b/environment.yml @@ -25,6 +25,7 @@ dependencies: - pandas - pre-commit - psutil + - pyirf - pytables - pytest - pytest-cov From 83351d8ec31632d3f8ab476fd6ea708191f6f2a5 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Wed, 3 May 2023 18:35:16 +0200 Subject: [PATCH 05/37] Big refactor to use components as intended --- ctapipe/irf/__init__.py | 5 +- ctapipe/irf/irf_classes.py | 19 +++++--- ctapipe/tools/make_irf.py | 99 +++++++++++++++++++++----------------- 3 files changed, 70 insertions(+), 53 deletions(-) diff --git a/ctapipe/irf/__init__.py b/ctapipe/irf/__init__.py index 362baef0e80..79bb3bb0c33 100644 --- a/ctapipe/irf/__init__.py +++ b/ctapipe/irf/__init__.py @@ -1,3 +1,4 @@ -from .irf_classes import DataBinning, IrfToolBase +from .irf_classes import DataBinning, ToolConfig,EventPreSelector -__all__ = ["IrfToolBase", "DataBinning"] + +__all__ = ["DataBinning","ToolConfig","EventPreSelector"] diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index 82785a5b329..c87c3c188d2 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -6,11 +6,10 @@ from astropy.table import QTable from pyirf.binning import create_bins_per_decade -from ..core import Component, QualityQuery, Tool, traits +from ..core import Component, QualityQuery, traits from ..core.traits import Bool, Float, Integer, List, Unicode - -class IrfToolBase(Tool): +class ToolConfig(Component): gamma_file = traits.Path( default_value=None, directory_ok=False, help="Gamma input filename and path" @@ -61,6 +60,9 @@ class IrfToolBase(Tool): alpha = Float( default_value=5.0, help="Ratio between size of on and off regions" ).tag(config=True) + ON_radius = Float( + default_value=1.0, help="Radius of ON region in degrees" + ).tag(config=True) max_bg_radius = Float( default_value=5.0, help="Radius used to calculate background rate in degrees" @@ -90,9 +92,10 @@ class IrfToolBase(Tool): help="Prefix of the classifier `_prediction` column", ).tag(config=True) +class EventPreSelector(Component): preselect_criteria = List( default_value=[ - ("multiplicity 4", "np.count_nonzero(tels,axis=1) >= 4"), +# ("multiplicity 4", "np.count_nonzero(tels,axis=1) >= 4"), ("valid classifier", "valid_classer"), ("valid geom reco", "valid_geom"), ("valid energy reco", "valid_energy"), @@ -120,10 +123,10 @@ def _preselect_events(self, events): "true_alt", ] rename_from = [ - f"{self.energy_reconstructor}_energy", - f"{self.geometry_reconstructor}_az", - f"{self.geometry_reconstructor}_alt", - f"{self.gammaness_classifier}_prediction", + f"{self.tc.energy_reconstructor}_energy", + f"{self.tc.geometry_reconstructor}_az", + f"{self.tc.geometry_reconstructor}_alt", + f"{self.tc.gammaness_classifier}_prediction", ] rename_to = ["reco_energy", "reco_az", "reco_alt", "gh_score"] diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index 593466622cc..7a27ea7fe6f 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -37,9 +37,9 @@ ) from pyirf.utils import calculate_source_fov_offset, calculate_theta -from ..core import Provenance +from ..core import Provenance, Tool from ..io import TableLoader -from ..irf import DataBinning, IrfToolBase +from ..irf import DataBinning, ToolConfig,EventPreSelector PYIRF_SPECTRA = { "CRAB_HEGRA": CRAB_HEGRA, @@ -48,11 +48,13 @@ } -class IrfTool(IrfToolBase, DataBinning): +class IrfTool(Tool): name = "ctapipe-make-irfs" description = "Tool to create IRF files in GAD format" - def make_derived_columns(self, events, spectrum, target_spectrum): + classes = [ DataBinning, ToolConfig,EventPreSelector] + + def make_derived_columns(self,kind, events, spectrum, target_spectrum): events["pointing_az"] = 0 * u.deg events["pointing_alt"] = 70 * u.deg @@ -68,6 +70,9 @@ def make_derived_columns(self, events, spectrum, target_spectrum): events["reco_source_fov_offset"] = calculate_source_fov_offset( events, prefix="reco" ) + # Gamma source is assumed to be pointlike + if kind == "gamma": + spectrum = spectrum.integrate_cone(0*u.deg,self.tc.ON_radius*u.deg) events["weights"] = calculate_event_weights( events["true_energy"], target_spectrum=target_spectrum, @@ -79,6 +84,8 @@ def make_derived_columns(self, events, spectrum, target_spectrum): def get_sim_info_and_spectrum(self, loader): sim = loader.read_simulation_configuration() + # These sims better have the same viewcone! + assert( sim["max_viewcone_radius"].std() == 0) sim_info = SimulatedEventsInfo( n_showers=sum(sim["n_showers"] * sim["shower_reuse"]), energy_min=sim["energy_range_min"].quantity[0], @@ -89,37 +96,43 @@ def get_sim_info_and_spectrum(self, loader): ) return sim_info, PowerLaw.from_simulation( - sim_info, obstime=self.obs_time * u.Unit(self.obs_time_unit) + sim_info, obstime=self.tc.obs_time * u.Unit(self.tc.obs_time_unit) ) def setup(self): + self.tc = ToolConfig(parent=self) + self.bins = DataBinning(parent=self) + self.eps = EventPreSelector(parent=self) + opts = dict(load_dl2=True, load_simulated=True, load_dl1_parameters=False) reduced_events = dict() for kind, file, target_spectrum in [ - ("gamma", self.gamma_file, PYIRF_SPECTRA[self.gamma_sim_spectrum]), - ("proton", self.proton_file, PYIRF_SPECTRA[self.proton_sim_spectrum]), - ("electron", self.electron_file, PYIRF_SPECTRA[self.electron_sim_spectrum]), + ("gamma", self.tc.gamma_file, PYIRF_SPECTRA[self.tc.gamma_sim_spectrum]), + ("proton", self.tc.proton_file, PYIRF_SPECTRA[self.tc.proton_sim_spectrum]), + ("electron", self.tc.electron_file, PYIRF_SPECTRA[self.tc.electron_sim_spectrum]), ]: with TableLoader(file, **opts) as load: Provenance().add_input_file(file) - table = self._make_empty_table() + table = self.eps._make_empty_table() sim_info, spectrum = self.get_sim_info_and_spectrum(load) if kind == "gamma": self.sim_info = sim_info self.spectrum = spectrum for start, stop, events in load.read_subarray_events_chunked( - self.chunk_size + self.tc.chunk_size ): - selected = self._preselect_events(events) - selected = self.make_derived_columns( + selected = self.eps._preselect_events(events) + selected = self.eps.make_derived_columns( + kind, selected, spectrum, target_spectrum ) - table = vstack(table, selected) + table = vstack([table, selected]) reduced_events[kind] = table - self.signal = reduced_events["gamma"] - self.background = vstack(reduced_events["proton"], reduced_events["electron"]) + select_ON = reduced_events["gamma"]["theta"] <= self.ON_radius*u.deg + self.signal = reduced_events["gamma"][select_ON] + self.background = vstack([reduced_events["proton"], reduced_events["electron"]]) self.theta_bins = add_overflow_bins( create_bins_per_decade( @@ -127,16 +140,16 @@ def setup(self): ) ) - self.energy_reco_bins = self.reco_energy_bins() - self.energy_true_bins = self.true_energy_bins() - self.source_offset_bins = self.source_offset_bins() - self.fov_offset_bins = self.fov_offset_bins() - self.energy_migration_bins = self.energy_migration_bins() + self.reco_energy_bins = self.bins.reco_energy_bins() + self.true_energy_bins = self.bins.true_energy_bins() + self.source_offset_bins = self.bins.source_offset_bins() + self.fov_offset_bins = self.bins.fov_offset_bins() + self.energy_migration_bins = self.bins.energy_migration_bins() def start(self): - + breakpoint() INITIAL_GH_CUT = np.quantile( - self.signal["gh_score"], (1 - self.initial_gh_cut_efficency) + self.signal["gh_score"], (1 - self.tc.initial_gh_cut_efficency) ) self.log.info( f"Using fixed G/H cut of {INITIAL_GH_CUT} to calculate theta cuts" @@ -148,29 +161,29 @@ def start(self): self.signal["theta"][mask_theta_cuts], self.signal["reco_energy"][mask_theta_cuts], bins=self.theta_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, + 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, percentile=68, ) self.log.info("Optimizing G/H separation cut for best sensitivity") 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, + self.tc.gh_cut_efficiency_step, + self.tc.max_gh_cut_efficiency + self.tc.gh_cut_efficiency_step / 2, + self.tc.gh_cut_efficiency_step, ) sens2, self.gh_cuts = optimize_gh_cut( self.signal, self.background, - reco_energy_bins=self.energy_reco_bins, + reco_energy_bins=self.bins.reco_energy_bins, gh_cut_efficiencies=gh_cut_efficiencies, op=operator.ge, theta_cuts=theta_cuts, - alpha=self.alpha, - background_radius=self.max_bg_radius * u.deg, + alpha=self.tc.alpha, + background_radius=self.tc.max_bg_radius * u.deg, ) # now that we have the optimized gh cuts, we recalculate the theta @@ -186,10 +199,10 @@ def start(self): self.signal[self.signal["selected_gh"]]["reco_energy"], self.theta_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, + 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"], @@ -201,20 +214,21 @@ def start(self): self.signal["selected_theta"] & self.signal["selected_gh"] ) + breakpoint() # calculate sensitivity signal_hist = create_histogram_table( - self.signal[self.signal["selected"]], bins=self.energy_reco_bins + self.signal[self.signal["selected"]], bins=self.reco_energy_bins ) background_hist = estimate_background( self.background[self.background["selected_gh"]], - reco_energy_bins=self.energy_reco_bins, + reco_energy_bins=self.reco_energy_bins, theta_cuts=self.theta_cuts_opt, - alpha=self.alpha, - fov_offset_min=self.fov_offset_min, - fov_offset_max=self.fov_offset_max, + alpha=self.tc.alpha, + fov_offset_min=self.bins.fov_offset_min, + fov_offset_max=self.bins.fov_offset_max, ) self.sensitivity = calculate_sensitivity( - signal_hist, background_hist, alpha=self.alpha + signal_hist, background_hist, alpha=self.tc.alpha ) # scale relative sensitivity by Crab flux to get the flux sensitivity @@ -224,7 +238,6 @@ def start(self): ) def finish(self): - masks = { "": self.signal["selected"], "_NO_CUTS": slice(None), From b17948cdf8efe6b9f30602d4f32ec7bdc6bb53f7 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Thu, 4 May 2023 14:54:58 +0200 Subject: [PATCH 06/37] Got it all to run --- ctapipe/irf/__init__.py | 4 ++-- ctapipe/irf/irf_classes.py | 18 ++++++++++-------- ctapipe/tools/make_irf.py | 34 +++++++++++++++++----------------- 3 files changed, 29 insertions(+), 27 deletions(-) diff --git a/ctapipe/irf/__init__.py b/ctapipe/irf/__init__.py index 79bb3bb0c33..08839dfd523 100644 --- a/ctapipe/irf/__init__.py +++ b/ctapipe/irf/__init__.py @@ -1,4 +1,4 @@ -from .irf_classes import DataBinning, ToolConfig,EventPreSelector +from .irf_classes import DataBinning, ToolConfig, EventPreProcessor -__all__ = ["DataBinning","ToolConfig","EventPreSelector"] +__all__ = ["DataBinning","ToolConfig","EventPreProcessor"] diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index c87c3c188d2..346fd0ae5f2 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -72,13 +72,15 @@ class ToolConfig(Component): default_value=0.8, help="Maximum gamma purity requested" ).tag(config=True) gh_cut_efficiency_step = Float( - default_value=0.01, + 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 purity before optimisation" ).tag(config=True) + +class EventPreProcessor(Component): energy_reconstructor = Unicode( default_value="RandomForestRegressor", help="Prefix of the reco `_energy` column", @@ -92,7 +94,6 @@ class ToolConfig(Component): help="Prefix of the classifier `_prediction` column", ).tag(config=True) -class EventPreSelector(Component): preselect_criteria = List( default_value=[ # ("multiplicity 4", "np.count_nonzero(tels,axis=1) >= 4"), @@ -115,6 +116,7 @@ class EventPreSelector(Component): ) def _preselect_events(self, events): + tc = self.parent.tc keep_columns = [ "obs_id", "event_id", @@ -123,10 +125,10 @@ def _preselect_events(self, events): "true_alt", ] rename_from = [ - f"{self.tc.energy_reconstructor}_energy", - f"{self.tc.geometry_reconstructor}_az", - f"{self.tc.geometry_reconstructor}_alt", - f"{self.tc.gammaness_classifier}_prediction", + f"{self.energy_reconstructor}_energy", + f"{self.geometry_reconstructor}_az", + f"{self.geometry_reconstructor}_alt", + f"{self.gammaness_classifier}_prediction", ] rename_to = ["reco_energy", "reco_az", "reco_alt", "gh_score"] @@ -268,12 +270,12 @@ class DataBinning(Component): ).tag(config=True) fov_offset_min = Float( - help="Minimum value for FoV Offset bins", + help="Minimum value for FoV Offset bins in degrees", default_value=0.1, ).tag(config=True) fov_offset_max = Float( - help="Maximum value for FoV offset bins", + help="Maximum value for FoV offset bins in degrees", default_value=1.1, ).tag(config=True) diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index 7a27ea7fe6f..b820317b4fe 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -39,7 +39,7 @@ from ..core import Provenance, Tool from ..io import TableLoader -from ..irf import DataBinning, ToolConfig,EventPreSelector +from ..irf import DataBinning, ToolConfig, EventPreProcessor PYIRF_SPECTRA = { "CRAB_HEGRA": CRAB_HEGRA, @@ -52,7 +52,7 @@ class IrfTool(Tool): name = "ctapipe-make-irfs" description = "Tool to create IRF files in GAD format" - classes = [ DataBinning, ToolConfig,EventPreSelector] + classes = [ DataBinning, ToolConfig,EventPreProcessor] def make_derived_columns(self,kind, events, spectrum, target_spectrum): events["pointing_az"] = 0 * u.deg @@ -73,7 +73,7 @@ def make_derived_columns(self,kind, events, spectrum, target_spectrum): # Gamma source is assumed to be pointlike if kind == "gamma": spectrum = spectrum.integrate_cone(0*u.deg,self.tc.ON_radius*u.deg) - events["weights"] = calculate_event_weights( + events["weight"] = calculate_event_weights( events["true_energy"], target_spectrum=target_spectrum, simulated_spectrum=spectrum, @@ -102,7 +102,7 @@ def get_sim_info_and_spectrum(self, loader): def setup(self): self.tc = ToolConfig(parent=self) self.bins = DataBinning(parent=self) - self.eps = EventPreSelector(parent=self) + self.eps = EventPreProcessor(parent=self) opts = dict(load_dl2=True, load_simulated=True, load_dl1_parameters=False) reduced_events = dict() @@ -122,7 +122,7 @@ def setup(self): self.tc.chunk_size ): selected = self.eps._preselect_events(events) - selected = self.eps.make_derived_columns( + selected = self.make_derived_columns( kind, selected, spectrum, target_spectrum ) @@ -130,7 +130,7 @@ def setup(self): reduced_events[kind] = table - select_ON = reduced_events["gamma"]["theta"] <= self.ON_radius*u.deg + select_ON = reduced_events["gamma"]["theta"] <= self.tc.ON_radius*u.deg self.signal = reduced_events["gamma"][select_ON] self.background = vstack([reduced_events["proton"], reduced_events["electron"]]) @@ -144,10 +144,10 @@ def setup(self): self.true_energy_bins = self.bins.true_energy_bins() self.source_offset_bins = self.bins.source_offset_bins() self.fov_offset_bins = self.bins.fov_offset_bins() + self.bkg_fov_offset_bins = self.bins.bkg_fov_offset_bins() self.energy_migration_bins = self.bins.energy_migration_bins() def start(self): - breakpoint() INITIAL_GH_CUT = np.quantile( self.signal["gh_score"], (1 - self.tc.initial_gh_cut_efficency) ) @@ -178,12 +178,12 @@ def start(self): sens2, self.gh_cuts = optimize_gh_cut( self.signal, self.background, - reco_energy_bins=self.bins.reco_energy_bins, + reco_energy_bins=self.reco_energy_bins, gh_cut_efficiencies=gh_cut_efficiencies, op=operator.ge, theta_cuts=theta_cuts, alpha=self.tc.alpha, - background_radius=self.tc.max_bg_radius * u.deg, + fov_offset_max=self.tc.max_bg_radius * u.deg, ) # now that we have the optimized gh cuts, we recalculate the theta @@ -214,18 +214,18 @@ def start(self): self.signal["selected_theta"] & self.signal["selected_gh"] ) - breakpoint() # calculate sensitivity signal_hist = create_histogram_table( self.signal[self.signal["selected"]], bins=self.reco_energy_bins ) + background_hist = estimate_background( self.background[self.background["selected_gh"]], reco_energy_bins=self.reco_energy_bins, theta_cuts=self.theta_cuts_opt, alpha=self.tc.alpha, - fov_offset_min=self.bins.fov_offset_min, - fov_offset_max=self.bins.fov_offset_max, + fov_offset_min=self.bins.fov_offset_min*u.deg, + fov_offset_max=self.bins.fov_offset_max*u.deg, ) self.sensitivity = calculate_sensitivity( signal_hist, background_hist, alpha=self.tc.alpha @@ -308,15 +308,15 @@ def finish(self): background_rate = background_2d( self.background[self.background["selected_gh"]], self.reco_energy_bins, - fov_offset_bins=np.arange(0, 11) * u.deg, - t_obs=self.obs_time * u.Unit(self.obs_time_unit), + fov_offset_bins=self.bkg_fov_offset_bins, + t_obs=self.tc.obs_time * u.Unit(self.tc.obs_time_unit), ) hdus.append( create_background_2d_hdu( background_rate, self.reco_energy_bins, - fov_offset_bins=np.arange(0, 11) * u.deg, + fov_offset_bins=self.bkg_fov_offset_bins, ) ) @@ -338,9 +338,9 @@ def finish(self): hdus.append(fits.BinTableHDU(ang_res, name="ANGULAR_RESOLUTION")) hdus.append(fits.BinTableHDU(bias_resolution, name="ENERGY_BIAS_RESOLUTION")) - self.log.info("Writing outputfile '%s'" % self.output_path) + self.log.info("Writing outputfile '%s'" % self.tc.output_path) fits.HDUList(hdus).writeto( - self.output_path, + self.tc.output_path, overwrite=self.overwrite, ) From 8b4607bb43de8f1aab951d68fe3714ca01bfa63a Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Thu, 4 May 2023 14:57:07 +0200 Subject: [PATCH 07/37] Formatting fixes --- ctapipe/irf/__init__.py | 5 ++--- ctapipe/irf/irf_classes.py | 10 +++++----- ctapipe/tools/make_irf.py | 27 +++++++++++++++------------ 3 files changed, 22 insertions(+), 20 deletions(-) diff --git a/ctapipe/irf/__init__.py b/ctapipe/irf/__init__.py index 08839dfd523..19b3ac99230 100644 --- a/ctapipe/irf/__init__.py +++ b/ctapipe/irf/__init__.py @@ -1,4 +1,3 @@ -from .irf_classes import DataBinning, ToolConfig, EventPreProcessor +from .irf_classes import DataBinning, EventPreProcessor, ToolConfig - -__all__ = ["DataBinning","ToolConfig","EventPreProcessor"] +__all__ = ["DataBinning", "ToolConfig", "EventPreProcessor"] diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index 346fd0ae5f2..7f95c4b9a59 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -9,6 +9,7 @@ from ..core import Component, QualityQuery, traits from ..core.traits import Bool, Float, Integer, List, Unicode + class ToolConfig(Component): gamma_file = traits.Path( @@ -60,9 +61,9 @@ class ToolConfig(Component): alpha = Float( default_value=5.0, help="Ratio between size of on and off regions" ).tag(config=True) - ON_radius = Float( - default_value=1.0, help="Radius of ON region in degrees" - ).tag(config=True) + ON_radius = Float(default_value=1.0, help="Radius of ON region in degrees").tag( + config=True + ) max_bg_radius = Float( default_value=5.0, help="Radius used to calculate background rate in degrees" @@ -96,7 +97,7 @@ class EventPreProcessor(Component): preselect_criteria = List( default_value=[ -# ("multiplicity 4", "np.count_nonzero(tels,axis=1) >= 4"), + # ("multiplicity 4", "np.count_nonzero(tels,axis=1) >= 4"), ("valid classifier", "valid_classer"), ("valid geom reco", "valid_geom"), ("valid energy reco", "valid_energy"), @@ -116,7 +117,6 @@ class EventPreProcessor(Component): ) def _preselect_events(self, events): - tc = self.parent.tc keep_columns = [ "obs_id", "event_id", diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index b820317b4fe..010d2e748e6 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -39,7 +39,7 @@ from ..core import Provenance, Tool from ..io import TableLoader -from ..irf import DataBinning, ToolConfig, EventPreProcessor +from ..irf import DataBinning, EventPreProcessor, ToolConfig PYIRF_SPECTRA = { "CRAB_HEGRA": CRAB_HEGRA, @@ -52,9 +52,9 @@ class IrfTool(Tool): name = "ctapipe-make-irfs" description = "Tool to create IRF files in GAD format" - classes = [ DataBinning, ToolConfig,EventPreProcessor] + classes = [DataBinning, ToolConfig, EventPreProcessor] - def make_derived_columns(self,kind, events, spectrum, target_spectrum): + def make_derived_columns(self, kind, events, spectrum, target_spectrum): events["pointing_az"] = 0 * u.deg events["pointing_alt"] = 70 * u.deg @@ -70,9 +70,9 @@ def make_derived_columns(self,kind, events, spectrum, target_spectrum): events["reco_source_fov_offset"] = calculate_source_fov_offset( events, prefix="reco" ) - # Gamma source is assumed to be pointlike + # Gamma source is assumed to be pointlike if kind == "gamma": - spectrum = spectrum.integrate_cone(0*u.deg,self.tc.ON_radius*u.deg) + spectrum = spectrum.integrate_cone(0 * u.deg, self.tc.ON_radius * u.deg) events["weight"] = calculate_event_weights( events["true_energy"], target_spectrum=target_spectrum, @@ -85,7 +85,7 @@ def get_sim_info_and_spectrum(self, loader): sim = loader.read_simulation_configuration() # These sims better have the same viewcone! - assert( sim["max_viewcone_radius"].std() == 0) + assert sim["max_viewcone_radius"].std() == 0 sim_info = SimulatedEventsInfo( n_showers=sum(sim["n_showers"] * sim["shower_reuse"]), energy_min=sim["energy_range_min"].quantity[0], @@ -109,7 +109,11 @@ def setup(self): for kind, file, target_spectrum in [ ("gamma", self.tc.gamma_file, PYIRF_SPECTRA[self.tc.gamma_sim_spectrum]), ("proton", self.tc.proton_file, PYIRF_SPECTRA[self.tc.proton_sim_spectrum]), - ("electron", self.tc.electron_file, PYIRF_SPECTRA[self.tc.electron_sim_spectrum]), + ( + "electron", + self.tc.electron_file, + PYIRF_SPECTRA[self.tc.electron_sim_spectrum], + ), ]: with TableLoader(file, **opts) as load: Provenance().add_input_file(file) @@ -123,14 +127,13 @@ def setup(self): ): selected = self.eps._preselect_events(events) selected = self.make_derived_columns( - kind, - selected, spectrum, target_spectrum + kind, selected, spectrum, target_spectrum ) table = vstack([table, selected]) reduced_events[kind] = table - select_ON = reduced_events["gamma"]["theta"] <= self.tc.ON_radius*u.deg + select_ON = reduced_events["gamma"]["theta"] <= self.tc.ON_radius * u.deg self.signal = reduced_events["gamma"][select_ON] self.background = vstack([reduced_events["proton"], reduced_events["electron"]]) @@ -224,8 +227,8 @@ def start(self): reco_energy_bins=self.reco_energy_bins, theta_cuts=self.theta_cuts_opt, alpha=self.tc.alpha, - fov_offset_min=self.bins.fov_offset_min*u.deg, - fov_offset_max=self.bins.fov_offset_max*u.deg, + fov_offset_min=self.bins.fov_offset_min * u.deg, + fov_offset_max=self.bins.fov_offset_max * u.deg, ) self.sensitivity = calculate_sensitivity( signal_hist, background_hist, alpha=self.tc.alpha From 616d66cc777705610022256ba719de78b77b1aac Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Thu, 4 May 2023 16:39:15 +0200 Subject: [PATCH 08/37] Adjusting defaults to make only one bin in offset? --- ctapipe/irf/irf_classes.py | 2 +- ctapipe/tools/make_irf.py | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index 7f95c4b9a59..84d778bb66c 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -281,7 +281,7 @@ class DataBinning(Component): fov_offset_n_edges = Integer( help="Number of edges for FoV offset bins", - default_value=9, + default_value=2, ).tag(config=True) bkg_fov_offset_min = Float( diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index 010d2e748e6..023e2d3d34c 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -262,6 +262,8 @@ def finish(self): self.sim_info, true_energy_bins=self.true_energy_bins, ) + self.log.debug(self.true_energy_bins) + self.log.debug(self.fov_offset_bins) hdus.append( create_aeff2d_hdu( effective_area[..., np.newaxis], # +1 dimension for FOV offset From a2a65a7cccac57ca49189c541cd5e861c6dc87cf Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Wed, 10 May 2023 14:03:58 +0200 Subject: [PATCH 09/37] Moved preselection to start() --- ctapipe/tools/make_irf.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index 023e2d3d34c..b600984450a 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -99,11 +99,7 @@ def get_sim_info_and_spectrum(self, loader): sim_info, obstime=self.tc.obs_time * u.Unit(self.tc.obs_time_unit) ) - def setup(self): - self.tc = ToolConfig(parent=self) - self.bins = DataBinning(parent=self) - self.eps = EventPreProcessor(parent=self) - + def load_preselected_events(self): opts = dict(load_dl2=True, load_simulated=True, load_dl1_parameters=False) reduced_events = dict() for kind, file, target_spectrum in [ @@ -137,6 +133,11 @@ def setup(self): self.signal = reduced_events["gamma"][select_ON] self.background = vstack([reduced_events["proton"], reduced_events["electron"]]) + def setup(self): + self.tc = ToolConfig(parent=self) + self.bins = DataBinning(parent=self) + self.eps = EventPreProcessor(parent=self) + self.theta_bins = add_overflow_bins( create_bins_per_decade( self.sim_info.energy_min, self.sim_info.energy_max, 50 @@ -151,6 +152,8 @@ def setup(self): self.energy_migration_bins = self.bins.energy_migration_bins() def start(self): + self.load_preselected_events() + INITIAL_GH_CUT = np.quantile( self.signal["gh_score"], (1 - self.tc.initial_gh_cut_efficency) ) From d850a76710906a71503e19a00efecf04c62cf562 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Wed, 10 May 2023 14:57:35 +0200 Subject: [PATCH 10/37] Use only one true energy axis --- ctapipe/tools/make_irf.py | 14 ++------------ 1 file changed, 2 insertions(+), 12 deletions(-) diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index b600984450a..5611a23d8e6 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -6,11 +6,7 @@ from astropy.io import fits from astropy.table import vstack from pyirf.benchmarks import angular_resolution, energy_bias_resolution -from pyirf.binning import ( - add_overflow_bins, - create_bins_per_decade, - create_histogram_table, -) +from pyirf.binning import create_histogram_table from pyirf.cut_optimization import optimize_gh_cut from pyirf.cuts import calculate_percentile_cut, evaluate_binned_cut from pyirf.io import ( @@ -138,12 +134,6 @@ def setup(self): self.bins = DataBinning(parent=self) self.eps = EventPreProcessor(parent=self) - self.theta_bins = add_overflow_bins( - create_bins_per_decade( - self.sim_info.energy_min, self.sim_info.energy_max, 50 - ) - ) - self.reco_energy_bins = self.bins.reco_energy_bins() self.true_energy_bins = self.bins.true_energy_bins() self.source_offset_bins = self.bins.source_offset_bins() @@ -166,7 +156,7 @@ def start(self): theta_cuts = calculate_percentile_cut( self.signal["theta"][mask_theta_cuts], self.signal["reco_energy"][mask_theta_cuts], - bins=self.theta_bins, + bins=self.true_energy_bins, 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, From 1781cf78acb544a506a6f1b18cdcd6b102aa5d8d Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Wed, 10 May 2023 15:09:18 +0200 Subject: [PATCH 11/37] Use only one true energy axis --- ctapipe/tools/make_irf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index 5611a23d8e6..ad38dfec940 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -193,7 +193,7 @@ def start(self): self.theta_cuts_opt = calculate_percentile_cut( self.signal[self.signal["selected_gh"]]["theta"], self.signal[self.signal["selected_gh"]]["reco_energy"], - self.theta_bins, + self.true_energy_bins, percentile=68, min_value=self.bins.theta_min_angle * u.deg, max_value=self.bins.theta_max_angle * u.deg, From e06d375b4e682b2b476e8890aac3fc50fbccdb9d Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Wed, 10 May 2023 15:38:02 +0200 Subject: [PATCH 12/37] Use only one true energy axis --- ctapipe/tools/make_irf.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index ad38dfec940..cb30a54bba6 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -326,10 +326,11 @@ def finish(self): self.fov_offset_bins, ) ) + hdus.append( create_rad_max_hdu( self.theta_cuts_opt["cut"][:, np.newaxis], - self.theta_bins, + self.true_energy_bins, self.fov_offset_bins, ) ) From 77c4056378138512ca35201c88e8da9144fdca0e Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Wed, 10 May 2023 15:38:45 +0200 Subject: [PATCH 13/37] Refactoring to group things in the same manner --- ctapipe/tools/make_irf.py | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index cb30a54bba6..c9c59c78dd8 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -287,6 +287,7 @@ def finish(self): self.reco_energy_bins, energy_type="reco", ) + hdus.append(fits.BinTableHDU(bias_resolution, name="ENERGY_BIAS_RESOLUTION")) # Here we use reconstructed energy instead of true energy for the sake of # current pipelines comparisons @@ -295,13 +296,7 @@ def finish(self): self.reco_energy_bins, energy_type="reco", ) - - psf = psf_table( - self.signal[self.signal["selected_gh"]], - self.true_energy_bins, - fov_offset_bins=self.fov_offset_bins, - source_offset_bins=self.source_offset_bins, - ) + hdus.append(fits.BinTableHDU(ang_res, name="ANGULAR_RESOLUTION")) background_rate = background_2d( self.background[self.background["selected_gh"]], @@ -309,7 +304,6 @@ def finish(self): fov_offset_bins=self.bkg_fov_offset_bins, t_obs=self.tc.obs_time * u.Unit(self.tc.obs_time_unit), ) - hdus.append( create_background_2d_hdu( background_rate, @@ -318,6 +312,12 @@ def finish(self): ) ) + psf = psf_table( + self.signal[self.signal["selected_gh"]], + self.true_energy_bins, + fov_offset_bins=self.fov_offset_bins, + source_offset_bins=self.source_offset_bins, + ) hdus.append( create_psf_table_hdu( psf, @@ -334,8 +334,6 @@ def finish(self): self.fov_offset_bins, ) ) - hdus.append(fits.BinTableHDU(ang_res, name="ANGULAR_RESOLUTION")) - hdus.append(fits.BinTableHDU(bias_resolution, name="ENERGY_BIAS_RESOLUTION")) self.log.info("Writing outputfile '%s'" % self.tc.output_path) fits.HDUList(hdus).writeto( From 8efa3b9448667693677217e597ec05c0775929a3 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Fri, 12 May 2023 16:51:13 +0200 Subject: [PATCH 14/37] Split preselct into remaning and quality selecting, made EventPreProcessor a quality query --- ctapipe/irf/irf_classes.py | 32 +++++++++++++------------------- ctapipe/tools/make_irf.py | 5 +++-- 2 files changed, 16 insertions(+), 21 deletions(-) diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index 84d778bb66c..d6dc6fac134 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -81,7 +81,7 @@ class ToolConfig(Component): ).tag(config=True) -class EventPreProcessor(Component): +class EventPreProcessor(QualityQuery): energy_reconstructor = Unicode( default_value="RandomForestRegressor", help="Prefix of the reco `_energy` column", @@ -97,26 +97,22 @@ class EventPreProcessor(Component): preselect_criteria = List( default_value=[ - # ("multiplicity 4", "np.count_nonzero(tels,axis=1) >= 4"), - ("valid classifier", "valid_classer"), - ("valid geom reco", "valid_geom"), - ("valid energy reco", "valid_energy"), + ("multiplicity 4", "subarray.multiplicity(tels_with_trigger) >= 4"), + ("valid classifier", "RandomForestClassifier_is_valid"), + ("valid geom reco", "HillasReconstructor_is_valid"), + ("valid energy reco", "RandomForestRegressor_is_valid"), ], help=QualityQuery.quality_criteria.help, ).tag(config=True) rename_columns = List( - help="List containing translation pairs of quality columns" - "used for quality filters and their names as given in the input file used." + help="List containing translation pairs new and old column names" + "used when processing input with names differing from the CTA prod5b format" "Ex: [('valid_geom','HillasReconstructor_is_valid')]", - default_value=[ - ("valid_geom", "HillasReconstructor_is_valid"), - ("valid_energy", "RandomForestRegressor_is_valid"), - ("valid_classer", "RandomForestClassifier_is_valid"), - ], + default_value=[], ) - def _preselect_events(self, events): + def normalise_column_names(self, events): keep_columns = [ "obs_id", "event_id", @@ -132,6 +128,7 @@ def _preselect_events(self, events): ] rename_to = ["reco_energy", "reco_az", "reco_alt", "gh_score"] + # We never enter the loop if rename_columns is empty for new, old in self.rename_columns: rename_from.append(old) rename_to.append(new) @@ -139,13 +136,10 @@ def _preselect_events(self, events): keep_columns.extend(rename_from) events = QTable(events[keep_columns], copy=False) events.rename_columns(rename_from, rename_to) - keep = QualityQuery(quality_criteria=self.preselect_criteria).get_table_mask( - events - ) - - return events[keep] + return events - def _make_empty_table(self): + def make_empty_table(self): + """This function defines the columns later functions expect to be present in the event table""" columns = [ "obs_id", "event_id", diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index c9c59c78dd8..1a57f323524 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -109,7 +109,7 @@ def load_preselected_events(self): ]: with TableLoader(file, **opts) as load: Provenance().add_input_file(file) - table = self.eps._make_empty_table() + table = self.eps.make_empty_table() sim_info, spectrum = self.get_sim_info_and_spectrum(load) if kind == "gamma": self.sim_info = sim_info @@ -117,7 +117,8 @@ def load_preselected_events(self): for start, stop, events in load.read_subarray_events_chunked( self.tc.chunk_size ): - selected = self.eps._preselect_events(events) + selected = self.eps.normalise_column_names(events) + selected = selected[self.eps.get_table_mask(selected)] selected = self.make_derived_columns( kind, selected, spectrum, target_spectrum ) From b7eca8d20ae1bfd53b61dbe21bdc110d65e7be02 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Wed, 23 Aug 2023 18:39:37 +0200 Subject: [PATCH 15/37] Fixed serveral problems --- ctapipe/irf/irf_classes.py | 6 ++++-- ctapipe/tools/make_irf.py | 28 ++++++++++++++++++---------- 2 files changed, 22 insertions(+), 12 deletions(-) diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index d6dc6fac134..19d5ab100a4 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -152,9 +152,10 @@ def make_empty_table(self): "gh_score", "pointing_az", "pointing_alt", + "theta", "true_source_fov_offset", "reco_source_fov_offset", - "weights", + "weight", ] units = [ None, @@ -170,6 +171,7 @@ def make_empty_table(self): u.deg, u.deg, u.deg, + u.deg, None, ] @@ -265,7 +267,7 @@ class DataBinning(Component): fov_offset_min = Float( help="Minimum value for FoV Offset bins in degrees", - default_value=0.1, + default_value=0.0, ).tag(config=True) fov_offset_max = Float( diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index 1a57f323524..ee7449dc97c 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -50,16 +50,21 @@ class IrfTool(Tool): classes = [DataBinning, ToolConfig, EventPreProcessor] - def make_derived_columns(self, kind, events, spectrum, target_spectrum): - events["pointing_az"] = 0 * u.deg - events["pointing_alt"] = 70 * u.deg + def make_derived_columns(self, kind, events, spectrum, target_spectrum, obs_conf): + + if obs_conf["subarray_pointing_lat"].std() < 1e-3: + assert obs_conf["subarray_pointing_frame"] == 0 + # Lets suppose 0 means ALTAZ + events["pointing_alt"] = obs["subarray_pointing_lat"][0] + events["pointing_az"] = obs["subarray_pointing_lon"][0] + else: + raise NotImplemented("No support for making irfs from varying pointings yet") events["theta"] = calculate_theta( events, assumed_source_az=events["true_az"], assumed_source_alt=events["true_alt"], ) - events["true_source_fov_offset"] = calculate_source_fov_offset( events, prefix="true" ) @@ -77,7 +82,8 @@ def make_derived_columns(self, kind, events, spectrum, target_spectrum): return events - def get_sim_info_and_spectrum(self, loader): + def get_metadata(self, loader): + obs = loader.read_observation_information() sim = loader.read_simulation_configuration() # These sims better have the same viewcone! @@ -93,7 +99,7 @@ def get_sim_info_and_spectrum(self, loader): return sim_info, PowerLaw.from_simulation( sim_info, obstime=self.tc.obs_time * u.Unit(self.tc.obs_time_unit) - ) + ), obs def load_preselected_events(self): opts = dict(load_dl2=True, load_simulated=True, load_dl1_parameters=False) @@ -109,21 +115,23 @@ def load_preselected_events(self): ]: with TableLoader(file, **opts) as load: Provenance().add_input_file(file) - table = self.eps.make_empty_table() - sim_info, spectrum = self.get_sim_info_and_spectrum(load) + header = self.eps.make_empty_table() + sim_info, spectrum, obs_conf = self.get_metadata(load) if kind == "gamma": self.sim_info = sim_info self.spectrum = spectrum + bits = [header] for start, stop, events in load.read_subarray_events_chunked( self.tc.chunk_size ): selected = self.eps.normalise_column_names(events) selected = selected[self.eps.get_table_mask(selected)] selected = self.make_derived_columns( - kind, selected, spectrum, target_spectrum + kind, selected, spectrum, target_spectrum, obs_conf ) - table = vstack([table, selected]) + bits.append(selected) + table = vstack(bits,join_type="exact") reduced_events[kind] = table select_ON = reduced_events["gamma"]["theta"] <= self.tc.ON_radius * u.deg From b2fd6167dcf9e53650dd6337b4c9f32932e5f955 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Wed, 23 Aug 2023 22:07:58 +0200 Subject: [PATCH 16/37] Fixed unit issue --- ctapipe/tools/make_irf.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index ee7449dc97c..25905904193 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -53,10 +53,10 @@ class IrfTool(Tool): def make_derived_columns(self, kind, events, spectrum, target_spectrum, obs_conf): if obs_conf["subarray_pointing_lat"].std() < 1e-3: - assert obs_conf["subarray_pointing_frame"] == 0 + assert all( obs_conf["subarray_pointing_frame"] == 0) # Lets suppose 0 means ALTAZ - events["pointing_alt"] = obs["subarray_pointing_lat"][0] - events["pointing_az"] = obs["subarray_pointing_lon"][0] + events["pointing_alt"] = obs_conf["subarray_pointing_lat"][0]*u.deg + events["pointing_az"] = obs_conf["subarray_pointing_lon"][0]*u.deg else: raise NotImplemented("No support for making irfs from varying pointings yet") From 0bb8a80f33485174b728b45c931d0fa62d05bc3f Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Wed, 23 Aug 2023 22:35:54 +0200 Subject: [PATCH 17/37] Updated some defaults --- ctapipe/irf/irf_classes.py | 8 ++++---- ctapipe/tools/make_irf.py | 25 ++++++++++++++++--------- 2 files changed, 20 insertions(+), 13 deletions(-) diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index 19d5ab100a4..ba4ac9f5a03 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -20,14 +20,14 @@ class ToolConfig(Component): help="Name of the pyrif spectra used for the simulated gamma spectrum", ).tag(config=True) proton_file = traits.Path( - default_value=None, directory_ok=False, help="Gamma input filename and path" + default_value=None, directory_ok=False, help="Proton input filename and path" ).tag(config=True) proton_sim_spectrum = traits.Unicode( default_value="IRFDOC_PROTON_SPECTRUM", help="Name of the pyrif spectra used for the simulated proton spectrum", ).tag(config=True) electron_file = traits.Path( - default_value=None, directory_ok=False, help="Gamma input filename and path" + default_value=None, directory_ok=False, help="Electron input filename and path" ).tag(config=True) electron_sim_spectrum = traits.Unicode( default_value="IRFDOC_ELECTRON_SPECTRUM", @@ -59,7 +59,7 @@ class ToolConfig(Component): ).tag(config=True) alpha = Float( - default_value=5.0, help="Ratio between size of on and off regions" + default_value=0.2, help="Ratio between size of on and off regions" ).tag(config=True) ON_radius = Float(default_value=1.0, help="Radius of ON region in degrees").tag( config=True @@ -215,7 +215,7 @@ class DataBinning(Component): true_energy_n_bins_per_decade = Float( help="Number of edges per decade for True Energy bins", - default_value=5, + default_value=10, ).tag(config=True) reco_energy_min = Float( diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index 25905904193..067d61c16b7 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -53,12 +53,14 @@ class IrfTool(Tool): def make_derived_columns(self, kind, events, spectrum, target_spectrum, obs_conf): if obs_conf["subarray_pointing_lat"].std() < 1e-3: - assert all( obs_conf["subarray_pointing_frame"] == 0) + assert all(obs_conf["subarray_pointing_frame"] == 0) # Lets suppose 0 means ALTAZ - events["pointing_alt"] = obs_conf["subarray_pointing_lat"][0]*u.deg - events["pointing_az"] = obs_conf["subarray_pointing_lon"][0]*u.deg + events["pointing_alt"] = obs_conf["subarray_pointing_lat"][0] * u.deg + events["pointing_az"] = obs_conf["subarray_pointing_lon"][0] * u.deg else: - raise NotImplemented("No support for making irfs from varying pointings yet") + raise NotImplementedError( + "No support for making irfs from varying pointings yet" + ) events["theta"] = calculate_theta( events, @@ -97,9 +99,13 @@ def get_metadata(self, loader): viewcone=sim["max_viewcone_radius"].quantity[0], ) - return sim_info, PowerLaw.from_simulation( - sim_info, obstime=self.tc.obs_time * u.Unit(self.tc.obs_time_unit) - ), obs + return ( + sim_info, + PowerLaw.from_simulation( + sim_info, obstime=self.tc.obs_time * u.Unit(self.tc.obs_time_unit) + ), + obs, + ) def load_preselected_events(self): opts = dict(load_dl2=True, load_simulated=True, load_dl1_parameters=False) @@ -131,7 +137,7 @@ def load_preselected_events(self): ) bits.append(selected) - table = vstack(bits,join_type="exact") + table = vstack(bits, join_type="exact") reduced_events[kind] = table select_ON = reduced_events["gamma"]["theta"] <= self.tc.ON_radius * u.deg @@ -280,6 +286,7 @@ def finish(self): fov_offset_bins=self.fov_offset_bins, migration_bins=self.energy_migration_bins, ) + breakpoint() hdus.append( create_energy_dispersion_hdu( edisp, @@ -338,7 +345,7 @@ def finish(self): hdus.append( create_rad_max_hdu( - self.theta_cuts_opt["cut"][:, np.newaxis], + self.theta_cuts_opt["cut"].reshape(-1, 1), self.true_energy_bins, self.fov_offset_bins, ) From 791e6fbfe40aee5fa5bcbdaf7624de3165bd3250 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Mon, 28 Aug 2023 15:20:03 +0200 Subject: [PATCH 18/37] Updated some comments to avoid problems with the doctest --- ctapipe/irf/irf_classes.py | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index ba4ac9f5a03..2122ad79d12 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -313,13 +313,6 @@ class DataBinning(Component): def true_energy_bins(self): """ Creates bins per decade for true MC energy using pyirf function. - The overflow binning added is not needed at the current stage. - - Examples - -------- - It can be used as: - - >>> add_overflow_bins(***)[1:-1] """ true_energy = create_bins_per_decade( self.true_energy_min * u.TeV, @@ -331,13 +324,6 @@ def true_energy_bins(self): def reco_energy_bins(self): """ Creates bins per decade for reconstructed MC energy using pyirf function. - The overflow binning added is not needed at the current stage. - - Examples - -------- - It can be used as: - - >>> add_overflow_bins(***)[1:-1] """ reco_energy = create_bins_per_decade( self.reco_energy_min * u.TeV, From ac410a9c43c1c031cc97b4441b24ee0a2fac4c23 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Mon, 28 Aug 2023 15:30:02 +0200 Subject: [PATCH 19/37] Add changelog --- docs/changes/2315.irf-maker.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 docs/changes/2315.irf-maker.rst diff --git a/docs/changes/2315.irf-maker.rst b/docs/changes/2315.irf-maker.rst new file mode 100644 index 00000000000..37a898ba437 --- /dev/null +++ b/docs/changes/2315.irf-maker.rst @@ -0,0 +1 @@ +Add a `make-irf tool` able to produce irfs given a gamma, proton and electron DL2 input files. From 15588cdf3f5f8669fedc9fcc8643bcf49f615d55 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Mon, 28 Aug 2023 15:41:39 +0200 Subject: [PATCH 20/37] Fixed small bug, trying to pass doctests --- ctapipe/irf/irf_classes.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index 2122ad79d12..432d83be716 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -110,7 +110,7 @@ class EventPreProcessor(QualityQuery): "used when processing input with names differing from the CTA prod5b format" "Ex: [('valid_geom','HillasReconstructor_is_valid')]", default_value=[], - ) + ).tag(config=True) def normalise_column_names(self, events): keep_columns = [ From e2b1e02c39df6f90d48bc5db0c18db2cad39a9e5 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Tue, 5 Sep 2023 12:01:32 +0200 Subject: [PATCH 21/37] Refactored after feedback --- ctapipe/irf/__init__.py | 4 +- ctapipe/irf/irf_classes.py | 184 +++++++++++++------------------------ ctapipe/tools/make_irf.py | 119 ++++++++++++++++++------ 3 files changed, 156 insertions(+), 151 deletions(-) diff --git a/ctapipe/irf/__init__.py b/ctapipe/irf/__init__.py index 19b3ac99230..b1056dbb807 100644 --- a/ctapipe/irf/__init__.py +++ b/ctapipe/irf/__init__.py @@ -1,3 +1,3 @@ -from .irf_classes import DataBinning, EventPreProcessor, ToolConfig +from .irf_classes import CutOptimising, DataBinning, EnergyBinning, EventPreProcessor -__all__ = ["DataBinning", "ToolConfig", "EventPreProcessor"] +__all__ = ["CutOptimising", "DataBinning", "EnergyBinning", "EventPreProcessor"] diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index 432d83be716..a524093c75b 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -6,82 +6,28 @@ from astropy.table import QTable from pyirf.binning import create_bins_per_decade -from ..core import Component, QualityQuery, traits -from ..core.traits import Bool, Float, Integer, List, Unicode +from ..core import Component, QualityQuery +from ..core.traits import Float, Integer, List, Unicode -class ToolConfig(Component): - - gamma_file = traits.Path( - default_value=None, directory_ok=False, help="Gamma input filename and path" - ).tag(config=True) - gamma_sim_spectrum = traits.Unicode( - default_value="CRAB_HEGRA", - help="Name of the pyrif spectra used for the simulated gamma spectrum", - ).tag(config=True) - proton_file = traits.Path( - default_value=None, directory_ok=False, help="Proton input filename and path" - ).tag(config=True) - proton_sim_spectrum = traits.Unicode( - default_value="IRFDOC_PROTON_SPECTRUM", - help="Name of the pyrif spectra used for the simulated proton spectrum", - ).tag(config=True) - electron_file = traits.Path( - default_value=None, directory_ok=False, help="Electron input filename and path" - ).tag(config=True) - electron_sim_spectrum = traits.Unicode( - default_value="IRFDOC_ELECTRON_SPECTRUM", - help="Name of the pyrif spectra used for the simulated electron spectrum", - ).tag(config=True) - - chunk_size = Integer( - default_value=100000, - allow_none=True, - help="How many subarray events to load at once for making predictions.", - ).tag(config=True) - - output_path = traits.Path( - default_value="./IRF.fits.gz", - allow_none=False, - directory_ok=False, - help="Output file", - ).tag(config=True) - - overwrite = Bool( - False, - help="Overwrite the output file if it exists", - ).tag(config=True) - - obs_time = Float(default_value=50.0, help="Observation time").tag(config=True) - obs_time_unit = Unicode( - default_value="hour", - help="Unit used to specify observation time as an astropy unit string.", - ).tag(config=True) - - alpha = Float( - default_value=0.2, help="Ratio between size of on and off regions" - ).tag(config=True) - ON_radius = Float(default_value=1.0, help="Radius of ON region in degrees").tag( - config=True - ) - - max_bg_radius = Float( - default_value=5.0, help="Radius used to calculate background rate in degrees" - ).tag(config=True) +class CutOptimising(Component): + """Collects settings related to the cut configuration""" max_gh_cut_efficiency = Float( - default_value=0.8, help="Maximum gamma purity requested" + 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 purity before optimisation" + default_value=0.4, help="Start value of gamma efficiency before optimisation" ).tag(config=True) class EventPreProcessor(QualityQuery): + """Defines preselection cuts and the necessary renaming of columns""" + energy_reconstructor = Unicode( default_value="RandomForestRegressor", help="Prefix of the reco `_energy` column", @@ -157,23 +103,19 @@ def make_empty_table(self): "reco_source_fov_offset", "weight", ] - units = [ - None, - None, - u.TeV, - u.deg, - u.deg, - u.TeV, - u.deg, - u.deg, - None, - u.deg, - u.deg, - u.deg, - u.deg, - u.deg, - None, - ] + units = { + "true_energy": u.TeV, + "true_az": u.deg, + "true_alt": u.deg, + "reco_energy": u.TeV, + "reco_az": u.deg, + "reco_alt": u.deg, + "pointing_az": u.deg, + "pointing_alt": u.deg, + "theta": u.deg, + "true_source_fov_offset": u.deg, + "reco_source_fov_offset": u.deg, + } return QTable(names=columns, units=units) @@ -195,13 +137,8 @@ class ThetaSettings(Component): ).tag(config=True) -class DataBinning(Component): - """ - Collects information on generating energy and angular bins for - generating IRFs as per pyIRF requirements. - - Stolen from LSTChain - """ +class EnergyBinning(Component): + """Collects energy binning settings""" true_energy_min = Float( help="Minimum value for True Energy bins in TeV units", @@ -248,6 +185,48 @@ class DataBinning(Component): default_value=31, ).tag(config=True) + def true_energy_bins(self): + """ + Creates bins per decade for true MC energy using pyirf function. + """ + true_energy = create_bins_per_decade( + self.true_energy_min * u.TeV, + self.true_energy_max * u.TeV, + self.true_energy_n_bins_per_decade, + ) + return true_energy + + 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 energy_migration_bins(self): + """ + Creates bins for energy migration. + """ + energy_migration = np.geomspace( + self.energy_migration_min, + self.energy_migration_max, + self.energy_migration_n_bins, + ) + return energy_migration + + +class DataBinning(Component): + """ + Collects information on generating energy and angular bins for + generating IRFs as per pyIRF requirements. + + Stolen from LSTChain + """ + theta_min_angle = Float( default_value=0.05, help="Smallest angular cut value allowed" ).tag(config=True) @@ -310,39 +289,6 @@ class DataBinning(Component): default_value=101, ).tag(config=True) - def true_energy_bins(self): - """ - Creates bins per decade for true MC energy using pyirf function. - """ - true_energy = create_bins_per_decade( - self.true_energy_min * u.TeV, - self.true_energy_max * u.TeV, - self.true_energy_n_bins_per_decade, - ) - return true_energy - - 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 energy_migration_bins(self): - """ - Creates bins for energy migration. - """ - energy_migration = np.geomspace( - self.energy_migration_min, - self.energy_migration_max, - self.energy_migration_n_bins, - ) - return energy_migration - def fov_offset_bins(self): """ Creates bins for single/multiple FoV offset. diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index 067d61c16b7..2fbc50969d8 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -33,9 +33,10 @@ ) from pyirf.utils import calculate_source_fov_offset, calculate_theta -from ..core import Provenance, Tool +from ..core import Provenance, Tool, traits +from ..core.traits import Bool, Float, Integer, Unicode from ..io import TableLoader -from ..irf import DataBinning, EventPreProcessor, ToolConfig +from ..irf import CutOptimising, DataBinning, EnergyBinning, EventPreProcessor PYIRF_SPECTRA = { "CRAB_HEGRA": CRAB_HEGRA, @@ -48,7 +49,63 @@ class IrfTool(Tool): name = "ctapipe-make-irfs" description = "Tool to create IRF files in GAD format" - classes = [DataBinning, ToolConfig, EventPreProcessor] + gamma_file = traits.Path( + default_value=None, directory_ok=False, help="Gamma input filename and path" + ).tag(config=True) + gamma_sim_spectrum = traits.Unicode( + default_value="CRAB_HEGRA", + help="Name of the pyrif spectra used for the simulated gamma spectrum", + ).tag(config=True) + proton_file = traits.Path( + default_value=None, directory_ok=False, help="Proton input filename and path" + ).tag(config=True) + proton_sim_spectrum = traits.Unicode( + default_value="IRFDOC_PROTON_SPECTRUM", + help="Name of the pyrif spectra used for the simulated proton spectrum", + ).tag(config=True) + electron_file = traits.Path( + default_value=None, directory_ok=False, help="Electron input filename and path" + ).tag(config=True) + electron_sim_spectrum = traits.Unicode( + default_value="IRFDOC_ELECTRON_SPECTRUM", + help="Name of the pyrif spectra used for the simulated electron spectrum", + ).tag(config=True) + + chunk_size = Integer( + default_value=100000, + allow_none=True, + help="How many subarray events to load at once for making predictions.", + ).tag(config=True) + + output_path = traits.Path( + default_value="./IRF.fits.gz", + allow_none=False, + directory_ok=False, + help="Output file", + ).tag(config=True) + + overwrite = Bool( + False, + help="Overwrite the output file if it exists", + ).tag(config=True) + + obs_time = Float(default_value=50.0, help="Observation time").tag(config=True) + obs_time_unit = Unicode( + default_value="hour", + help="Unit used to specify observation time as an astropy unit string.", + ).tag(config=True) + + alpha = Float( + default_value=0.2, help="Ratio between size of on and off regions" + ).tag(config=True) + ON_radius = Float(default_value=1.0, help="Radius of ON region in degrees").tag( + config=True + ) + max_bg_radius = Float( + default_value=3.0, help="Radius used to calculate background rate in degrees" + ).tag(config=True) + + classes = [CutOptimising, DataBinning, EnergyBinning, EventPreProcessor] def make_derived_columns(self, kind, events, spectrum, target_spectrum, obs_conf): @@ -75,7 +132,7 @@ def make_derived_columns(self, kind, events, spectrum, target_spectrum, obs_conf ) # Gamma source is assumed to be pointlike if kind == "gamma": - spectrum = spectrum.integrate_cone(0 * u.deg, self.tc.ON_radius * u.deg) + spectrum = spectrum.integrate_cone(0 * u.deg, self.ON_radius * u.deg) events["weight"] = calculate_event_weights( events["true_energy"], target_spectrum=target_spectrum, @@ -102,7 +159,7 @@ def get_metadata(self, loader): return ( sim_info, PowerLaw.from_simulation( - sim_info, obstime=self.tc.obs_time * u.Unit(self.tc.obs_time_unit) + sim_info, obstime=self.obs_time * u.Unit(self.obs_time_unit) ), obs, ) @@ -111,27 +168,27 @@ def load_preselected_events(self): opts = dict(load_dl2=True, load_simulated=True, load_dl1_parameters=False) reduced_events = dict() for kind, file, target_spectrum in [ - ("gamma", self.tc.gamma_file, PYIRF_SPECTRA[self.tc.gamma_sim_spectrum]), - ("proton", self.tc.proton_file, PYIRF_SPECTRA[self.tc.proton_sim_spectrum]), + ("gamma", self.gamma_file, PYIRF_SPECTRA[self.gamma_sim_spectrum]), + ("proton", self.proton_file, PYIRF_SPECTRA[self.proton_sim_spectrum]), ( "electron", - self.tc.electron_file, - PYIRF_SPECTRA[self.tc.electron_sim_spectrum], + self.electron_file, + PYIRF_SPECTRA[self.electron_sim_spectrum], ), ]: with TableLoader(file, **opts) as load: Provenance().add_input_file(file) - header = self.eps.make_empty_table() + header = self.epp.make_empty_table() sim_info, spectrum, obs_conf = self.get_metadata(load) if kind == "gamma": self.sim_info = sim_info self.spectrum = spectrum bits = [header] for start, stop, events in load.read_subarray_events_chunked( - self.tc.chunk_size + self.chunk_size ): - selected = self.eps.normalise_column_names(events) - selected = selected[self.eps.get_table_mask(selected)] + selected = self.epp.normalise_column_names(events) + selected = selected[self.epp.get_table_mask(selected)] selected = self.make_derived_columns( kind, selected, spectrum, target_spectrum, obs_conf ) @@ -140,27 +197,29 @@ def load_preselected_events(self): table = vstack(bits, join_type="exact") reduced_events[kind] = table - select_ON = reduced_events["gamma"]["theta"] <= self.tc.ON_radius * u.deg + select_ON = reduced_events["gamma"]["theta"] <= self.ON_radius * u.deg self.signal = reduced_events["gamma"][select_ON] self.background = vstack([reduced_events["proton"], reduced_events["electron"]]) def setup(self): - self.tc = ToolConfig(parent=self) + self.co = CutOptimising(parent=self) + self.e_bins = EnergyBinning(parent=self) self.bins = DataBinning(parent=self) - self.eps = EventPreProcessor(parent=self) + self.epp = EventPreProcessor(parent=self) + + self.reco_energy_bins = self.e_bins.reco_energy_bins() + self.true_energy_bins = self.e_bins.true_energy_bins() + self.energy_migration_bins = self.e_bins.energy_migration_bins() - self.reco_energy_bins = self.bins.reco_energy_bins() - self.true_energy_bins = self.bins.true_energy_bins() self.source_offset_bins = self.bins.source_offset_bins() self.fov_offset_bins = self.bins.fov_offset_bins() self.bkg_fov_offset_bins = self.bins.bkg_fov_offset_bins() - self.energy_migration_bins = self.bins.energy_migration_bins() def start(self): self.load_preselected_events() INITIAL_GH_CUT = np.quantile( - self.signal["gh_score"], (1 - self.tc.initial_gh_cut_efficency) + self.signal["gh_score"], (1 - self.co.initial_gh_cut_efficency) ) self.log.info( f"Using fixed G/H cut of {INITIAL_GH_CUT} to calculate theta cuts" @@ -181,9 +240,9 @@ def start(self): self.log.info("Optimizing G/H separation cut for best sensitivity") gh_cut_efficiencies = np.arange( - self.tc.gh_cut_efficiency_step, - self.tc.max_gh_cut_efficiency + self.tc.gh_cut_efficiency_step / 2, - self.tc.gh_cut_efficiency_step, + self.co.gh_cut_efficiency_step, + self.co.max_gh_cut_efficiency + self.co.gh_cut_efficiency_step / 2, + self.co.gh_cut_efficiency_step, ) sens2, self.gh_cuts = optimize_gh_cut( @@ -193,8 +252,8 @@ def start(self): gh_cut_efficiencies=gh_cut_efficiencies, op=operator.ge, theta_cuts=theta_cuts, - alpha=self.tc.alpha, - fov_offset_max=self.tc.max_bg_radius * u.deg, + alpha=self.alpha, + fov_offset_max=self.max_bg_radius * u.deg, ) # now that we have the optimized gh cuts, we recalculate the theta @@ -234,12 +293,12 @@ def start(self): self.background[self.background["selected_gh"]], reco_energy_bins=self.reco_energy_bins, theta_cuts=self.theta_cuts_opt, - alpha=self.tc.alpha, + alpha=self.alpha, fov_offset_min=self.bins.fov_offset_min * u.deg, fov_offset_max=self.bins.fov_offset_max * u.deg, ) self.sensitivity = calculate_sensitivity( - signal_hist, background_hist, alpha=self.tc.alpha + signal_hist, background_hist, alpha=self.alpha ) # scale relative sensitivity by Crab flux to get the flux sensitivity @@ -318,7 +377,7 @@ def finish(self): self.background[self.background["selected_gh"]], self.reco_energy_bins, fov_offset_bins=self.bkg_fov_offset_bins, - t_obs=self.tc.obs_time * u.Unit(self.tc.obs_time_unit), + t_obs=self.obs_time * u.Unit(self.obs_time_unit), ) hdus.append( create_background_2d_hdu( @@ -351,9 +410,9 @@ def finish(self): ) ) - self.log.info("Writing outputfile '%s'" % self.tc.output_path) + self.log.info("Writing outputfile '%s'" % self.output_path) fits.HDUList(hdus).writeto( - self.tc.output_path, + self.output_path, overwrite=self.overwrite, ) From 972e2cf3afa7012395c715df9ec06f18afea1189 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Tue, 5 Sep 2023 13:42:42 +0200 Subject: [PATCH 22/37] Refactored after feedback --- ctapipe/irf/irf_classes.py | 46 ++++++++++++++++++++++++++++++++++++++ ctapipe/tools/make_irf.py | 42 +++++----------------------------- 2 files changed, 52 insertions(+), 36 deletions(-) diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index a524093c75b..54afe51e2cd 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -1,10 +1,14 @@ """ Define a parent IrfTool class to hold all the options """ +import operator + import astropy.units as u import numpy as np 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 ..core import Component, QualityQuery from ..core.traits import Float, Integer, List, Unicode @@ -24,6 +28,48 @@ class CutOptimising(Component): 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 + ): + INITIAL_GH_CUT = np.quantile( + signal["gh_score"], (1 - self.initial_gh_cut_efficency) + ) + self.log.info( + f"Using fixed G/H cut of {INITIAL_GH_CUT} to calculate theta cuts" + ) + + mask_theta_cuts = signal["gh_score"] >= INITIAL_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, + percentile=68, + ) + + self.log.info("Optimizing G/H separation cut for best sensitivity") + 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, + ) + + sens2, gh_cuts = optimize_gh_cut( + signal, + background, + reco_energy_bins=Er_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 + class EventPreProcessor(QualityQuery): """Defines preselection cuts and the necessary renaming of columns""" diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index 2fbc50969d8..f055041c8d4 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -7,7 +7,6 @@ from astropy.table import vstack from pyirf.benchmarks import angular_resolution, energy_bias_resolution from pyirf.binning import create_histogram_table -from pyirf.cut_optimization import optimize_gh_cut from pyirf.cuts import calculate_percentile_cut, evaluate_binned_cut from pyirf.io import ( create_aeff2d_hdu, @@ -217,43 +216,14 @@ def setup(self): def start(self): self.load_preselected_events() - - INITIAL_GH_CUT = np.quantile( - self.signal["gh_score"], (1 - self.co.initial_gh_cut_efficency) - ) - self.log.info( - f"Using fixed G/H cut of {INITIAL_GH_CUT} to calculate theta cuts" - ) - - mask_theta_cuts = self.signal["gh_score"] >= INITIAL_GH_CUT - - theta_cuts = calculate_percentile_cut( - self.signal["theta"][mask_theta_cuts], - self.signal["reco_energy"][mask_theta_cuts], - bins=self.true_energy_bins, - 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, - percentile=68, - ) - - self.log.info("Optimizing G/H separation cut for best sensitivity") - gh_cut_efficiencies = np.arange( - self.co.gh_cut_efficiency_step, - self.co.max_gh_cut_efficiency + self.co.gh_cut_efficiency_step / 2, - self.co.gh_cut_efficiency_step, - ) - - sens2, self.gh_cuts = optimize_gh_cut( + self.gh_cuts, sens2 = self.co.optimise_gh_cut( self.signal, self.background, - reco_energy_bins=self.reco_energy_bins, - gh_cut_efficiencies=gh_cut_efficiencies, - op=operator.ge, - theta_cuts=theta_cuts, - alpha=self.alpha, - fov_offset_max=self.max_bg_radius * u.deg, + 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 From 8ba9c68cf7b44f208f0c39c95df0dc9ec88c4209 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Tue, 5 Sep 2023 13:45:01 +0200 Subject: [PATCH 23/37] Refactored after feedback --- ctapipe/irf/irf_classes.py | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index 54afe51e2cd..6c4658d0658 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -166,23 +166,6 @@ def make_empty_table(self): return QTable(names=columns, units=units) -class ThetaSettings(Component): - - min_angle = Float( - default_value=0.05, help="Smallest angular cut value allowed" - ).tag(config=True) - max_angle = Float(default_value=0.32, help="Largest angular cut value allowed").tag( - config=True - ) - min_counts = Integer( - default_value=10, - help="Minimum number of events in a bin to attempt to find a cut value", - ).tag(config=True) - fill_value = Float( - default_value=0.32, help="Angular cut value used for bins with too few events" - ).tag(config=True) - - class EnergyBinning(Component): """Collects energy binning settings""" From bfdad602e57fd4ad0cc8f65d66c9d1ff855a9bce Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Tue, 5 Sep 2023 14:07:08 +0200 Subject: [PATCH 24/37] Made spectra into enums --- ctapipe/tools/make_irf.py | 30 ++++++++++++++++++++---------- 1 file changed, 20 insertions(+), 10 deletions(-) diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index f055041c8d4..cbf4ae84ce4 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -1,5 +1,6 @@ """Tool to generate IRFs""" import operator +from enum import Enum import astropy.units as u import numpy as np @@ -37,10 +38,17 @@ from ..io import TableLoader from ..irf import CutOptimising, DataBinning, EnergyBinning, EventPreProcessor + +class Spectra(Enum): + CRAB_HEGRA = 1 + IRFDOC_ELECTRON_SPECTRUM = 2 + IRFDOC_PROTON_SPECTRUM = 3 + + PYIRF_SPECTRA = { - "CRAB_HEGRA": CRAB_HEGRA, - "IRFDOC_ELECTRON_SPECTRUM": IRFDOC_ELECTRON_SPECTRUM, - "IRFDOC_PROTON_SPECTRUM": IRFDOC_PROTON_SPECTRUM, + Spectra.CRAB_HEGRA: CRAB_HEGRA, + Spectra.IRFDOC_ELECTRON_SPECTRUM: IRFDOC_ELECTRON_SPECTRUM, + Spectra.IRFDOC_PROTON_SPECTRUM: IRFDOC_PROTON_SPECTRUM, } @@ -51,22 +59,25 @@ class IrfTool(Tool): gamma_file = traits.Path( default_value=None, directory_ok=False, help="Gamma input filename and path" ).tag(config=True) - gamma_sim_spectrum = traits.Unicode( - default_value="CRAB_HEGRA", + gamma_sim_spectrum = traits.UseEnum( + Spectra, + default_value=Spectra.CRAB_HEGRA, help="Name of the pyrif spectra used for the simulated gamma spectrum", ).tag(config=True) proton_file = traits.Path( default_value=None, directory_ok=False, help="Proton input filename and path" ).tag(config=True) - proton_sim_spectrum = traits.Unicode( - default_value="IRFDOC_PROTON_SPECTRUM", + proton_sim_spectrum = traits.UseEnum( + Spectra, + default_value=Spectra.IRFDOC_PROTON_SPECTRUM, help="Name of the pyrif spectra used for the simulated proton spectrum", ).tag(config=True) electron_file = traits.Path( default_value=None, directory_ok=False, help="Electron input filename and path" ).tag(config=True) - electron_sim_spectrum = traits.Unicode( - default_value="IRFDOC_ELECTRON_SPECTRUM", + electron_sim_spectrum = traits.UseEnum( + Spectra, + default_value=Spectra.IRFDOC_ELECTRON_SPECTRUM, help="Name of the pyrif spectra used for the simulated electron spectrum", ).tag(config=True) @@ -315,7 +326,6 @@ def finish(self): fov_offset_bins=self.fov_offset_bins, migration_bins=self.energy_migration_bins, ) - breakpoint() hdus.append( create_energy_dispersion_hdu( edisp, From 29b3f196d9a4d1a3069322963657fe4411c1b4df Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Thu, 7 Sep 2023 14:54:48 +0200 Subject: [PATCH 25/37] 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"], From 6972b080cad1655676eaec014456c496a2ecbd1b Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Fri, 8 Sep 2023 15:39:49 +0200 Subject: [PATCH 26/37] Refactored based on feedback: new tests that simulations stay consistent and some renaming --- ctapipe/irf/__init__.py | 4 +-- ctapipe/irf/irf_classes.py | 2 +- ctapipe/tools/make_irf.py | 74 +++++++++++++++++++++++--------------- 3 files changed, 49 insertions(+), 31 deletions(-) diff --git a/ctapipe/irf/__init__.py b/ctapipe/irf/__init__.py index 6cd70eb6a8c..825fcf698cd 100644 --- a/ctapipe/irf/__init__.py +++ b/ctapipe/irf/__init__.py @@ -1,8 +1,8 @@ from .irf_classes import ( - CutOptimising, + CutOptimizer, DataBinning, EventPreProcessor, OutputEnergyBinning, ) -__all__ = ["CutOptimising", "DataBinning", "OutputEnergyBinning", "EventPreProcessor"] +__all__ = ["CutOptimizer", "DataBinning", "OutputEnergyBinning", "EventPreProcessor"] diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index 88a1e470e61..f4a8b90d5af 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -14,7 +14,7 @@ from ..core.traits import Float, Integer, List, Unicode -class CutOptimising(Component): +class CutOptimizer(Component): """Performs cut optimisation""" max_gh_cut_efficiency = Float( diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index f9ad3871ca6..d20f8ccafae 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -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, EventPreProcessor, OutputEnergyBinning +from ..irf import CutOptimizer, 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, OutputEnergyBinning, EventPreProcessor] + classes = [CutOptimizer, DataBinning, OutputEnergyBinning, EventPreProcessor] def make_derived_columns(self, kind, events, spectrum, target_spectrum, obs_conf): @@ -154,11 +154,25 @@ def make_derived_columns(self, kind, events, spectrum, target_spectrum, obs_conf def get_metadata(self, loader): obs = loader.read_observation_information() sim = loader.read_simulation_configuration() + show = loader.read_shower_distribution() # These sims better have the same viewcone! + if not np.diff(sim["energy_range_max"]).sum() == 0: + raise NotImplementedError( + "Unsupported: 'energy_range_max' differs across simulation runs" + ) + if not np.diff(sim["energy_range_min"]).sum() == 0: + raise NotImplementedError( + "Unsupported: 'energy_range_min' differs across simulation runs" + ) + if not np.diff(sim["spectral_index"]).sum() == 0: + raise NotImplementedError( + "Unsupported: 'spectral_index' differs across simulation runs" + ) + assert sim["max_viewcone_radius"].std() == 0 sim_info = SimulatedEventsInfo( - n_showers=sum(sim["n_showers"] * sim["shower_reuse"]), + n_showers=show["n_entries"].sum(), energy_min=sim["energy_range_min"].quantity[0], energy_max=sim["energy_range_max"].quantity[0], max_impact=sim["max_scatter_range"].quantity[0], @@ -208,11 +222,13 @@ def load_preselected_events(self): reduced_events[kind] = table select_ON = reduced_events["gamma"]["theta"] <= self.ON_radius * u.deg - self.signal = reduced_events["gamma"][select_ON] - self.background = vstack([reduced_events["proton"], reduced_events["electron"]]) + self.signal_events = reduced_events["gamma"][select_ON] + self.background_events = vstack( + [reduced_events["proton"], reduced_events["electron"]] + ) def setup(self): - self.co = CutOptimising(parent=self) + self.co = CutOptimizer(parent=self) self.e_bins = OutputEnergyBinning(parent=self) self.bins = DataBinning(parent=self) self.epp = EventPreProcessor(parent=self) @@ -227,30 +243,31 @@ def setup(self): def start(self): self.load_preselected_events() - self.gh_cuts, self.theta_cuts_opt, sens2 = self.co.optimise_gh_cut( - self.signal, - self.background, + self.gh_cuts, self.theta_cuts_opt, self.sens2 = self.co.optimise_gh_cut( + self.signal_events, + self.background_events, self.alpha, self.max_bg_radius, ) - self.signal["selected_theta"] = evaluate_binned_cut( - self.signal["theta"], - self.signal["reco_energy"], + self.signal_events["selected_theta"] = evaluate_binned_cut( + self.signal_events["theta"], + self.signal_events["reco_energy"], self.theta_cuts_opt, operator.le, ) - self.signal["selected"] = ( - self.signal["selected_theta"] & self.signal["selected_gh"] + self.signal_events["selected"] = ( + self.signal_events["selected_theta"] & self.signal_events["selected_gh"] ) # calculate sensitivity signal_hist = create_histogram_table( - self.signal[self.signal["selected"]], bins=self.reco_energy_bins + self.signal_events[self.signal_events["selected"]], + bins=self.reco_energy_bins, ) background_hist = estimate_background( - self.background[self.background["selected_gh"]], + self.background_events[self.background_events["selected_gh"]], reco_energy_bins=self.reco_energy_bins, theta_cuts=self.theta_cuts_opt, alpha=self.alpha, @@ -262,30 +279,30 @@ def start(self): ) # scale relative sensitivity by Crab flux to get the flux sensitivity - for s in (sens2, self.sensitivity): + for s in (self.sens2, self.sensitivity): s["flux_sensitivity"] = s["relative_sensitivity"] * self.spectrum( s["reco_energy_center"] ) def finish(self): masks = { - "": self.signal["selected"], + "": self.signal_events["selected"], "_NO_CUTS": slice(None), - "_ONLY_GH": self.signal["selected_gh"], - "_ONLY_THETA": self.signal["selected_theta"], + "_ONLY_GH": self.signal_events["selected_gh"], + "_ONLY_THETA": self.signal_events["selected_theta"], } hdus = [ fits.PrimaryHDU(), fits.BinTableHDU(self.sensitivity, name="SENSITIVITY"), - # fits.BinTableHDU(sensitivity_step_2, name="SENSITIVITY_STEP_2"), - # fits.BinTableHDU(self.theta_cuts, name="THETA_CUTS"), + fits.BinTableHDU(self.sens2, name="SENSITIVITY_STEP_2"), + fits.BinTableHDU(self.theta_cuts, name="THETA_CUTS"), fits.BinTableHDU(self.theta_cuts_opt, name="THETA_CUTS_OPT"), fits.BinTableHDU(self.gh_cuts, name="GH_CUTS"), ] for label, mask in masks.items(): effective_area = effective_area_per_energy( - self.signal[mask], + self.signal_events[mask], self.sim_info, true_energy_bins=self.true_energy_bins, ) @@ -300,7 +317,7 @@ def finish(self): ) ) edisp = energy_dispersion( - self.signal[mask], + self.signal_events[mask], true_energy_bins=self.true_energy_bins, fov_offset_bins=self.fov_offset_bins, migration_bins=self.energy_migration_bins, @@ -317,7 +334,7 @@ def finish(self): # Here we use reconstructed energy instead of true energy for the sake of # current pipelines comparisons bias_resolution = energy_bias_resolution( - self.signal[self.signal["selected"]], + self.signal_events[self.signal_events["selected"]], self.reco_energy_bins, energy_type="reco", ) @@ -326,14 +343,14 @@ def finish(self): # Here we use reconstructed energy instead of true energy for the sake of # current pipelines comparisons ang_res = angular_resolution( - self.signal[self.signal["selected_gh"]], + self.signal_events[self.signal_events["selected_gh"]], self.reco_energy_bins, energy_type="reco", ) hdus.append(fits.BinTableHDU(ang_res, name="ANGULAR_RESOLUTION")) background_rate = background_2d( - self.background[self.background["selected_gh"]], + self.background_events[self.background_events["selected_gh"]], self.reco_energy_bins, fov_offset_bins=self.bkg_fov_offset_bins, t_obs=self.obs_time * u.Unit(self.obs_time_unit), @@ -347,7 +364,7 @@ def finish(self): ) psf = psf_table( - self.signal[self.signal["selected_gh"]], + self.signal_events[self.signal_events["selected_gh"]], self.true_energy_bins, fov_offset_bins=self.fov_offset_bins, source_offset_bins=self.source_offset_bins, @@ -374,6 +391,7 @@ def finish(self): self.output_path, overwrite=self.overwrite, ) + Provenance().add_output_file(self.output_path, role="IRF") def main(): From 4c9f681af3c5fa610851bf5f1b23ec7e23c1f69b Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Fri, 8 Sep 2023 15:47:48 +0200 Subject: [PATCH 27/37] Refactored based on feedback: cleaner test --- ctapipe/tools/make_irf.py | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index d20f8ccafae..ac1e9519199 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -118,7 +118,6 @@ class IrfTool(Tool): classes = [CutOptimizer, DataBinning, OutputEnergyBinning, EventPreProcessor] def make_derived_columns(self, kind, events, spectrum, target_spectrum, obs_conf): - if obs_conf["subarray_pointing_lat"].std() < 1e-3: assert all(obs_conf["subarray_pointing_frame"] == 0) # Lets suppose 0 means ALTAZ @@ -156,19 +155,11 @@ def get_metadata(self, loader): sim = loader.read_simulation_configuration() show = loader.read_shower_distribution() - # These sims better have the same viewcone! - if not np.diff(sim["energy_range_max"]).sum() == 0: - raise NotImplementedError( - "Unsupported: 'energy_range_max' differs across simulation runs" - ) - if not np.diff(sim["energy_range_min"]).sum() == 0: - raise NotImplementedError( - "Unsupported: 'energy_range_min' differs across simulation runs" - ) - if not np.diff(sim["spectral_index"]).sum() == 0: - raise NotImplementedError( - "Unsupported: 'spectral_index' differs across simulation runs" - ) + for itm in ["spectral_index", "energy_range_min", "energy_range_max"]: + if len(np.unique(sim[itm])) > 1: + raise NotImplementedError( + f"Unsupported: '{itm}' differs across simulation runs" + ) assert sim["max_viewcone_radius"].std() == 0 sim_info = SimulatedEventsInfo( From b6eb6ff488c93ff7038c05613421c7fc55853599 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Mon, 11 Sep 2023 18:42:06 +0200 Subject: [PATCH 28/37] Various changes to match reference script. --- ctapipe/irf/irf_classes.py | 6 ++-- ctapipe/tools/make_irf.py | 57 ++++++++++++++++++++++++++++---------- 2 files changed, 46 insertions(+), 17 deletions(-) diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index f4a8b90d5af..bf2409135ce 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -149,9 +149,9 @@ class EventPreProcessor(QualityQuery): help="Prefix of the classifier `_prediction` column", ).tag(config=True) - preselect_criteria = List( + quality_criteria = List( default_value=[ - ("multiplicity 4", "subarray.multiplicity(tels_with_trigger) >= 4"), + ("multiplicity 4", "np.count_nonzero(tels_with_trigger,axis=1) >= 4"), ("valid classifier", "RandomForestClassifier_is_valid"), ("valid geom reco", "HillasReconstructor_is_valid"), ("valid energy reco", "RandomForestRegressor_is_valid"), @@ -325,7 +325,7 @@ class DataBinning(Component): fov_offset_max = Float( help="Maximum value for FoV offset bins in degrees", - default_value=1.1, + default_value=2.0, ).tag(config=True) fov_offset_n_edges = Integer( diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index ac1e9519199..b02d0601349 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -18,7 +18,7 @@ ) from pyirf.irf import ( background_2d, - effective_area_per_energy, + effective_area_per_energy_and_fov, energy_dispersion, psf_table, ) @@ -108,9 +108,9 @@ class IrfTool(Tool): alpha = Float( default_value=0.2, help="Ratio between size of on and off regions" ).tag(config=True) - ON_radius = Float(default_value=1.0, help="Radius of ON region in degrees").tag( - config=True - ) + # ON_radius = Float(default_value=1.0, help="Radius of ON region in degrees").tag( + # config=True + # ) max_bg_radius = Float( default_value=3.0, help="Radius used to calculate background rate in degrees" ).tag(config=True) @@ -139,9 +139,12 @@ def make_derived_columns(self, kind, events, spectrum, target_spectrum, obs_conf events["reco_source_fov_offset"] = calculate_source_fov_offset( events, prefix="reco" ) - # Gamma source is assumed to be pointlike + # TODO: Honestly not sure why this integral is needed, nor what + # are correct bounds if kind == "gamma": - spectrum = spectrum.integrate_cone(0 * u.deg, self.ON_radius * u.deg) + spectrum = spectrum.integrate_cone( + self.bins.fov_offset_min * u.deg, self.bins.fov_offset_max * u.deg + ) events["weight"] = calculate_event_weights( events["true_energy"], target_spectrum=target_spectrum, @@ -161,7 +164,6 @@ def get_metadata(self, loader): f"Unsupported: '{itm}' differs across simulation runs" ) - assert sim["max_viewcone_radius"].std() == 0 sim_info = SimulatedEventsInfo( n_showers=show["n_entries"].sum(), energy_min=sim["energy_range_min"].quantity[0], @@ -199,21 +201,44 @@ def load_preselected_events(self): self.sim_info = sim_info self.spectrum = spectrum bits = [header] + n_raw_events = 0 for start, stop, events in load.read_subarray_events_chunked( self.chunk_size ): - selected = self.epp.normalise_column_names(events) - selected = selected[self.epp.get_table_mask(selected)] + selected = events[self.epp.get_table_mask(events)] + selected = self.epp.normalise_column_names(selected) selected = self.make_derived_columns( kind, selected, spectrum, target_spectrum, obs_conf ) bits.append(selected) + n_raw_events += len(events) table = vstack(bits, join_type="exact") reduced_events[kind] = table - - select_ON = reduced_events["gamma"]["theta"] <= self.ON_radius * u.deg - self.signal_events = reduced_events["gamma"][select_ON] + reduced_events[f"{kind}_count"] = n_raw_events + + self.log.debug( + "Loaded %d gammas, %d protons, %d electrons" + % ( + reduced_events["gamma_count"], + reduced_events["proton_count"], + reduced_events["electron_count"], + ) + ) + self.log.debug( + "Keeping %d gammas, %d protons, %d electrons" + % ( + len(reduced_events["gamma"]), + len(reduced_events["proton"]), + len(reduced_events["electron"]), + ) + ) + self.log.debug(self.epp.to_table()) + select_fov = ( + reduced_events["gamma"]["true_source_fov_offset"] + <= self.bins.fov_offset_max * u.deg + ) + self.signal_events = reduced_events["gamma"][select_fov] self.background_events = vstack( [reduced_events["proton"], reduced_events["electron"]] ) @@ -234,6 +259,10 @@ def setup(self): def start(self): self.load_preselected_events() + self.log.info( + "Optimising cuts using %d signal and %d background events" + % (len(self.signal_events), len(self.background_events)), + ) self.gh_cuts, self.theta_cuts_opt, self.sens2 = self.co.optimise_gh_cut( self.signal_events, self.background_events, @@ -286,16 +315,16 @@ def finish(self): fits.PrimaryHDU(), fits.BinTableHDU(self.sensitivity, name="SENSITIVITY"), fits.BinTableHDU(self.sens2, name="SENSITIVITY_STEP_2"), - fits.BinTableHDU(self.theta_cuts, name="THETA_CUTS"), fits.BinTableHDU(self.theta_cuts_opt, name="THETA_CUTS_OPT"), fits.BinTableHDU(self.gh_cuts, name="GH_CUTS"), ] for label, mask in masks.items(): - effective_area = effective_area_per_energy( + effective_area = effective_area_per_energy_and_fov( self.signal_events[mask], self.sim_info, true_energy_bins=self.true_energy_bins, + fov_offset_bins=self.fov_offset_bins, ) self.log.debug(self.true_energy_bins) self.log.debug(self.fov_offset_bins) From 2d5ce2251bbaf82d8434bd60cdb7dbca59f03306 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Tue, 19 Sep 2023 16:33:19 +0200 Subject: [PATCH 29/37] Reworked how initial theta cuts are calculated, changed some logging printouts --- ctapipe/irf/irf_classes.py | 26 +++++++++++++++++--------- ctapipe/tools/make_irf.py | 6 +++--- 2 files changed, 20 insertions(+), 12 deletions(-) diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index bf2409135ce..893dfea9e0d 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -74,18 +74,26 @@ def reco_energy_bins(self): 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) - ) - self.log.info( - f"Using fixed G/H cut of {INITIAL_GH_CUT} to calculate theta cuts" + initial_gh_cuts = calculate_percentile_cut( + signal["gh_score"], + signal["reco_energy"], + bins=self.reco_energy_bins(), + fill_value=0.0, + percentile=100 * (1 - self.initial_gh_cut_efficency), + min_events=25, + smoothing=1, ) - mask_theta_cuts = signal["gh_score"] >= INITIAL_GH_CUT + initial_gh_mask = evaluate_binned_cut( + signal["gh_score"], + signal["reco_energy"], + initial_gh_cuts, + op=operator.gt, + ) theta_cuts = calculate_percentile_cut( - signal["theta"][mask_theta_cuts], - signal["reco_energy"][mask_theta_cuts], + signal["theta"][initial_gh_mask], + signal["reco_energy"][initial_gh_mask], bins=self.reco_energy_bins(), min_value=self.theta_min_angle * u.deg, max_value=self.theta_max_angle * u.deg, @@ -325,7 +333,7 @@ class DataBinning(Component): fov_offset_max = Float( help="Maximum value for FoV offset bins in degrees", - default_value=2.0, + default_value=5.0, ).tag(config=True) fov_offset_n_edges = Integer( diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index b02d0601349..f88f177c235 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -233,11 +233,11 @@ def load_preselected_events(self): len(reduced_events["electron"]), ) ) - self.log.debug(self.epp.to_table()) select_fov = ( reduced_events["gamma"]["true_source_fov_offset"] <= self.bins.fov_offset_max * u.deg ) + # TODO: verify that this fov cut on only gamma is ok self.signal_events = reduced_events["gamma"][select_fov] self.background_events = vstack( [reduced_events["proton"], reduced_events["electron"]] @@ -319,6 +319,8 @@ def finish(self): fits.BinTableHDU(self.gh_cuts, name="GH_CUTS"), ] + self.log.debug("True Energy bins", self.true_energy_bins) + self.log.debug("FoV offset bins", self.fov_offset_bins) for label, mask in masks.items(): effective_area = effective_area_per_energy_and_fov( self.signal_events[mask], @@ -326,8 +328,6 @@ def finish(self): true_energy_bins=self.true_energy_bins, fov_offset_bins=self.fov_offset_bins, ) - self.log.debug(self.true_energy_bins) - self.log.debug(self.fov_offset_bins) hdus.append( create_aeff2d_hdu( effective_area[..., np.newaxis], # +1 dimension for FOV offset From 9e8b68817fa6e8931dc2910aea4effa708789b4e Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Wed, 20 Sep 2023 13:39:57 +0200 Subject: [PATCH 30/37] Update conf.py copy Max's fix to the doc config to ignore traitlets things. --- docs/conf.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/conf.py b/docs/conf.py index b848261c4d8..e6edb74981d 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -107,7 +107,13 @@ def setup(app): ("py:class", "traitlets.traitlets.Int"), ("py:class", "traitlets.config.application.Application"), ("py:class", "traitlets.utils.sentinel.Sentinel"), + ("py:class", "traitlets.traitlets.T"), + ("py:class", "re.Pattern[t.Any]"), + ("py:class", "Sentinel"), + ("py:class", "ObserveHandler"), ("py:class", "traitlets.traitlets.ObserveHandler"), + ("py:obj", "traitlets.traitlets.G"), + ("py:obj", "traitlets.traitlets.S"), ("py:obj", "traitlets.config.boolean_flag"), ("py:obj", "traitlets.TraitError"), ("py:obj", "-v"), # fix for wrong syntax in a traitlets docstring From c62a9ef7ee74ab3b719306dd10c6133ace791973 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Wed, 20 Sep 2023 16:07:44 +0200 Subject: [PATCH 31/37] Remove blankspace, maybe that fixes things? --- docs/conf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/conf.py b/docs/conf.py index 398410d9065..790d74e5cc2 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -133,7 +133,7 @@ def setup(app): ("py:class", "traitlets.traitlets.T"), ("py:class", "re.Pattern[t.Any]"), ("py:class", "Sentinel"), - ("py:class", "ObserveHandler"), + ("py:class", "ObserveHandler"), ("py:class", "traitlets.traitlets.ObserveHandler"), ("py:obj", "traitlets.traitlets.G"), ("py:obj", "traitlets.traitlets.S"), From abcbb38c77cdd6f180d4ef8f4bc08970e4dcde4f Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Tue, 26 Sep 2023 15:28:44 +0200 Subject: [PATCH 32/37] Moved PSF related quantities into its own component --- ctapipe/irf/__init__.py | 9 +++- ctapipe/irf/irf_classes.py | 91 ++++++++++++++++++++++++-------------- ctapipe/tools/make_irf.py | 11 ++++- 3 files changed, 76 insertions(+), 35 deletions(-) diff --git a/ctapipe/irf/__init__.py b/ctapipe/irf/__init__.py index 825fcf698cd..e637f2e47e6 100644 --- a/ctapipe/irf/__init__.py +++ b/ctapipe/irf/__init__.py @@ -3,6 +3,13 @@ DataBinning, EventPreProcessor, OutputEnergyBinning, + PointSpreadFunction, ) -__all__ = ["CutOptimizer", "DataBinning", "OutputEnergyBinning", "EventPreProcessor"] +__all__ = [ + "CutOptimizer", + "DataBinning", + "OutputEnergyBinning", + "EventPreProcessor", + "PointSpreadFunction", +] diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index 893dfea9e0d..b15482b3660 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -45,23 +45,6 @@ class CutOptimizer(Component): 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. @@ -73,7 +56,9 @@ def reco_energy_bins(self): ) return reco_energy - def optimise_gh_cut(self, signal, background, alpha, max_bg_radius): + def optimise_gh_cut( + self, signal, background, alpha, min_fov_radius, max_fov_radius, psf + ): initial_gh_cuts = calculate_percentile_cut( signal["gh_score"], signal["reco_energy"], @@ -91,15 +76,10 @@ def optimise_gh_cut(self, signal, background, alpha, max_bg_radius): op=operator.gt, ) - theta_cuts = calculate_percentile_cut( + theta_cuts = psf.calculate_theta_cuts( signal["theta"][initial_gh_mask], signal["reco_energy"][initial_gh_mask], - 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, + self.reco_energy_bins(), ) self.log.info("Optimizing G/H separation cut for best sensitivity") @@ -117,28 +97,75 @@ def optimise_gh_cut(self, signal, background, alpha, max_bg_radius): op=operator.ge, theta_cuts=theta_cuts, alpha=alpha, - fov_offset_max=max_bg_radius * u.deg, + fov_offset_max=max_fov_radius * u.deg, + fov_offset_min=min_fov_radius * u.deg, ) # 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 ) + self.log.info("Recalculating theta cut for optimized GH Cuts") - theta_cuts = calculate_percentile_cut( + theta_cuts = psf.calculate_theta_cuts( 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, + ) + + return gh_cuts, theta_cuts, sens2 + + +class PointSpreadFunction(Component): + theta_min_angle = Float( + default_value=-1, help="Smallest angular cut value allowed (-1 means no cut)" + ).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) + + theta_smoothing = Float( + default_value=-1, + help="When given, the width (in units of bins) of gaussian smoothing applied (-1)", + ).tag(config=True) + + target_percentile = Float( + default_value=68, + help="Percent of events in each energy bin keep after the theta cut", + ).tag(config=True) + + def calculate_theta_cuts(self, theta, reco_energy, energy_bins): + theta_min_angle = ( + None if self.theta_min_angle < 0 else self.theta_min_angle * u.deg + ) + theta_max_angle = ( + None if self.theta_max_angle < 0 else self.theta_max_angle * u.deg + ) + theta_smoothing = None if self.theta_smoothing < 0 else self.theta_smoothing + + return calculate_percentile_cut( + theta, + reco_energy, + energy_bins, + min_value=theta_min_angle, + max_value=theta_max_angle, + smoothing=theta_smoothing, + percentile=self.target_percentile, fill_value=self.theta_fill_value * u.deg, min_events=self.theta_min_counts, ) - return gh_cuts, theta_cuts, sens2 class EventPreProcessor(QualityQuery): diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index f88f177c235..d9008fb22fd 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -36,7 +36,13 @@ from ..core import Provenance, Tool, traits from ..core.traits import Bool, Float, Integer, Unicode from ..io import TableLoader -from ..irf import CutOptimizer, DataBinning, EventPreProcessor, OutputEnergyBinning +from ..irf import ( + CutOptimizer, + DataBinning, + EventPreProcessor, + OutputEnergyBinning, + PointSpreadFunction, +) class Spectra(Enum): @@ -244,10 +250,11 @@ def load_preselected_events(self): ) def setup(self): + self.epp = EventPreProcessor(parent=self) self.co = CutOptimizer(parent=self) + self.psf = PointSpreadFunction(parent=self) self.e_bins = OutputEnergyBinning(parent=self) self.bins = DataBinning(parent=self) - self.epp = EventPreProcessor(parent=self) self.reco_energy_bins = self.e_bins.reco_energy_bins() self.true_energy_bins = self.e_bins.true_energy_bins() From cc9c5d4668e40a711f0c9eed4fe91f6e35dcdcdc Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Tue, 26 Sep 2023 16:58:25 +0200 Subject: [PATCH 33/37] Renamed PSF component to ThetaCutsCalculator, other small refactors --- ctapipe/irf/__init__.py | 4 ++-- ctapipe/irf/irf_classes.py | 16 ++++++++-------- ctapipe/tools/make_irf.py | 8 +++++--- 3 files changed, 15 insertions(+), 13 deletions(-) diff --git a/ctapipe/irf/__init__.py b/ctapipe/irf/__init__.py index e637f2e47e6..e6d04e2cff0 100644 --- a/ctapipe/irf/__init__.py +++ b/ctapipe/irf/__init__.py @@ -3,7 +3,7 @@ DataBinning, EventPreProcessor, OutputEnergyBinning, - PointSpreadFunction, + ThetaCutsCalculator, ) __all__ = [ @@ -11,5 +11,5 @@ "DataBinning", "OutputEnergyBinning", "EventPreProcessor", - "PointSpreadFunction", + "ThetaCutsCalculator", ] diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index b15482b3660..0981fe625bf 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -17,6 +17,10 @@ class CutOptimizer(Component): """Performs cut optimisation""" + initial_gh_cut_efficency = Float( + default_value=0.4, help="Start value of gamma efficiency before optimisation" + ).tag(config=True) + max_gh_cut_efficiency = Float( default_value=0.8, help="Maximum gamma efficiency requested" ).tag(config=True) @@ -26,10 +30,6 @@ class CutOptimizer(Component): 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) - reco_energy_min = Float( help="Minimum value for Reco Energy bins in TeV units", default_value=0.005, @@ -57,7 +57,7 @@ def reco_energy_bins(self): return reco_energy def optimise_gh_cut( - self, signal, background, alpha, min_fov_radius, max_fov_radius, psf + self, signal, background, alpha, min_fov_radius, max_fov_radius, theta ): initial_gh_cuts = calculate_percentile_cut( signal["gh_score"], @@ -76,7 +76,7 @@ def optimise_gh_cut( op=operator.gt, ) - theta_cuts = psf.calculate_theta_cuts( + theta_cuts = theta.calculate_theta_cuts( signal["theta"][initial_gh_mask], signal["reco_energy"][initial_gh_mask], self.reco_energy_bins(), @@ -109,7 +109,7 @@ def optimise_gh_cut( ) self.log.info("Recalculating theta cut for optimized GH Cuts") - theta_cuts = psf.calculate_theta_cuts( + theta_cuts = theta.calculate_theta_cuts( signal[signal["selected_gh"]]["theta"], signal[signal["selected_gh"]]["reco_energy"], self.reco_energy_bins(), @@ -118,7 +118,7 @@ def optimise_gh_cut( return gh_cuts, theta_cuts, sens2 -class PointSpreadFunction(Component): +class ThetaCutsCalculator(Component): theta_min_angle = Float( default_value=-1, help="Smallest angular cut value allowed (-1 means no cut)" ).tag(config=True) diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index d9008fb22fd..cd34ee5cf51 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -41,7 +41,7 @@ DataBinning, EventPreProcessor, OutputEnergyBinning, - PointSpreadFunction, + ThetaCutsCalculator, ) @@ -252,7 +252,7 @@ def load_preselected_events(self): def setup(self): self.epp = EventPreProcessor(parent=self) self.co = CutOptimizer(parent=self) - self.psf = PointSpreadFunction(parent=self) + self.theta = ThetaCutsCalculator(parent=self) self.e_bins = OutputEnergyBinning(parent=self) self.bins = DataBinning(parent=self) @@ -274,7 +274,9 @@ def start(self): self.signal_events, self.background_events, self.alpha, - self.max_bg_radius, + self.bins.fov_offset_min, + self.bins.fov_offset_max, + self.theta, ) self.signal_events["selected_theta"] = evaluate_binned_cut( From 412900c24392104ba8f65e4e58bb3b7677e1a53e Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Thu, 5 Oct 2023 11:05:50 +0200 Subject: [PATCH 34/37] Update to support newest pyirf version --- ctapipe/tools/make_irf.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index cd34ee5cf51..5ddae092c6c 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -176,7 +176,8 @@ def get_metadata(self, loader): energy_max=sim["energy_range_max"].quantity[0], max_impact=sim["max_scatter_range"].quantity[0], spectral_index=sim["spectral_index"][0], - viewcone=sim["max_viewcone_radius"].quantity[0], + viewcone_max=sim["max_viewcone_radius"].quantity[0], + viewcone_min=sim["min_viewcone_radius"].quantity[0], ) return ( From 57523e65dc5446e3912a3a73d739b794d9cf1ae7 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Thu, 5 Oct 2023 11:06:33 +0200 Subject: [PATCH 35/37] Fix logging error --- ctapipe/tools/make_irf.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index 5ddae092c6c..fedde8032dd 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -329,8 +329,8 @@ def finish(self): fits.BinTableHDU(self.gh_cuts, name="GH_CUTS"), ] - self.log.debug("True Energy bins", self.true_energy_bins) - self.log.debug("FoV offset bins", self.fov_offset_bins) + self.log.debug(f"True Energy bins: {str(self.true_energy_bins.value)}") + self.log.debug(f"FoV offset bins: {str(self.fov_offset_bins.value)}") for label, mask in masks.items(): effective_area = effective_area_per_energy_and_fov( self.signal_events[mask], From 09d9f442f76218041d9559db175b8e14b42c9391 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Wed, 18 Oct 2023 18:50:42 +0200 Subject: [PATCH 36/37] Use consistent offset binning --- ctapipe/irf/irf_classes.py | 30 ------------------------------ ctapipe/tools/make_irf.py | 5 ++--- 2 files changed, 2 insertions(+), 33 deletions(-) diff --git a/ctapipe/irf/irf_classes.py b/ctapipe/irf/irf_classes.py index 0981fe625bf..bd7d99ff4f5 100644 --- a/ctapipe/irf/irf_classes.py +++ b/ctapipe/irf/irf_classes.py @@ -368,21 +368,6 @@ class DataBinning(Component): default_value=2, ).tag(config=True) - bkg_fov_offset_min = Float( - help="Minimum value for FoV offset bins for Background IRF", - default_value=0, - ).tag(config=True) - - bkg_fov_offset_max = Float( - help="Maximum value for FoV offset bins for Background IRF", - default_value=10, - ).tag(config=True) - - bkg_fov_offset_n_edges = Integer( - help="Number of edges for FoV offset bins for Background IRF", - default_value=21, - ).tag(config=True) - source_offset_min = Float( help="Minimum value for Source offset for PSF IRF", default_value=0, @@ -412,21 +397,6 @@ def fov_offset_bins(self): ) return fov_offset - def bkg_fov_offset_bins(self): - """ - Creates bins for FoV offset for Background IRF, - Using the same binning as in pyirf example. - """ - background_offset = ( - np.linspace( - self.bkg_fov_offset_min, - self.bkg_fov_offset_max, - self.bkg_fov_offset_n_edges, - ) - * u.deg - ) - return background_offset - def source_offset_bins(self): """ Creates bins for source offset for generating PSF IRF. diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index fedde8032dd..e5cb02bfe9b 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -263,7 +263,6 @@ def setup(self): self.source_offset_bins = self.bins.source_offset_bins() self.fov_offset_bins = self.bins.fov_offset_bins() - self.bkg_fov_offset_bins = self.bins.bkg_fov_offset_bins() def start(self): self.load_preselected_events() @@ -382,14 +381,14 @@ def finish(self): background_rate = background_2d( self.background_events[self.background_events["selected_gh"]], self.reco_energy_bins, - fov_offset_bins=self.bkg_fov_offset_bins, + fov_offset_bins=self.fov_offset_bins, t_obs=self.obs_time * u.Unit(self.obs_time_unit), ) hdus.append( create_background_2d_hdu( background_rate, self.reco_energy_bins, - fov_offset_bins=self.bkg_fov_offset_bins, + fov_offset_bins=self.fov_offset_bins, ) ) From 483b40a4460d69af5d7f555616b59e69561dd0f5 Mon Sep 17 00:00:00 2001 From: Tomas Bylund Date: Mon, 23 Oct 2023 17:04:48 +0200 Subject: [PATCH 37/37] Add theta cut to the background, change logging statments --- ctapipe/tools/make_irf.py | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/ctapipe/tools/make_irf.py b/ctapipe/tools/make_irf.py index e5cb02bfe9b..1c168c4ac4b 100644 --- a/ctapipe/tools/make_irf.py +++ b/ctapipe/tools/make_irf.py @@ -288,7 +288,12 @@ def start(self): self.signal_events["selected"] = ( self.signal_events["selected_theta"] & self.signal_events["selected_gh"] ) - + self.background_events["selected_theta"] = evaluate_binned_cut( + self.background_events["theta"], + self.background_events["reco_energy"], + self.theta_cuts_opt, + operator.le, + ) # calculate sensitivity signal_hist = create_histogram_table( self.signal_events[self.signal_events["selected"]], @@ -328,8 +333,8 @@ def finish(self): fits.BinTableHDU(self.gh_cuts, name="GH_CUTS"), ] - self.log.debug(f"True Energy bins: {str(self.true_energy_bins.value)}") - self.log.debug(f"FoV offset bins: {str(self.fov_offset_bins.value)}") + self.log.debug("True Energy bins: %s" % str(self.true_energy_bins.value)) + self.log.debug("FoV offset bins: %s" % str(self.fov_offset_bins.value)) for label, mask in masks.items(): effective_area = effective_area_per_energy_and_fov( self.signal_events[mask], @@ -364,8 +369,9 @@ def finish(self): # current pipelines comparisons bias_resolution = energy_bias_resolution( self.signal_events[self.signal_events["selected"]], - self.reco_energy_bins, - energy_type="reco", + self.true_energy_bins, + bias_function=np.mean, + energy_type="true", ) hdus.append(fits.BinTableHDU(bias_resolution, name="ENERGY_BIAS_RESOLUTION")) @@ -378,8 +384,11 @@ def finish(self): ) hdus.append(fits.BinTableHDU(ang_res, name="ANGULAR_RESOLUTION")) + sel = self.background_events["selected_gh"] + self.log.debug("%d background events selected" % sel.sum()) + self.log.debug("%f obs time" % self.obs_time) background_rate = background_2d( - self.background_events[self.background_events["selected_gh"]], + self.background_events[sel], self.reco_energy_bins, fov_offset_bins=self.fov_offset_bins, t_obs=self.obs_time * u.Unit(self.obs_time_unit),