From ac0d71843984b0197537dd4d098bcbd5577f3ab3 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 20 Feb 2025 10:19:55 +0100 Subject: [PATCH 01/58] Add function to compute selection cuts --- src/ctapipe/irf/cuts.py | 48 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 src/ctapipe/irf/cuts.py diff --git a/src/ctapipe/irf/cuts.py b/src/ctapipe/irf/cuts.py new file mode 100644 index 00000000000..22c6519c288 --- /dev/null +++ b/src/ctapipe/irf/cuts.py @@ -0,0 +1,48 @@ +import operator + +from astropy.table import Table +from pyirf.cuts import evaluate_binned_cut + +from ctapipe.irf import OptimizationResult + + +def calculate_selections( + events: Table, cuts: OptimizationResult, apply_spatial_selection: bool = False +) -> Table: + """ + Add the selection columns to the events + + Parameters + ---------- + events: Table + The table containing the events on which selection need to be applied + cuts: OptimizationResult + The cuts that need to be applied on the events + apply_spatial_selection: bool + True if the theta cuts should be applied + + Returns + ------- + Table + events with selection columns added. + """ + + events["selected_gh"] = evaluate_binned_cut( + events["gh_score"], + events["reco_energy"], + cuts.gh_cuts, + operator.ge, + ) + + if apply_spatial_selection: + events["selected_theta"] = evaluate_binned_cut( + events["theta"], + events["reco_energy"], + cuts.spatial_selection_table, + operator.le, + ) + events["selected"] = events["selected_theta"] & events["selected_gh"] + else: + events["selected"] = events["gammas"]["selected_gh"] + + return events From 17f0ed0d1ef8a513cd06a1d8709faf127da25746 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 20 Feb 2025 12:55:01 +0100 Subject: [PATCH 02/58] Make a component for event selection --- src/ctapipe/irf/cuts.py | 121 +++++++++++++++++++++++++++++----------- 1 file changed, 88 insertions(+), 33 deletions(-) diff --git a/src/ctapipe/irf/cuts.py b/src/ctapipe/irf/cuts.py index 22c6519c288..f6b2b66476f 100644 --- a/src/ctapipe/irf/cuts.py +++ b/src/ctapipe/irf/cuts.py @@ -3,46 +3,101 @@ from astropy.table import Table from pyirf.cuts import evaluate_binned_cut +from ctapipe.core import Component, QualityQuery +from ctapipe.core.traits import List, Path, Tuple, Unicode from ctapipe.irf import OptimizationResult +__all__ = ["EventQualityQuery", "EventSelection"] -def calculate_selections( - events: Table, cuts: OptimizationResult, apply_spatial_selection: bool = False -) -> Table: + +class EventQualityQuery(QualityQuery): """ - Add the selection columns to the events - - Parameters - ---------- - events: Table - The table containing the events on which selection need to be applied - cuts: OptimizationResult - The cuts that need to be applied on the events - apply_spatial_selection: bool - True if the theta cuts should be applied - - Returns - ------- - Table - events with selection columns added. + Event pre-selection quality criteria for IRF computation with different defaults. """ - events["selected_gh"] = evaluate_binned_cut( - events["gh_score"], - events["reco_energy"], - cuts.gh_cuts, - operator.ge, - ) + quality_criteria = List( + Tuple(Unicode(), Unicode()), + default_value=[ + ( + "multiplicity 4", + "np.count_nonzero(HillasReconstructor_telescopes,axis=1) >= 4", + ), + ("valid classifier", "RandomForestClassifier_is_valid"), + ("valid geom reco", "HillasReconstructor_is_valid"), + ("valid energy reco", "RandomForestRegressor_is_valid"), + ], + help=QualityQuery.quality_criteria.help, + ).tag(config=True) + + +class EventSelection(Component): + """ + Event selection + """ + + quality_query = EventQualityQuery() + cuts_file = Path( + default_value=None, + allow_none=False, + directory_ok=False, + help="Path to the cuts file to apply to the observation.", + ).tag(config=True) + + classes = [EventQualityQuery] - if apply_spatial_selection: - events["selected_theta"] = evaluate_binned_cut( - events["theta"], + def __init__(self, **kwargs): + super().__init__(**kwargs) + self.cuts = OptimizationResult.read(self.cuts_file) + + def calculate_selection(self, events): + events = self.calculate_quality_selection(events) + events = self.calculate_gamma_selection(events) + events["selected"] = events["selected_quality"] & events["selected_gamma"] + return events + + def calculate_quality_selection(self, events): + events["selected_quality"] = self.quality_query.get_table_mask(events) + events["selected"] = events["selected_quality"] + return events + + def calculate_gamma_selection( + self, events: Table, apply_spatial_selection: bool = False + ) -> Table: + """ + Add the selection columns to the events + + Parameters + ---------- + events: Table + The table containing the events on which selection need to be applied + cuts: OptimizationResult + The cuts that need to be applied on the events + apply_spatial_selection: bool + True if the theta cuts should be applied + + Returns + ------- + Table + events with selection columns added. + """ + + events["selected_gh"] = evaluate_binned_cut( + events["gh_score"], events["reco_energy"], - cuts.spatial_selection_table, - operator.le, + self.cuts.gh_cuts, + operator.ge, ) - events["selected"] = events["selected_theta"] & events["selected_gh"] - else: - events["selected"] = events["gammas"]["selected_gh"] - return events + if apply_spatial_selection: + events["selected_theta"] = evaluate_binned_cut( + events["theta"], + events["reco_energy"], + self.cuts.spatial_selection_table, + operator.le, + ) + events["selected_gamma"] = events["selected_theta"] & events["selected_gh"] + else: + events["selected_gamma"] = events["selected_gh"] + + events["selected"] = events["selected_gamma"] + return events From dbce02b5a5c318016575bc150a70f1cf93dc480d Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 20 Feb 2025 13:18:18 +0100 Subject: [PATCH 03/58] Generalise behaviour --- src/ctapipe/irf/cuts.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/ctapipe/irf/cuts.py b/src/ctapipe/irf/cuts.py index f6b2b66476f..3c6394016f1 100644 --- a/src/ctapipe/irf/cuts.py +++ b/src/ctapipe/irf/cuts.py @@ -3,14 +3,14 @@ from astropy.table import Table from pyirf.cuts import evaluate_binned_cut -from ctapipe.core import Component, QualityQuery +from ctapipe.core import QualityQuery from ctapipe.core.traits import List, Path, Tuple, Unicode from ctapipe.irf import OptimizationResult -__all__ = ["EventQualityQuery", "EventSelection"] +__all__ = ["EventQualitySelection", "EventSelection"] -class EventQualityQuery(QualityQuery): +class EventQualitySelection(QualityQuery): """ Event pre-selection quality criteria for IRF computation with different defaults. """ @@ -29,13 +29,20 @@ class EventQualityQuery(QualityQuery): help=QualityQuery.quality_criteria.help, ).tag(config=True) + def calculate_selection(self, events): + return self.calculate_quality_selection(events) + + def calculate_quality_selection(self, events): + events["selected_quality"] = self.get_table_mask(events) + events["selected"] = events["selected_quality"] + return events + -class EventSelection(Component): +class EventSelection(EventQualitySelection): """ Event selection """ - quality_query = EventQualityQuery() cuts_file = Path( default_value=None, allow_none=False, @@ -43,8 +50,6 @@ class EventSelection(Component): help="Path to the cuts file to apply to the observation.", ).tag(config=True) - classes = [EventQualityQuery] - def __init__(self, **kwargs): super().__init__(**kwargs) self.cuts = OptimizationResult.read(self.cuts_file) @@ -55,11 +60,6 @@ def calculate_selection(self, events): events["selected"] = events["selected_quality"] & events["selected_gamma"] return events - def calculate_quality_selection(self, events): - events["selected_quality"] = self.quality_query.get_table_mask(events) - events["selected"] = events["selected_quality"] - return events - def calculate_gamma_selection( self, events: Table, apply_spatial_selection: bool = False ) -> Table: From 7606d2217a4deb0248c1d7e86eecaf6fdc9f5e44 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 20 Feb 2025 13:19:13 +0100 Subject: [PATCH 04/58] Update event loader in the view of using it for DL3 production --- src/ctapipe/irf/optimize.py | 8 +- src/ctapipe/irf/preprocessing.py | 89 +++++++------------ src/ctapipe/irf/tests/test_benchmarks.py | 10 +-- src/ctapipe/irf/tests/test_optimize.py | 10 +-- src/ctapipe/irf/tests/test_preprocessing.py | 8 +- src/ctapipe/tools/compute_irf.py | 11 ++- src/ctapipe/tools/create_dl3.py | 87 ++++++++++++++++++ src/ctapipe/tools/optimize_event_selection.py | 7 +- 8 files changed, 144 insertions(+), 86 deletions(-) create mode 100644 src/ctapipe/tools/create_dl3.py diff --git a/src/ctapipe/irf/optimize.py b/src/ctapipe/irf/optimize.py index fd8f9e9ad52..24666e4e67d 100644 --- a/src/ctapipe/irf/optimize.py +++ b/src/ctapipe/irf/optimize.py @@ -14,7 +14,7 @@ from ..core import Component, QualityQuery from ..core.traits import AstroQuantity, Float, Integer, Path from .binning import DefaultRecoEnergyBins, ResultValidRange -from .preprocessing import EventQualityQuery +from .cuts import EventQualitySelection __all__ = [ "CutOptimizerBase", @@ -168,7 +168,7 @@ def _check_events(self, events: dict[str, QTable]): def __call__( self, events: dict[str, QTable], - quality_query: EventQualityQuery, + quality_query: EventQualitySelection, clf_prefix: str, ) -> OptimizationResult: """ @@ -305,7 +305,7 @@ def __init__(self, config=None, parent=None, **kwargs): def __call__( self, events: dict[str, QTable], - quality_query: EventQualityQuery, + quality_query: EventQualitySelection, clf_prefix: str, ) -> OptimizationResult: self._check_events(events) @@ -393,7 +393,7 @@ def __init__(self, config=None, parent=None, **kwargs): def __call__( self, events: dict[str, QTable], - quality_query: EventQualityQuery, + quality_query: EventQualitySelection, clf_prefix: str, ) -> OptimizationResult: self._check_events(events) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 21653693793..caf4e861996 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -19,38 +19,19 @@ from ..compat import COPY_IF_NEEDED from ..containers import CoordinateFrameType from ..coordinates import NominalFrame -from ..core import Component, QualityQuery -from ..core.traits import List, Tuple, Unicode +from ..core import Component +from ..core.traits import Unicode from ..io import TableLoader +from .cuts import EventQualitySelection from .spectra import SPECTRA, Spectra -__all__ = ["EventLoader", "EventPreprocessor", "EventQualityQuery"] - - -class EventQualityQuery(QualityQuery): - """ - Event pre-selection quality criteria for IRF computation with different defaults. - """ - - quality_criteria = List( - Tuple(Unicode(), Unicode()), - default_value=[ - ( - "multiplicity 4", - "np.count_nonzero(HillasReconstructor_telescopes,axis=1) >= 4", - ), - ("valid classifier", "RandomForestClassifier_is_valid"), - ("valid geom reco", "HillasReconstructor_is_valid"), - ("valid energy reco", "RandomForestRegressor_is_valid"), - ], - help=QualityQuery.quality_criteria.help, - ).tag(config=True) +__all__ = ["EventLoader", "EventPreprocessor"] class EventPreprocessor(Component): """Defines pre-selection cuts and the necessary renaming of columns.""" - classes = [EventQualityQuery] + classes = [EventQualitySelection] energy_reconstructor = Unicode( default_value="RandomForestRegressor", @@ -69,7 +50,7 @@ class EventPreprocessor(Component): def __init__(self, config=None, parent=None, **kwargs): super().__init__(config=config, parent=parent, **kwargs) - self.quality_query = EventQualityQuery(parent=self) + self.quality_query = EventQualitySelection(parent=self) def normalise_column_names(self, events: Table) -> QTable: if events["subarray_pointing_lat"].std() > 1e-3: @@ -213,52 +194,48 @@ def __init__(self, file: Path, target_spectrum: Spectra, **kwargs): self.target_spectrum = SPECTRA[target_spectrum] self.file = file - def load_preselected_events( - self, chunk_size: int, obs_time: u.Quantity - ) -> tuple[QTable, int, dict]: + def load_preselected_events(self, chunk_size: int) -> tuple[QTable, int, dict]: opts = dict(dl2=True, simulated=True, observation_info=True) with TableLoader(self.file, parent=self, **opts) as load: header = self.epp.make_empty_table() - sim_info, spectrum = self.get_simulation_information(load, obs_time) - meta = {"sim_info": sim_info, "spectrum": spectrum} bits = [header] - n_raw_events = 0 for _, _, events in load.read_subarray_events_chunked(chunk_size, **opts): selected = events[self.epp.quality_query.get_table_mask(events)] selected = self.epp.normalise_column_names(selected) selected = self.make_derived_columns(selected) bits.append(selected) - n_raw_events += len(events) bits.append(header) # Putting it last ensures the correct metadata is used table = vstack(bits, join_type="exact", metadata_conflicts="silent") - return table, n_raw_events, meta + return table def get_simulation_information( - self, loader: TableLoader, obs_time: u.Quantity + self, obs_time: u.Quantity ) -> tuple[SimulatedEventsInfo, PowerLaw]: - sim = loader.read_simulation_configuration() - try: - show = loader.read_shower_distribution() - except NoSuchNodeError: - # Fall back to using the run header - show = Table([sim["n_showers"]], names=["n_entries"], dtype=[np.int64]) - - 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" - ) - - sim_info = SimulatedEventsInfo( - 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], - spectral_index=sim["spectral_index"][0], - viewcone_max=sim["max_viewcone_radius"].quantity[0], - viewcone_min=sim["min_viewcone_radius"].quantity[0], - ) + opts = dict(dl2=True, simulated=True, observation_info=True) + with TableLoader(self.file, parent=self, **opts) as load: + sim = load.read_simulation_configuration() + try: + show = load.read_shower_distribution() + except NoSuchNodeError: + # Fall back to using the run header + show = Table([sim["n_showers"]], names=["n_entries"], dtype=[np.int64]) + + 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" + ) + + sim_info = SimulatedEventsInfo( + 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], + spectral_index=sim["spectral_index"][0], + viewcone_max=sim["max_viewcone_radius"].quantity[0], + viewcone_min=sim["min_viewcone_radius"].quantity[0], + ) return sim_info, PowerLaw.from_simulation(sim_info, obstime=obs_time) diff --git a/src/ctapipe/irf/tests/test_benchmarks.py b/src/ctapipe/irf/tests/test_benchmarks.py index 8264003c9ee..076fe0d42fc 100644 --- a/src/ctapipe/irf/tests/test_benchmarks.py +++ b/src/ctapipe/irf/tests/test_benchmarks.py @@ -91,19 +91,13 @@ def test_make_2d_sensitivity( file=gamma_diffuse_full_reco_file, target_spectrum=Spectra.CRAB_HEGRA, ) - gamma_events, _, _ = gamma_loader.load_preselected_events( - chunk_size=10000, - obs_time=u.Quantity(50, u.h), - ) + gamma_events = gamma_loader.load_preselected_events(chunk_size=10000) proton_loader = EventLoader( config=irf_event_loader_test_config, file=proton_full_reco_file, target_spectrum=Spectra.IRFDOC_PROTON_SPECTRUM, ) - proton_events, _, _ = proton_loader.load_preselected_events( - chunk_size=10000, - obs_time=u.Quantity(50, u.h), - ) + proton_events = proton_loader.load_preselected_events(chunk_size=10000) sens_maker = Sensitivity2dMaker( fov_offset_n_bins=3, diff --git a/src/ctapipe/irf/tests/test_optimize.py b/src/ctapipe/irf/tests/test_optimize.py index 366c7ca334a..2bb42b2f5ca 100644 --- a/src/ctapipe/irf/tests/test_optimize.py +++ b/src/ctapipe/irf/tests/test_optimize.py @@ -94,19 +94,13 @@ def test_cut_optimizer( file=gamma_diffuse_full_reco_file, target_spectrum=Spectra.CRAB_HEGRA, ) - gamma_events, _, _ = gamma_loader.load_preselected_events( - chunk_size=10000, - obs_time=u.Quantity(50, u.h), - ) + gamma_events = gamma_loader.load_preselected_events(chunk_size=10000) proton_loader = EventLoader( config=irf_event_loader_test_config, file=proton_full_reco_file, target_spectrum=Spectra.IRFDOC_PROTON_SPECTRUM, ) - proton_events, _, _ = proton_loader.load_preselected_events( - chunk_size=10000, - obs_time=u.Quantity(50, u.h), - ) + proton_events = proton_loader.load_preselected_events(chunk_size=10000) optimizer = Optimizer() result = optimizer( diff --git a/src/ctapipe/irf/tests/test_preprocessing.py b/src/ctapipe/irf/tests/test_preprocessing.py index 71cedc22ade..0897e4ed872 100644 --- a/src/ctapipe/irf/tests/test_preprocessing.py +++ b/src/ctapipe/irf/tests/test_preprocessing.py @@ -71,10 +71,10 @@ def test_event_loader(gamma_diffuse_full_reco_file, irf_event_loader_test_config file=gamma_diffuse_full_reco_file, target_spectrum=Spectra.CRAB_HEGRA, ) - events, count, meta = loader.load_preselected_events( - chunk_size=10000, - obs_time=u.Quantity(50, u.h), - ) + events = loader.load_preselected_events(chunk_size=10000) + count = len(events) + sim_info, spectrum = loader.get_simulation_information(obs_time=u.Quantity(50, u.h)) + meta = {"sim_info": sim_info, "spectrum": spectrum} columns = [ "obs_id", diff --git a/src/ctapipe/tools/compute_irf.py b/src/ctapipe/tools/compute_irf.py index f9465293ec7..a8cc929edfa 100644 --- a/src/ctapipe/tools/compute_irf.py +++ b/src/ctapipe/tools/compute_irf.py @@ -30,13 +30,13 @@ EnergyBiasResolutionMakerBase, SensitivityMakerBase, ) +from ..irf.cuts import EventQualitySelection from ..irf.irfs import ( BackgroundRateMakerBase, EffectiveAreaMakerBase, EnergyDispersionMakerBase, PSFMakerBase, ) -from ..irf.preprocessing import EventQualityQuery __all__ = ["IrfTool"] @@ -487,7 +487,7 @@ def start(self): ], ) ) - loader.epp.quality_query = EventQualityQuery( + loader.epp.quality_query = EventQualitySelection( parent=loader, quality_criteria=self.opt_result.quality_query.quality_criteria, ) @@ -496,9 +496,12 @@ def start(self): "%s Quality criteria: %s" % (particle_type, loader.epp.quality_query.quality_criteria) ) - events, count, meta = loader.load_preselected_events( - self.chunk_size, self.obs_time + events = loader.load_preselected_events(self.chunk_size) + count = len(events) + sim_info, spectrum = loader.get_simulation_information( + obs_time=u.Quantity(50, u.h) ) + meta = {"sim_info": sim_info, "spectrum": spectrum} # Only calculate event weights if background or sensitivity should be calculated. if self.do_background: # Sensitivity is only calculated, if do_background is true diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py new file mode 100644 index 00000000000..c0ea3c20b78 --- /dev/null +++ b/src/ctapipe/tools/create_dl3.py @@ -0,0 +1,87 @@ +from ctapipe.core import Tool, traits +from ctapipe.core.traits import ( + Integer, +) + +from ..irf import ( + EventLoader, + OptimizationResult, +) + +__all__ = ["DL3Tool"] + + +class DL3Tool(Tool): + name = "mytool" + description = "do some things and stuff" + aliases = dict( + infile="AdvancedComponent.infile", + outfile="AdvancedComponent.outfile", + iterations="MyTool.iterations", + ) + + dl2_file = traits.Path( + default_value=None, + allow_none=False, + directory_ok=False, + help="DL2 input filename and path.", + ).tag(config=True) + + output_path = traits.Path( + default_value=None, + allow_none=False, + directory_ok=False, + help="Output file", + ).tag(config=True) + + cuts_file = traits.Path( + default_value=None, + allow_none=False, + directory_ok=False, + help="Path to the cuts file to apply to the observation.", + ).tag(config=True) + + irfs_file = traits.Path( + default_value=None, + allow_none=False, + directory_ok=False, + help="Path to the IRFs describing the observation", + ).tag(config=True) + + chunk_size = Integer( + default_value=100000, + allow_none=True, + help="How many subarray events to load at once while selecting.", + ).tag(config=True) + + # Which classes are registered for configuration + classes = [ + EventLoader, + ] + + def setup(self): + """ + Initialize components from config and load g/h (and theta) cuts. + """ + self.log.info("Loading cuts") + self.opt_result = OptimizationResult.read(self.cuts_file) + self.log.info("Loading events from DL2") + self.event_loader = EventLoader( + parent=self, + file=self.dl2_file, + ) + + def start(self): + pass + + def finish(self): + self.log.warning("Shutting down.") + + +def main(): + tool = DL3Tool() + tool.run() + + +if __name__ == "main": + main() diff --git a/src/ctapipe/tools/optimize_event_selection.py b/src/ctapipe/tools/optimize_event_selection.py index e2066fa77e1..3044d89ea5d 100644 --- a/src/ctapipe/tools/optimize_event_selection.py +++ b/src/ctapipe/tools/optimize_event_selection.py @@ -147,9 +147,12 @@ def start(self): """ reduced_events = dict() for particle_type, loader in self.event_loaders.items(): - events, count, meta = loader.load_preselected_events( - self.chunk_size, self.obs_time + events = loader.load_preselected_events(self.chunk_size) + count = len(events) + sim_info, spectrum = loader.get_simulation_information( + obs_time=u.Quantity(50, u.h) ) + meta = {"sim_info": sim_info, "spectrum": spectrum} if self.optimizer.needs_background: events = loader.make_event_weights( events, From 5111d1f9854ff5753c5a974e2c383ccbbd8a50ea Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 20 Feb 2025 14:50:59 +0100 Subject: [PATCH 05/58] change the output format of get simulation --- src/ctapipe/irf/preprocessing.py | 6 +++++- src/ctapipe/irf/tests/test_preprocessing.py | 3 +-- src/ctapipe/tools/compute_irf.py | 5 +---- src/ctapipe/tools/optimize_event_selection.py | 5 +---- 4 files changed, 8 insertions(+), 11 deletions(-) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index caf4e861996..b27ddf3bbc3 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -237,7 +237,11 @@ def get_simulation_information( viewcone_min=sim["min_viewcone_radius"].quantity[0], ) - return sim_info, PowerLaw.from_simulation(sim_info, obstime=obs_time) + meta = { + "sim_info": sim_info, + "spectrum": PowerLaw.from_simulation(sim_info, obstime=obs_time), + } + return meta def make_derived_columns(self, events: QTable) -> QTable: events["weight"] = ( diff --git a/src/ctapipe/irf/tests/test_preprocessing.py b/src/ctapipe/irf/tests/test_preprocessing.py index 0897e4ed872..0244933a6e6 100644 --- a/src/ctapipe/irf/tests/test_preprocessing.py +++ b/src/ctapipe/irf/tests/test_preprocessing.py @@ -73,8 +73,7 @@ def test_event_loader(gamma_diffuse_full_reco_file, irf_event_loader_test_config ) events = loader.load_preselected_events(chunk_size=10000) count = len(events) - sim_info, spectrum = loader.get_simulation_information(obs_time=u.Quantity(50, u.h)) - meta = {"sim_info": sim_info, "spectrum": spectrum} + meta = loader.get_simulation_information(obs_time=u.Quantity(50, u.h)) columns = [ "obs_id", diff --git a/src/ctapipe/tools/compute_irf.py b/src/ctapipe/tools/compute_irf.py index a8cc929edfa..4f505aad506 100644 --- a/src/ctapipe/tools/compute_irf.py +++ b/src/ctapipe/tools/compute_irf.py @@ -498,10 +498,7 @@ def start(self): ) events = loader.load_preselected_events(self.chunk_size) count = len(events) - sim_info, spectrum = loader.get_simulation_information( - obs_time=u.Quantity(50, u.h) - ) - meta = {"sim_info": sim_info, "spectrum": spectrum} + meta = loader.get_simulation_information(obs_time=u.Quantity(50, u.h)) # Only calculate event weights if background or sensitivity should be calculated. if self.do_background: # Sensitivity is only calculated, if do_background is true diff --git a/src/ctapipe/tools/optimize_event_selection.py b/src/ctapipe/tools/optimize_event_selection.py index 3044d89ea5d..431c26eb337 100644 --- a/src/ctapipe/tools/optimize_event_selection.py +++ b/src/ctapipe/tools/optimize_event_selection.py @@ -149,10 +149,7 @@ def start(self): for particle_type, loader in self.event_loaders.items(): events = loader.load_preselected_events(self.chunk_size) count = len(events) - sim_info, spectrum = loader.get_simulation_information( - obs_time=u.Quantity(50, u.h) - ) - meta = {"sim_info": sim_info, "spectrum": spectrum} + meta = loader.get_simulation_information(obs_time=u.Quantity(50, u.h)) if self.optimizer.needs_background: events = loader.make_event_weights( events, From 47fc14d4ca27f304566a25a28c9f24e295e3073d Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 20 Feb 2025 16:04:32 +0100 Subject: [PATCH 06/58] Fix cicular depedency --- src/ctapipe/irf/cuts.py | 6 +- src/ctapipe/irf/optimize/__init__.py | 15 ++ .../{optimize.py => optimize/algorithm.py} | 131 +----------------- src/ctapipe/irf/optimize/results.py | 127 +++++++++++++++++ .../{tests => optimize/test}/test_optimize.py | 2 +- src/ctapipe/tools/optimize_event_selection.py | 3 +- 6 files changed, 154 insertions(+), 130 deletions(-) create mode 100644 src/ctapipe/irf/optimize/__init__.py rename src/ctapipe/irf/{optimize.py => optimize/algorithm.py} (69%) create mode 100644 src/ctapipe/irf/optimize/results.py rename src/ctapipe/irf/{tests => optimize/test}/test_optimize.py (98%) diff --git a/src/ctapipe/irf/cuts.py b/src/ctapipe/irf/cuts.py index 3c6394016f1..610c21755c3 100644 --- a/src/ctapipe/irf/cuts.py +++ b/src/ctapipe/irf/cuts.py @@ -3,9 +3,9 @@ from astropy.table import Table from pyirf.cuts import evaluate_binned_cut -from ctapipe.core import QualityQuery -from ctapipe.core.traits import List, Path, Tuple, Unicode -from ctapipe.irf import OptimizationResult +from ..core import QualityQuery +from ..core.traits import List, Path, Tuple, Unicode +from .optimize.results import OptimizationResult __all__ = ["EventQualitySelection", "EventSelection"] diff --git a/src/ctapipe/irf/optimize/__init__.py b/src/ctapipe/irf/optimize/__init__.py new file mode 100644 index 00000000000..e1ccdd75b83 --- /dev/null +++ b/src/ctapipe/irf/optimize/__init__.py @@ -0,0 +1,15 @@ +from .algorithm import ( + GhPercentileCutCalculator, + PercentileCuts, + PointSourceSensitivityOptimizer, + ThetaPercentileCutCalculator, +) +from .results import OptimizationResult + +__all__ = [ + "OptimizationResult", + "GhPercentileCutCalculator", + "PercentileCuts", + "PointSourceSensitivityOptimizer", + "ThetaPercentileCutCalculator", +] diff --git a/src/ctapipe/irf/optimize.py b/src/ctapipe/irf/optimize/algorithm.py similarity index 69% rename from src/ctapipe/irf/optimize.py rename to src/ctapipe/irf/optimize/algorithm.py index 24666e4e67d..c4fc6a7aa2a 100644 --- a/src/ctapipe/irf/optimize.py +++ b/src/ctapipe/irf/optimize/algorithm.py @@ -2,147 +2,28 @@ import operator from abc import abstractmethod -from collections.abc import Sequence import astropy.units as u import numpy as np -from astropy.io import fits -from astropy.table import QTable, Table +from astropy.table import QTable from pyirf.cut_optimization import optimize_gh_cut from pyirf.cuts import calculate_percentile_cut, evaluate_binned_cut -from ..core import Component, QualityQuery -from ..core.traits import AstroQuantity, Float, Integer, Path -from .binning import DefaultRecoEnergyBins, ResultValidRange -from .cuts import EventQualitySelection +from ...core import Component +from ...core.traits import AstroQuantity, Float, Integer +from ..binning import DefaultRecoEnergyBins +from ..cuts import EventQualitySelection +from .results import OptimizationResult __all__ = [ "CutOptimizerBase", "GhPercentileCutCalculator", - "OptimizationResult", "PercentileCuts", "PointSourceSensitivityOptimizer", "ThetaPercentileCutCalculator", ] -class OptimizationResult: - """Result of an optimization of G/H and theta cuts or only G/H cuts.""" - - def __init__( - self, - valid_energy_min: u.Quantity, - valid_energy_max: u.Quantity, - valid_offset_min: u.Quantity, - valid_offset_max: u.Quantity, - gh_cuts: QTable, - clf_prefix: str, - spatial_selection_table: QTable | None = None, - quality_query: QualityQuery | Sequence | None = None, - ) -> None: - if quality_query: - if isinstance(quality_query, QualityQuery): - if len(quality_query.quality_criteria) == 0: - quality_query.quality_criteria = [ - (" ", " ") - ] # Ensures table serialises properly - - self.quality_query = quality_query - elif isinstance(quality_query, list): - self.quality_query = QualityQuery(quality_criteria=quality_query) - else: - self.quality_query = QualityQuery(quality_criteria=list(quality_query)) - else: - self.quality_query = QualityQuery(quality_criteria=[(" ", " ")]) - - self.valid_energy = ResultValidRange(min=valid_energy_min, max=valid_energy_max) - self.valid_offset = ResultValidRange(min=valid_offset_min, max=valid_offset_max) - self.gh_cuts = gh_cuts - self.clf_prefix = clf_prefix - self.spatial_selection_table = spatial_selection_table - - def __repr__(self): - if self.spatial_selection_table is not None: - return ( - f"" - ) - else: - return ( - f"" - ) - - def write(self, output_name: Path | str, overwrite: bool = False) -> None: - """Write an ``OptimizationResult`` to a file in FITS format.""" - - cut_expr_tab = Table( - rows=self.quality_query.quality_criteria, - names=["name", "cut_expr"], - dtype=[np.str_, np.str_], - ) - cut_expr_tab.meta["EXTNAME"] = "QUALITY_CUTS_EXPR" - - self.gh_cuts.meta["EXTNAME"] = "GH_CUTS" - self.gh_cuts.meta["CLFNAME"] = self.clf_prefix - - energy_lim_tab = QTable( - rows=[[self.valid_energy.min, self.valid_energy.max]], - names=["energy_min", "energy_max"], - ) - energy_lim_tab.meta["EXTNAME"] = "VALID_ENERGY" - - offset_lim_tab = QTable( - rows=[[self.valid_offset.min, self.valid_offset.max]], - names=["offset_min", "offset_max"], - ) - offset_lim_tab.meta["EXTNAME"] = "VALID_OFFSET" - - results = [cut_expr_tab, self.gh_cuts, energy_lim_tab, offset_lim_tab] - - if self.spatial_selection_table is not None: - self.spatial_selection_table.meta["EXTNAME"] = "RAD_MAX" - results.append(self.spatial_selection_table) - - # Overwrite if needed and allowed - results[0].write(output_name, format="fits", overwrite=overwrite) - - for table in results[1:]: - table.write(output_name, format="fits", append=True) - - @classmethod - def read(cls, file_name): - """Read an ``OptimizationResult`` from a file in FITS format.""" - - with fits.open(file_name) as hdul: - cut_expr_tab = Table.read(hdul[1]) - cut_expr_lst = [(name, expr) for name, expr in cut_expr_tab.iterrows()] - if (" ", " ") in cut_expr_lst: - cut_expr_lst.remove((" ", " ")) - - quality_query = QualityQuery(quality_criteria=cut_expr_lst) - gh_cuts = QTable.read(hdul[2]) - valid_energy = QTable.read(hdul[3]) - valid_offset = QTable.read(hdul[4]) - spatial_selection_table = QTable.read(hdul[5]) if len(hdul) > 5 else None - - return cls( - quality_query=quality_query, - valid_energy_min=valid_energy["energy_min"], - valid_energy_max=valid_energy["energy_max"], - valid_offset_min=valid_offset["offset_min"], - valid_offset_max=valid_offset["offset_max"], - gh_cuts=gh_cuts, - clf_prefix=gh_cuts.meta["CLFNAME"], - spatial_selection_table=spatial_selection_table, - ) - - class CutOptimizerBase(DefaultRecoEnergyBins): """Base class for cut optimization algorithms.""" diff --git a/src/ctapipe/irf/optimize/results.py b/src/ctapipe/irf/optimize/results.py new file mode 100644 index 00000000000..19de853e085 --- /dev/null +++ b/src/ctapipe/irf/optimize/results.py @@ -0,0 +1,127 @@ +from typing import Sequence + +import numpy as np +from astropy import units as u +from astropy.io import fits +from astropy.table import QTable, Table + +from ...core import QualityQuery +from ...core.traits import Path +from ...irf import ResultValidRange + + +class OptimizationResult: + """Result of an optimization of G/H and theta cuts or only G/H cuts.""" + + def __init__( + self, + valid_energy_min: u.Quantity, + valid_energy_max: u.Quantity, + valid_offset_min: u.Quantity, + valid_offset_max: u.Quantity, + gh_cuts: QTable, + clf_prefix: str, + spatial_selection_table: QTable | None = None, + quality_query: QualityQuery | Sequence | None = None, + ) -> None: + if quality_query: + if isinstance(quality_query, QualityQuery): + if len(quality_query.quality_criteria) == 0: + quality_query.quality_criteria = [ + (" ", " ") + ] # Ensures table serialises properly + + self.quality_query = quality_query + elif isinstance(quality_query, list): + self.quality_query = QualityQuery(quality_criteria=quality_query) + else: + self.quality_query = QualityQuery(quality_criteria=list(quality_query)) + else: + self.quality_query = QualityQuery(quality_criteria=[(" ", " ")]) + + self.valid_energy = ResultValidRange(min=valid_energy_min, max=valid_energy_max) + self.valid_offset = ResultValidRange(min=valid_offset_min, max=valid_offset_max) + self.gh_cuts = gh_cuts + self.clf_prefix = clf_prefix + self.spatial_selection_table = spatial_selection_table + + def __repr__(self): + if self.spatial_selection_table is not None: + return ( + f"" + ) + else: + return ( + f"" + ) + + def write(self, output_name: Path | str, overwrite: bool = False) -> None: + """Write an ``OptimizationResult`` to a file in FITS format.""" + + cut_expr_tab = Table( + rows=self.quality_query.quality_criteria, + names=["name", "cut_expr"], + dtype=[np.str_, np.str_], + ) + cut_expr_tab.meta["EXTNAME"] = "QUALITY_CUTS_EXPR" + + self.gh_cuts.meta["EXTNAME"] = "GH_CUTS" + self.gh_cuts.meta["CLFNAME"] = self.clf_prefix + + energy_lim_tab = QTable( + rows=[[self.valid_energy.min, self.valid_energy.max]], + names=["energy_min", "energy_max"], + ) + energy_lim_tab.meta["EXTNAME"] = "VALID_ENERGY" + + offset_lim_tab = QTable( + rows=[[self.valid_offset.min, self.valid_offset.max]], + names=["offset_min", "offset_max"], + ) + offset_lim_tab.meta["EXTNAME"] = "VALID_OFFSET" + + results = [cut_expr_tab, self.gh_cuts, energy_lim_tab, offset_lim_tab] + + if self.spatial_selection_table is not None: + self.spatial_selection_table.meta["EXTNAME"] = "RAD_MAX" + results.append(self.spatial_selection_table) + + # Overwrite if needed and allowed + results[0].write(output_name, format="fits", overwrite=overwrite) + + for table in results[1:]: + table.write(output_name, format="fits", append=True) + + @classmethod + def read(cls, file_name): + """Read an ``OptimizationResult`` from a file in FITS format.""" + + with fits.open(file_name) as hdul: + cut_expr_tab = Table.read(hdul[1]) + cut_expr_lst = [(name, expr) for name, expr in cut_expr_tab.iterrows()] + if (" ", " ") in cut_expr_lst: + cut_expr_lst.remove((" ", " ")) + + quality_query = QualityQuery(quality_criteria=cut_expr_lst) + gh_cuts = QTable.read(hdul[2]) + valid_energy = QTable.read(hdul[3]) + valid_offset = QTable.read(hdul[4]) + spatial_selection_table = QTable.read(hdul[5]) if len(hdul) > 5 else None + + return cls( + quality_query=quality_query, + valid_energy_min=valid_energy["energy_min"], + valid_energy_max=valid_energy["energy_max"], + valid_offset_min=valid_offset["offset_min"], + valid_offset_max=valid_offset["offset_max"], + gh_cuts=gh_cuts, + clf_prefix=gh_cuts.meta["CLFNAME"], + spatial_selection_table=spatial_selection_table, + ) diff --git a/src/ctapipe/irf/tests/test_optimize.py b/src/ctapipe/irf/optimize/test/test_optimize.py similarity index 98% rename from src/ctapipe/irf/tests/test_optimize.py rename to src/ctapipe/irf/optimize/test/test_optimize.py index 2bb42b2f5ca..bdd99ecf3f9 100644 --- a/src/ctapipe/irf/tests/test_optimize.py +++ b/src/ctapipe/irf/optimize/test/test_optimize.py @@ -4,7 +4,7 @@ from astropy.table import QTable from ctapipe.core import QualityQuery, non_abstract_children -from ctapipe.irf.optimize import CutOptimizerBase +from ctapipe.irf.optimize.algorithm import CutOptimizerBase def test_optimization_result(tmp_path, irf_event_loader_test_config): diff --git a/src/ctapipe/tools/optimize_event_selection.py b/src/ctapipe/tools/optimize_event_selection.py index 431c26eb337..6865a5db402 100644 --- a/src/ctapipe/tools/optimize_event_selection.py +++ b/src/ctapipe/tools/optimize_event_selection.py @@ -3,10 +3,11 @@ import astropy.units as u from astropy.table import vstack +from ctapipe.irf.optimize.algorithm import CutOptimizerBase + from ..core import Provenance, Tool, traits from ..core.traits import AstroQuantity, Integer, classes_with_traits from ..irf import EventLoader, Spectra -from ..irf.optimize import CutOptimizerBase __all__ = ["EventSelectionOptimizer"] From 1b2601710e2821403515fad115f7bbfb89afe070 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 20 Feb 2025 16:30:15 +0100 Subject: [PATCH 07/58] Update documentation of cuts --- src/ctapipe/irf/cuts.py | 61 +++++++++++++++++++++++++++++++++-------- 1 file changed, 50 insertions(+), 11 deletions(-) diff --git a/src/ctapipe/irf/cuts.py b/src/ctapipe/irf/cuts.py index 610c21755c3..59a1cf136a9 100644 --- a/src/ctapipe/irf/cuts.py +++ b/src/ctapipe/irf/cuts.py @@ -12,7 +12,7 @@ class EventQualitySelection(QualityQuery): """ - Event pre-selection quality criteria for IRF computation with different defaults. + Event pre-selection quality criteria for IRF and DL3 computation with different defaults. """ quality_criteria = List( @@ -29,10 +29,36 @@ class EventQualitySelection(QualityQuery): help=QualityQuery.quality_criteria.help, ).tag(config=True) - def calculate_selection(self, events): + def calculate_selection(self, events: Table): + """ + Add the selection columns to the events, will only compute quality selection + + Parameters + ---------- + events: Table + The table containing the events on which selection need to be applied + + Returns + ------- + Table + events with selection columns added. + """ return self.calculate_quality_selection(events) - def calculate_quality_selection(self, events): + def calculate_quality_selection(self, events: Table): + """ + Add the selection columns to the events, will only compute quality selection + + Parameters + ---------- + events: Table + The table containing the events on which selection need to be applied + + Returns + ------- + Table + events with selection columns added. + """ events["selected_quality"] = self.get_table_mask(events) events["selected"] = events["selected_quality"] return events @@ -40,7 +66,7 @@ def calculate_quality_selection(self, events): class EventSelection(EventQualitySelection): """ - Event selection + Event selection quality and gammaness criteria for IRF and DL3 """ cuts_file = Path( @@ -50,13 +76,28 @@ class EventSelection(EventQualitySelection): help="Path to the cuts file to apply to the observation.", ).tag(config=True) - def __init__(self, **kwargs): - super().__init__(**kwargs) + def __init__(self, config=None, parent=None, **kwargs): + super().__init__(config=config, parent=parent, **kwargs) self.cuts = OptimizationResult.read(self.cuts_file) - def calculate_selection(self, events): + def calculate_selection(self, events: Table, apply_spatial_selection: bool = False): + """ + Add the selection columns to the events + + Parameters + ---------- + events: Table + The table containing the events on which selection need to be applied + apply_spatial_selection: bool + True if the theta cuts should be applied + + Returns + ------- + Table + events with selection columns added. + """ events = self.calculate_quality_selection(events) - events = self.calculate_gamma_selection(events) + events = self.calculate_gamma_selection(events, apply_spatial_selection) events["selected"] = events["selected_quality"] & events["selected_gamma"] return events @@ -64,14 +105,12 @@ def calculate_gamma_selection( self, events: Table, apply_spatial_selection: bool = False ) -> Table: """ - Add the selection columns to the events + Add the selection columns to the events, will compute only gamma criteria Parameters ---------- events: Table The table containing the events on which selection need to be applied - cuts: OptimizationResult - The cuts that need to be applied on the events apply_spatial_selection: bool True if the theta cuts should be applied From 1af409405edc49ab8f4e8b4a2bc142c5475d696a Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 20 Feb 2025 17:31:34 +0100 Subject: [PATCH 08/58] Update event processor to support gamaness event selection --- src/ctapipe/irf/optimize/test/test_optimize.py | 6 +++--- src/ctapipe/irf/preprocessing.py | 16 +++++++++++----- src/ctapipe/tools/compute_irf.py | 12 ++++++------ src/ctapipe/tools/optimize_event_selection.py | 2 +- 4 files changed, 21 insertions(+), 15 deletions(-) diff --git a/src/ctapipe/irf/optimize/test/test_optimize.py b/src/ctapipe/irf/optimize/test/test_optimize.py index bdd99ecf3f9..5f6c5ea3214 100644 --- a/src/ctapipe/irf/optimize/test/test_optimize.py +++ b/src/ctapipe/irf/optimize/test/test_optimize.py @@ -15,13 +15,13 @@ def test_optimization_result(tmp_path, irf_event_loader_test_config): ) result_path = tmp_path / "result.h5" - epp = EventPreprocessor(irf_event_loader_test_config) + epp = EventPreprocessor(config=irf_event_loader_test_config) gh_cuts = QTable( data=[[0.2, 0.8, 1.5] * u.TeV, [0.8, 1.5, 10] * u.TeV, [0.82, 0.91, 0.88]], names=["low", "high", "cut"], ) result = OptimizationResult( - quality_query=epp.quality_query, + quality_query=epp.event_selection, gh_cuts=gh_cuts, clf_prefix="ExtraTreesClassifier", valid_energy_min=0.2 * u.TeV, @@ -105,7 +105,7 @@ def test_cut_optimizer( optimizer = Optimizer() result = optimizer( events={"signal": gamma_events, "background": proton_events}, - quality_query=gamma_loader.epp.quality_query, # identical qualityquery for all particle types + quality_query=gamma_loader.epp.event_selection, # identical qualityquery for all particle types clf_prefix="ExtraTreesClassifier", ) assert isinstance(result, OptimizationResult) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index b27ddf3bbc3..1d45b2e42c7 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -22,7 +22,7 @@ from ..core import Component from ..core.traits import Unicode from ..io import TableLoader -from .cuts import EventQualitySelection +from .cuts import EventQualitySelection, EventSelection from .spectra import SPECTRA, Spectra __all__ = ["EventLoader", "EventPreprocessor"] @@ -31,7 +31,7 @@ class EventPreprocessor(Component): """Defines pre-selection cuts and the necessary renaming of columns.""" - classes = [EventQualitySelection] + classes = [EventQualitySelection, EventSelection] energy_reconstructor = Unicode( default_value="RandomForestRegressor", @@ -48,9 +48,14 @@ class EventPreprocessor(Component): help="Prefix of the classifier `_prediction` column", ).tag(config=True) - def __init__(self, config=None, parent=None, **kwargs): + def __init__( + self, quality_selection_only: bool = True, config=None, parent=None, **kwargs + ): super().__init__(config=config, parent=parent, **kwargs) - self.quality_query = EventQualitySelection(parent=self) + if quality_selection_only: + self.event_selection = EventQualitySelection(parent=self) + else: + self.event_selection = EventSelection(parent=self) def normalise_column_names(self, events: Table) -> QTable: if events["subarray_pointing_lat"].std() > 1e-3: @@ -200,7 +205,8 @@ def load_preselected_events(self, chunk_size: int) -> tuple[QTable, int, dict]: header = self.epp.make_empty_table() bits = [header] for _, _, events in load.read_subarray_events_chunked(chunk_size, **opts): - selected = events[self.epp.quality_query.get_table_mask(events)] + events = self.epp.event_selection.calculate_selection(events) + selected = events[events["selected"]] selected = self.epp.normalise_column_names(selected) selected = self.make_derived_columns(selected) bits.append(selected) diff --git a/src/ctapipe/tools/compute_irf.py b/src/ctapipe/tools/compute_irf.py index 4f505aad506..7028ff2242d 100644 --- a/src/ctapipe/tools/compute_irf.py +++ b/src/ctapipe/tools/compute_irf.py @@ -471,7 +471,7 @@ def start(self): ) if ( - loader.epp.quality_query.quality_criteria + loader.epp.event_selection.quality_criteria != self.opt_result.quality_query.quality_criteria ): self.log.warning( @@ -479,22 +479,22 @@ def start(self): "calculating g/h / theta cuts. Provided quality criteria:\n%s. " "\nUsing the same quality criteria as g/h / theta cuts:\n%s. " % ( - loader.epp.quality_query.to_table(functions=True)[ + loader.epp.event_selection.to_table(functions=True)[ "criteria", "func" ], - self.opt_result.quality_query.to_table(functions=True)[ + self.opt_result.event_selection.to_table(functions=True)[ "criteria", "func" ], ) ) - loader.epp.quality_query = EventQualitySelection( + loader.epp.event_selection = EventQualitySelection( parent=loader, - quality_criteria=self.opt_result.quality_query.quality_criteria, + quality_criteria=self.opt_result.event_selection.quality_criteria, ) self.log.debug( "%s Quality criteria: %s" - % (particle_type, loader.epp.quality_query.quality_criteria) + % (particle_type, loader.epp.event_selection.quality_criteria) ) events = loader.load_preselected_events(self.chunk_size) count = len(events) diff --git a/src/ctapipe/tools/optimize_event_selection.py b/src/ctapipe/tools/optimize_event_selection.py index 6865a5db402..445c9cf3ac9 100644 --- a/src/ctapipe/tools/optimize_event_selection.py +++ b/src/ctapipe/tools/optimize_event_selection.py @@ -206,7 +206,7 @@ def start(self): self.result = self.optimizer( events={"signal": self.signal_events, "background": self.background_events}, # identical quality_query for all particle types - quality_query=self.event_loaders["gammas"].epp.quality_query, + quality_query=self.event_loaders["gammas"].epp.event_selection, clf_prefix=self.event_loaders["gammas"].epp.gammaness_classifier, ) From 62df52454b3cbbe84ec73843ee365ab64fa92c67 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 20 Feb 2025 17:40:32 +0100 Subject: [PATCH 09/58] Update event loader for gamma event selection --- src/ctapipe/irf/optimize/test/test_optimize.py | 4 ++-- src/ctapipe/irf/preprocessing.py | 12 ++++++++++-- src/ctapipe/irf/tests/test_benchmarks.py | 4 ++-- src/ctapipe/irf/tests/test_preprocessing.py | 2 +- 4 files changed, 15 insertions(+), 7 deletions(-) diff --git a/src/ctapipe/irf/optimize/test/test_optimize.py b/src/ctapipe/irf/optimize/test/test_optimize.py index 5f6c5ea3214..0cb032a12bf 100644 --- a/src/ctapipe/irf/optimize/test/test_optimize.py +++ b/src/ctapipe/irf/optimize/test/test_optimize.py @@ -90,15 +90,15 @@ def test_cut_optimizer( from ctapipe.irf import EventLoader, OptimizationResult, Spectra gamma_loader = EventLoader( - config=irf_event_loader_test_config, file=gamma_diffuse_full_reco_file, target_spectrum=Spectra.CRAB_HEGRA, + config=irf_event_loader_test_config, ) gamma_events = gamma_loader.load_preselected_events(chunk_size=10000) proton_loader = EventLoader( - config=irf_event_loader_test_config, file=proton_full_reco_file, target_spectrum=Spectra.IRFDOC_PROTON_SPECTRUM, + config=irf_event_loader_test_config, ) proton_events = proton_loader.load_preselected_events(chunk_size=10000) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 1d45b2e42c7..00a88e4fe3d 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -192,10 +192,18 @@ class EventLoader(Component): classes = [EventPreprocessor] - def __init__(self, file: Path, target_spectrum: Spectra, **kwargs): + def __init__( + self, + file: Path, + target_spectrum: Spectra, + quality_selection_only: bool = True, + **kwargs, + ): super().__init__(**kwargs) - self.epp = EventPreprocessor(parent=self) + self.epp = EventPreprocessor( + quality_selection_only=quality_selection_only, parent=self + ) self.target_spectrum = SPECTRA[target_spectrum] self.file = file diff --git a/src/ctapipe/irf/tests/test_benchmarks.py b/src/ctapipe/irf/tests/test_benchmarks.py index 076fe0d42fc..4dbd2f6a37f 100644 --- a/src/ctapipe/irf/tests/test_benchmarks.py +++ b/src/ctapipe/irf/tests/test_benchmarks.py @@ -87,15 +87,15 @@ def test_make_2d_sensitivity( from ctapipe.irf.tests.test_irfs import _check_boundaries_in_hdu gamma_loader = EventLoader( - config=irf_event_loader_test_config, file=gamma_diffuse_full_reco_file, target_spectrum=Spectra.CRAB_HEGRA, + config=irf_event_loader_test_config, ) gamma_events = gamma_loader.load_preselected_events(chunk_size=10000) proton_loader = EventLoader( - config=irf_event_loader_test_config, file=proton_full_reco_file, target_spectrum=Spectra.IRFDOC_PROTON_SPECTRUM, + config=irf_event_loader_test_config, ) proton_events = proton_loader.load_preselected_events(chunk_size=10000) diff --git a/src/ctapipe/irf/tests/test_preprocessing.py b/src/ctapipe/irf/tests/test_preprocessing.py index 0244933a6e6..8daa90ebe23 100644 --- a/src/ctapipe/irf/tests/test_preprocessing.py +++ b/src/ctapipe/irf/tests/test_preprocessing.py @@ -67,9 +67,9 @@ def test_event_loader(gamma_diffuse_full_reco_file, irf_event_loader_test_config from ctapipe.irf import EventLoader, Spectra loader = EventLoader( - config=irf_event_loader_test_config, file=gamma_diffuse_full_reco_file, target_spectrum=Spectra.CRAB_HEGRA, + config=irf_event_loader_test_config, ) events = loader.load_preselected_events(chunk_size=10000) count = len(events) From 50c1af4799cec8e1988226662da2a9add79365bf Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 20 Feb 2025 17:44:22 +0100 Subject: [PATCH 10/58] Update path criteria for input --- src/ctapipe/irf/cuts.py | 1 + src/ctapipe/tools/create_dl3.py | 4 ++++ 2 files changed, 5 insertions(+) diff --git a/src/ctapipe/irf/cuts.py b/src/ctapipe/irf/cuts.py index 59a1cf136a9..2e1458d5174 100644 --- a/src/ctapipe/irf/cuts.py +++ b/src/ctapipe/irf/cuts.py @@ -73,6 +73,7 @@ class EventSelection(EventQualitySelection): default_value=None, allow_none=False, directory_ok=False, + exists=True, help="Path to the cuts file to apply to the observation.", ).tag(config=True) diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index c0ea3c20b78..b4c602ab0e9 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -24,6 +24,7 @@ class DL3Tool(Tool): default_value=None, allow_none=False, directory_ok=False, + exists=True, help="DL2 input filename and path.", ).tag(config=True) @@ -31,6 +32,7 @@ class DL3Tool(Tool): default_value=None, allow_none=False, directory_ok=False, + exists=True, help="Output file", ).tag(config=True) @@ -38,6 +40,7 @@ class DL3Tool(Tool): default_value=None, allow_none=False, directory_ok=False, + exists=True, help="Path to the cuts file to apply to the observation.", ).tag(config=True) @@ -45,6 +48,7 @@ class DL3Tool(Tool): default_value=None, allow_none=False, directory_ok=False, + exists=True, help="Path to the IRFs describing the observation", ).tag(config=True) From aa64f93f4e9123fa49a5755853d9c354805e381d Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Fri, 21 Feb 2025 10:59:36 +0100 Subject: [PATCH 11/58] Further update to prepare for DL3 --- src/ctapipe/irf/cuts.py | 1 - src/ctapipe/irf/preprocessing.py | 45 +++++++++++++++++++++++++++----- 2 files changed, 38 insertions(+), 8 deletions(-) diff --git a/src/ctapipe/irf/cuts.py b/src/ctapipe/irf/cuts.py index 2e1458d5174..8e0425995da 100644 --- a/src/ctapipe/irf/cuts.py +++ b/src/ctapipe/irf/cuts.py @@ -70,7 +70,6 @@ class EventSelection(EventQualitySelection): """ cuts_file = Path( - default_value=None, allow_none=False, directory_ok=False, exists=True, diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 00a88e4fe3d..a54c3612d43 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -98,9 +98,31 @@ def normalise_column_names(self, events: Table) -> QTable: f"Required column {c} is missing." ) - events = QTable(events[keep_columns], copy=COPY_IF_NEEDED) - events.rename_columns(rename_from, rename_to) - return events + renamed_events = QTable(events, copy=COPY_IF_NEEDED) + renamed_events.rename_columns(rename_from, rename_to) + return renamed_events + + def keep_nescessary_columns_only(self, events: Table) -> QTable: + keep_columns = [ + "obs_id", + "event_id", + "true_energy", + "true_az", + "true_alt", + "reco_energy", + "reco_az", + "reco_alt", + "reco_fov_lat", + "reco_fov_lon", + "pointing_az", + "pointing_alt", + "theta", + "true_source_fov_offset", + "reco_source_fov_offset", + "gh_score", + "weight", + ] + return QTable(events[keep_columns], copy=COPY_IF_NEEDED) def make_empty_table(self) -> QTable: """ @@ -195,7 +217,7 @@ class EventLoader(Component): def __init__( self, file: Path, - target_spectrum: Spectra, + target_spectrum: Spectra = None, quality_selection_only: bool = True, **kwargs, ): @@ -204,7 +226,10 @@ def __init__( self.epp = EventPreprocessor( quality_selection_only=quality_selection_only, parent=self ) - self.target_spectrum = SPECTRA[target_spectrum] + if target_spectrum is not None: + self.target_spectrum = SPECTRA[target_spectrum] + else: + self.target_spectrum = None self.file = file def load_preselected_events(self, chunk_size: int) -> tuple[QTable, int, dict]: @@ -213,10 +238,11 @@ def load_preselected_events(self, chunk_size: int) -> tuple[QTable, int, dict]: header = self.epp.make_empty_table() bits = [header] for _, _, events in load.read_subarray_events_chunked(chunk_size, **opts): + events = self.epp.normalise_column_names(events) + events = self.make_derived_columns(events) events = self.epp.event_selection.calculate_selection(events) selected = events[events["selected"]] - selected = self.epp.normalise_column_names(selected) - selected = self.make_derived_columns(selected) + selected = self.epp.keep_nescessary_columns_only(selected) bits.append(selected) bits.append(header) # Putting it last ensures the correct metadata is used @@ -296,6 +322,11 @@ def make_event_weights( kind: str, fov_offset_bins: u.Quantity | None = None, ) -> QTable: + if self.target_spectrum is None: + raise Exception( + "No spectrum is defined, need a spectrum for events weighting" + ) + if ( kind == "gammas" and self.target_spectrum.normalization.unit.is_equivalent( From 207cb7c8945a422dfd5cd0fe43bcd3dd5520c353 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Fri, 21 Feb 2025 11:02:24 +0100 Subject: [PATCH 12/58] Skeleton for creating DL3 --- pyproject.toml | 1 + src/ctapipe/tools/create_dl3.py | 48 ++++++++++++--------------------- 2 files changed, 18 insertions(+), 31 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 42f7ea9d9c7..ee2880dd6bc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -108,6 +108,7 @@ ctapipe-process = "ctapipe.tools.process:main" ctapipe-merge = "ctapipe.tools.merge:main" ctapipe-optimize-event-selection = "ctapipe.tools.optimize_event_selection:main" ctapipe-compute-irf = "ctapipe.tools.compute_irf:main" +ctapipe-create-dl3 = "ctapipe.tools.create_dl3:main" ctapipe-fileinfo = "ctapipe.tools.fileinfo:main" ctapipe-quickstart = "ctapipe.tools.quickstart:main" ctapipe-calculate-pixel-statistics = "ctapipe.tools.calculate_pixel_stats:main" diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index b4c602ab0e9..194c176d4b2 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -1,27 +1,18 @@ from ctapipe.core import Tool, traits -from ctapipe.core.traits import ( - Integer, -) +from ctapipe.core.traits import Integer, classes_with_traits -from ..irf import ( - EventLoader, - OptimizationResult, -) +from ..irf import EventLoader __all__ = ["DL3Tool"] +from ..irf.cuts import EventSelection + class DL3Tool(Tool): - name = "mytool" - description = "do some things and stuff" - aliases = dict( - infile="AdvancedComponent.infile", - outfile="AdvancedComponent.outfile", - iterations="MyTool.iterations", - ) + name = "ctapipe-create-dl3" + description = "Create DL3 file from DL2 observation file" dl2_file = traits.Path( - default_value=None, allow_none=False, directory_ok=False, exists=True, @@ -29,23 +20,12 @@ class DL3Tool(Tool): ).tag(config=True) output_path = traits.Path( - default_value=None, allow_none=False, directory_ok=False, - exists=True, help="Output file", ).tag(config=True) - cuts_file = traits.Path( - default_value=None, - allow_none=False, - directory_ok=False, - exists=True, - help="Path to the cuts file to apply to the observation.", - ).tag(config=True) - irfs_file = traits.Path( - default_value=None, allow_none=False, directory_ok=False, exists=True, @@ -61,19 +41,25 @@ class DL3Tool(Tool): # Which classes are registered for configuration classes = [ EventLoader, - ] + ] + classes_with_traits(EventSelection) + + aliases = { + "cuts": "EventSelection.cuts_file", + "dl2-file": "DL3Tool.dl2_file", + "irfs-file": "DL3Tool.irfs_file", + "output": "DL3Tool.output_path", + "chunk-size": "DL3Tool.chunk_size", + } def setup(self): """ Initialize components from config and load g/h (and theta) cuts. """ - self.log.info("Loading cuts") - self.opt_result = OptimizationResult.read(self.cuts_file) self.log.info("Loading events from DL2") self.event_loader = EventLoader( - parent=self, - file=self.dl2_file, + parent=self, file=self.dl2_file, quality_selection_only=True ) + print(self.event_loader.load_preselected_events(self.chunk_size)) def start(self): pass From 2ea79cc210b103e3700fb2e1a3a43b03124fa944 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Fri, 21 Feb 2025 11:54:39 +0100 Subject: [PATCH 13/58] Basic DL3 structure --- src/ctapipe/tools/create_dl3.py | 35 +++++++++++++++++++++++++-------- 1 file changed, 27 insertions(+), 8 deletions(-) diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index 194c176d4b2..0220f3589eb 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -1,5 +1,7 @@ +from astropy.io import fits + from ctapipe.core import Tool, traits -from ctapipe.core.traits import Integer, classes_with_traits +from ctapipe.core.traits import Bool, Integer, classes_with_traits from ..irf import EventLoader @@ -19,7 +21,7 @@ class DL3Tool(Tool): help="DL2 input filename and path.", ).tag(config=True) - output_path = traits.Path( + output_file = traits.Path( allow_none=False, directory_ok=False, help="Output file", @@ -38,6 +40,11 @@ class DL3Tool(Tool): help="How many subarray events to load at once while selecting.", ).tag(config=True) + overwrite = Bool( + default_value=False, + help="If true, allow to overwrite already existing output file", + ).tag(config=True) + # Which classes are registered for configuration classes = [ EventLoader, @@ -47,25 +54,37 @@ class DL3Tool(Tool): "cuts": "EventSelection.cuts_file", "dl2-file": "DL3Tool.dl2_file", "irfs-file": "DL3Tool.irfs_file", - "output": "DL3Tool.output_path", + "output": "DL3Tool.output_file", "chunk-size": "DL3Tool.chunk_size", + "overwrite": "DL3Tool.overwrite", } def setup(self): """ Initialize components from config and load g/h (and theta) cuts. """ + pass + + def start(self): self.log.info("Loading events from DL2") self.event_loader = EventLoader( - parent=self, file=self.dl2_file, quality_selection_only=True + parent=self, file=self.dl2_file, quality_selection_only=False ) - print(self.event_loader.load_preselected_events(self.chunk_size)) + events = self.event_loader.load_preselected_events(self.chunk_size) - def start(self): - pass + hdu_dl3 = fits.HDUList([fits.PrimaryHDU()]) + hdu_dl3.append(fits.BinTableHDU(data=events, name="EVENTS")) + + self.log.info("Loading IRFs") + hdu_irfs = fits.open(self.irfs_file, checksum=True) + for i in range(1, len(hdu_irfs)): + hdu_dl3.append(hdu_irfs[i]) + + self.log.info("Writing DL3 File") + hdu_dl3.writeto(self.output_file, checksum=True, overwrite=self.overwrite) def finish(self): - self.log.warning("Shutting down.") + pass def main(): From 7811ab139d8f60cb8b175edf1a1c3b88f6d69027 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Mon, 3 Mar 2025 12:17:28 +0100 Subject: [PATCH 14/58] Add option for optional DL3 columns --- src/ctapipe/irf/dl3.py | 2 ++ src/ctapipe/irf/preprocessing.py | 6 +++++- src/ctapipe/tools/create_dl3.py | 16 ++++++++++++++-- 3 files changed, 21 insertions(+), 3 deletions(-) create mode 100644 src/ctapipe/irf/dl3.py diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py new file mode 100644 index 00000000000..0e684bd8c50 --- /dev/null +++ b/src/ctapipe/irf/dl3.py @@ -0,0 +1,2 @@ +def get_hdu_header_events(): + return {"HDUCLASS": "GADF", "HDUCLAS1": "EVENTS"} diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index a54c3612d43..b7bfa018ffb 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -20,7 +20,7 @@ from ..containers import CoordinateFrameType from ..coordinates import NominalFrame from ..core import Component -from ..core.traits import Unicode +from ..core.traits import Bool, Unicode from ..io import TableLoader from .cuts import EventQualitySelection, EventSelection from .spectra import SPECTRA, Spectra @@ -48,6 +48,10 @@ class EventPreprocessor(Component): help="Prefix of the classifier `_prediction` column", ).tag(config=True) + optional_dl3_columns = Bool( + default_value=False, help="If true add optional columns to produce file" + ).tag(config=True) + def __init__( self, quality_selection_only: bool = True, config=None, parent=None, **kwargs ): diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index 0220f3589eb..bebfcaa195b 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -1,3 +1,5 @@ +from datetime import datetime + from astropy.io import fits from ctapipe.core import Tool, traits @@ -8,6 +10,7 @@ __all__ = ["DL3Tool"] from ..irf.cuts import EventSelection +from ..irf.dl3 import get_hdu_header_events class DL3Tool(Tool): @@ -55,6 +58,7 @@ class DL3Tool(Tool): "dl2-file": "DL3Tool.dl2_file", "irfs-file": "DL3Tool.irfs_file", "output": "DL3Tool.output_file", + "optional-column": "EventPreprocessor.optional_dl3_columns", "chunk-size": "DL3Tool.chunk_size", "overwrite": "DL3Tool.overwrite", } @@ -72,8 +76,16 @@ def start(self): ) events = self.event_loader.load_preselected_events(self.chunk_size) - hdu_dl3 = fits.HDUList([fits.PrimaryHDU()]) - hdu_dl3.append(fits.BinTableHDU(data=events, name="EVENTS")) + hdu_dl3 = fits.HDUList( + [ + fits.PrimaryHDU( + header={"CREATED": datetime.now(tz=datetime.UTC).isoformat()} + ) + ] + ) + hdu_dl3.append( + fits.BinTableHDU(data=events, name="EVENTS", header=get_hdu_header_events()) + ) self.log.info("Loading IRFs") hdu_irfs = fits.open(self.irfs_file, checksum=True) From 861d2fa6d523591d91cbd23e58c1acfc8eefa4f9 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Mon, 3 Mar 2025 12:21:15 +0100 Subject: [PATCH 15/58] Change option for flag --- src/ctapipe/tools/create_dl3.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index bebfcaa195b..55336258939 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -3,7 +3,7 @@ from astropy.io import fits from ctapipe.core import Tool, traits -from ctapipe.core.traits import Bool, Integer, classes_with_traits +from ctapipe.core.traits import Bool, Integer, classes_with_traits, flag from ..irf import EventLoader @@ -63,6 +63,15 @@ class DL3Tool(Tool): "overwrite": "DL3Tool.overwrite", } + flags = { + **flag( + "optional-column", + "EventPreprocessor.optional_dl3_columns", + "Add optional columns for events in the DL3 file", + "Do not add optional column for events in the DL3 file", + ), + } + def setup(self): """ Initialize components from config and load g/h (and theta) cuts. From bf0cd867bc64cf9624d42d282d48baa75d1f5f2a Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Mon, 3 Mar 2025 12:22:43 +0100 Subject: [PATCH 16/58] Fix the option --- src/ctapipe/tools/create_dl3.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index 55336258939..4ddaff2be68 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -5,7 +5,7 @@ from ctapipe.core import Tool, traits from ctapipe.core.traits import Bool, Integer, classes_with_traits, flag -from ..irf import EventLoader +from ..irf import EventLoader, EventPreprocessor __all__ = ["DL3Tool"] @@ -49,16 +49,19 @@ class DL3Tool(Tool): ).tag(config=True) # Which classes are registered for configuration - classes = [ - EventLoader, - ] + classes_with_traits(EventSelection) + classes = ( + [ + EventLoader, + ] + + classes_with_traits(EventSelection) + + classes_with_traits(EventPreprocessor) + ) aliases = { "cuts": "EventSelection.cuts_file", "dl2-file": "DL3Tool.dl2_file", "irfs-file": "DL3Tool.irfs_file", "output": "DL3Tool.output_file", - "optional-column": "EventPreprocessor.optional_dl3_columns", "chunk-size": "DL3Tool.chunk_size", "overwrite": "DL3Tool.overwrite", } From 4b3765c0dd907ae4442ff55c049ef9e72674a20a Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Mon, 3 Mar 2025 19:53:49 +0100 Subject: [PATCH 17/58] Add DL3 and IRF mode for preprocessor --- src/ctapipe/irf/preprocessing.py | 5 +++++ src/ctapipe/tools/compute_irf.py | 5 +++++ src/ctapipe/tools/create_dl3.py | 4 +++- 3 files changed, 13 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index b7bfa018ffb..77b1c954f69 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -48,6 +48,11 @@ class EventPreprocessor(Component): help="Prefix of the classifier `_prediction` column", ).tag(config=True) + irf_pre_processing = Bool( + default_value=True, + help="If true the pre processing assume the purpose is for IRF production, if false, for DL3 production", + ).tag(config=False) + optional_dl3_columns = Bool( default_value=False, help="If true add optional columns to produce file" ).tag(config=True) diff --git a/src/ctapipe/tools/compute_irf.py b/src/ctapipe/tools/compute_irf.py index 7028ff2242d..15e2e713a02 100644 --- a/src/ctapipe/tools/compute_irf.py +++ b/src/ctapipe/tools/compute_irf.py @@ -21,6 +21,7 @@ from ..core.traits import AstroQuantity, Bool, Integer, classes_with_traits, flag from ..irf import ( EventLoader, + EventPreprocessor, OptimizationResult, Spectra, check_bins_in_range, @@ -231,6 +232,10 @@ def setup(self): """ Initialize components from config and load g/h (and theta) cuts. """ + + # Force the preprocessing for IRF + EventPreprocessor.irf_pre_processing = True + self.opt_result = OptimizationResult.read(self.cuts_file) if ( self.spatial_selection_applied diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index 4ddaff2be68..49dd8f5448f 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -79,7 +79,9 @@ def setup(self): """ Initialize components from config and load g/h (and theta) cuts. """ - pass + + # Force the preprocessing for DL3 + EventPreprocessor.irf_pre_processing = False def start(self): self.log.info("Loading events from DL2") From 9ff2f21e9eff85078104caff9e088d4e94403897 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Mon, 3 Mar 2025 20:08:00 +0100 Subject: [PATCH 18/58] Adapt preprocessor columns request --- src/ctapipe/irf/preprocessing.py | 40 +++++++++++++++++++++++++++++--- 1 file changed, 37 insertions(+), 3 deletions(-) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 77b1c954f69..a268ffb87a8 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -79,10 +79,22 @@ def normalise_column_names(self, events: Table) -> QTable: keep_columns = [ "obs_id", "event_id", - "true_energy", - "true_az", - "true_alt", ] + if self.irf_pre_processing: + keep_columns += [ + "true_energy", + "true_az", + "true_alt", + ] + else: + keep_columns += [ + "time", + ] + if self.optional_dl3_columns: + keep_columns += [ + "tels_with_trigger", + ] + rename_from = [ f"{self.energy_reconstructor}_energy", f"{self.geometry_reconstructor}_az", @@ -99,6 +111,28 @@ def normalise_column_names(self, events: Table) -> QTable: "pointing_alt", "pointing_az", ] + if self.optional_dl3_columns: + rename_from = [ + f"{self.energy_reconstructor}_energy_uncert", + f"{self.geometry_reconstructor}_ang_distance_uncert", + f"{self.geometry_reconstructor}_core_x", + f"{self.geometry_reconstructor}_core_y", + f"{self.geometry_reconstructor}_core_uncert_x", + f"{self.geometry_reconstructor}_core_uncert_y", + f"{self.geometry_reconstructor}_h_max", + f"{self.geometry_reconstructor}_h_max_uncert", + ] + rename_to = [ + "reco_energy_uncert", + "reco_dir_uncert", + "reco_core_x", + "reco_core_y", + "reco_core_uncert_x", + "reco_core_uncert_y", + "reco_h_max", + "reco_h_max_uncert", + ] + keep_columns.extend(rename_from) for c in keep_columns: if c not in events.colnames: From 912d871d0b619729436606f117ff4b5fe0fa38d2 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Mon, 3 Mar 2025 20:17:54 +0100 Subject: [PATCH 19/58] Add option to skip missing optional column --- src/ctapipe/irf/preprocessing.py | 21 +++++++++++++++++++-- src/ctapipe/tools/create_dl3.py | 6 ++++++ 2 files changed, 25 insertions(+), 2 deletions(-) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index a268ffb87a8..a12ad1c314a 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -57,6 +57,11 @@ class EventPreprocessor(Component): default_value=False, help="If true add optional columns to produce file" ).tag(config=True) + raise_error_for_optional = Bool( + default_value=True, + help="If true will raise error in the case optional column are missing", + ).tag(config=True) + def __init__( self, quality_selection_only: bool = True, config=None, parent=None, **kwargs ): @@ -112,7 +117,7 @@ def normalise_column_names(self, events: Table) -> QTable: "pointing_az", ] if self.optional_dl3_columns: - rename_from = [ + rename_from_optional = [ f"{self.energy_reconstructor}_energy_uncert", f"{self.geometry_reconstructor}_ang_distance_uncert", f"{self.geometry_reconstructor}_core_x", @@ -122,7 +127,7 @@ def normalise_column_names(self, events: Table) -> QTable: f"{self.geometry_reconstructor}_h_max", f"{self.geometry_reconstructor}_h_max_uncert", ] - rename_to = [ + rename_to_optional = [ "reco_energy_uncert", "reco_dir_uncert", "reco_core_x", @@ -132,6 +137,18 @@ def normalise_column_names(self, events: Table) -> QTable: "reco_h_max", "reco_h_max_uncert", ] + if not self.raise_error_for_optional: + for i, c in enumerate(rename_from_optional): + if c not in events.colnames: + self.log.warning( + f"Optional column {c} is missing from the DL2 file." + ) + else: + rename_from.append(rename_from_optional[i]) + rename_to.append(rename_to_optional[i]) + else: + rename_from += rename_from_optional + rename_to += rename_to_optional keep_columns.extend(rename_from) for c in keep_columns: diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index 49dd8f5448f..a44f2e69ddb 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -73,6 +73,12 @@ class DL3Tool(Tool): "Add optional columns for events in the DL3 file", "Do not add optional column for events in the DL3 file", ), + **flag( + "raise-error-for-optional", + "EventPreprocessor.raise_error_for_optional", + "Raise an error if an optional column is missing", + "Only display a warning if an optional column is missing, it will lead to optional columns missing in the DL3 file", + ), } def setup(self): From 8336606f4d62551379c8d72ebc8dfaa63bc11e5a Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Mon, 3 Mar 2025 20:19:30 +0100 Subject: [PATCH 20/58] Allow for variation of the pointing for DL3 --- src/ctapipe/irf/preprocessing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index a12ad1c314a..1ab77e335cf 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -72,7 +72,7 @@ def __init__( self.event_selection = EventSelection(parent=self) def normalise_column_names(self, events: Table) -> QTable: - if events["subarray_pointing_lat"].std() > 1e-3: + if self.irf_pre_processing and events["subarray_pointing_lat"].std() > 1e-3: raise NotImplementedError( "No support for making irfs from varying pointings yet" ) From 636fbae8ad3a09c3baf9030be66c504668147502 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 6 Mar 2025 13:00:07 +0100 Subject: [PATCH 21/58] Move the derived columns function --- src/ctapipe/irf/preprocessing.py | 66 ++++++++++++++++---------------- 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 1ab77e335cf..b91af478698 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -184,6 +184,38 @@ def keep_nescessary_columns_only(self, events: Table) -> QTable: ] return QTable(events[keep_columns], copy=COPY_IF_NEEDED) + def make_derived_columns(self, events: QTable) -> QTable: + events["weight"] = ( + 1.0 * u.dimensionless_unscaled + ) # defer calculation of proper weights to later + events["gh_score"].unit = u.dimensionless_unscaled + 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" + ) + + altaz = AltAz() + pointing = SkyCoord( + alt=events["pointing_alt"], az=events["pointing_az"], frame=altaz + ) + reco = SkyCoord( + alt=events["reco_alt"], + az=events["reco_az"], + frame=altaz, + ) + nominal = NominalFrame(origin=pointing) + reco_nominal = reco.transform_to(nominal) + events["reco_fov_lon"] = u.Quantity(-reco_nominal.fov_lon) # minus for GADF + events["reco_fov_lat"] = u.Quantity(reco_nominal.fov_lat) + return events + def make_empty_table(self) -> QTable: """ This function defines the columns later functions expect to be present @@ -299,7 +331,7 @@ def load_preselected_events(self, chunk_size: int) -> tuple[QTable, int, dict]: bits = [header] for _, _, events in load.read_subarray_events_chunked(chunk_size, **opts): events = self.epp.normalise_column_names(events) - events = self.make_derived_columns(events) + events = self.epp.make_derived_columns(events) events = self.epp.event_selection.calculate_selection(events) selected = events[events["selected"]] selected = self.epp.keep_nescessary_columns_only(selected) @@ -343,38 +375,6 @@ def get_simulation_information( } return meta - def make_derived_columns(self, events: QTable) -> QTable: - events["weight"] = ( - 1.0 * u.dimensionless_unscaled - ) # defer calculation of proper weights to later - events["gh_score"].unit = u.dimensionless_unscaled - 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" - ) - - altaz = AltAz() - pointing = SkyCoord( - alt=events["pointing_alt"], az=events["pointing_az"], frame=altaz - ) - reco = SkyCoord( - alt=events["reco_alt"], - az=events["reco_az"], - frame=altaz, - ) - nominal = NominalFrame(origin=pointing) - reco_nominal = reco.transform_to(nominal) - events["reco_fov_lon"] = u.Quantity(-reco_nominal.fov_lon) # minus for GADF - events["reco_fov_lat"] = u.Quantity(reco_nominal.fov_lat) - return events - def make_event_weights( self, events: QTable, From 15e089f465f4add3a29ac8b69830ee9606806a22 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 6 Mar 2025 13:46:32 +0100 Subject: [PATCH 22/58] Update make_derived_columns function for producing columns needed for DL3 --- src/ctapipe/irf/preprocessing.py | 84 ++++++++++++++++++++++---------- 1 file changed, 58 insertions(+), 26 deletions(-) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index b91af478698..f3f41c61af7 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -4,8 +4,9 @@ import astropy.units as u import numpy as np -from astropy.coordinates import AltAz, SkyCoord +from astropy.coordinates import ICRS, AltAz, EarthLocation, Galactic, SkyCoord from astropy.table import Column, QTable, Table, vstack +from astropy.time import Time from pyirf.simulations import SimulatedEventsInfo from pyirf.spectral import ( DIFFUSE_FLUX_UNIT, @@ -184,36 +185,65 @@ def keep_nescessary_columns_only(self, events: Table) -> QTable: ] return QTable(events[keep_columns], copy=COPY_IF_NEEDED) - def make_derived_columns(self, events: QTable) -> QTable: - events["weight"] = ( - 1.0 * u.dimensionless_unscaled - ) # defer calculation of proper weights to later + def make_derived_columns( + self, events: QTable, location_subarray: EarthLocation + ) -> QTable: + """ + This function compute all the derived columns necessary either to irf production or DL3 file + """ + events["gh_score"].unit = u.dimensionless_unscaled - 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" - ) - altaz = AltAz() - pointing = SkyCoord( - alt=events["pointing_alt"], az=events["pointing_az"], frame=altaz - ) + if self.irf_pre_processing: + frame_subarray = AltAz() + else: + frame_subarray = AltAz( + location=location_subarray, obstime=Time(events["time"]) + ) reco = SkyCoord( alt=events["reco_alt"], az=events["reco_az"], - frame=altaz, + frame=frame_subarray, ) - nominal = NominalFrame(origin=pointing) - reco_nominal = reco.transform_to(nominal) - events["reco_fov_lon"] = u.Quantity(-reco_nominal.fov_lon) # minus for GADF - events["reco_fov_lat"] = u.Quantity(reco_nominal.fov_lat) + + if self.irf_pre_processing or self.optional_dl3_columns: + pointing = SkyCoord( + alt=events["pointing_alt"], + az=events["pointing_az"], + frame=frame_subarray, + ) + nominal = NominalFrame(origin=pointing) + reco_nominal = reco.transform_to(nominal) + events["reco_fov_lon"] = u.Quantity(-reco_nominal.fov_lon) # minus for GADF + events["reco_fov_lat"] = u.Quantity(reco_nominal.fov_lat) + + if self.irf_pre_processing: + events["weight"] = ( + 1.0 * u.dimensionless_unscaled + ) # defer calculation of proper weights to later + 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" + ) + + else: + reco_icrs = reco.transform_to(ICRS()) + events["reco_ra"] = reco_icrs.ra + events["reco_dec"] = reco_icrs.dec + if self.optional_dl3_columns: + reco_gal = reco_icrs.transform_to(Galactic()) + events["reco_glon"] = reco_gal.l + events["reco_glat"] = reco_gal.b + + events["multiplicity"] = np.sum(events["tels_with_trigger"], axis=1) + return events def make_empty_table(self) -> QTable: @@ -331,7 +361,9 @@ def load_preselected_events(self, chunk_size: int) -> tuple[QTable, int, dict]: bits = [header] for _, _, events in load.read_subarray_events_chunked(chunk_size, **opts): events = self.epp.normalise_column_names(events) - events = self.epp.make_derived_columns(events) + events = self.epp.make_derived_columns( + events, load.subarray.reference_location + ) events = self.epp.event_selection.calculate_selection(events) selected = events[events["selected"]] selected = self.epp.keep_nescessary_columns_only(selected) From df6a9c63f2c23ba73a1dc7ac5662b6e7f815e286 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 6 Mar 2025 15:00:01 +0100 Subject: [PATCH 23/58] Update columns selection for DL3 production --- src/ctapipe/irf/preprocessing.py | 209 +++++++++++++++++++------------ 1 file changed, 126 insertions(+), 83 deletions(-) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index f3f41c61af7..6020959f8cd 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -14,7 +14,11 @@ PowerLaw, calculate_event_weights, ) -from pyirf.utils import calculate_source_fov_offset, calculate_theta +from pyirf.utils import ( + calculate_source_fov_offset, + calculate_source_fov_position_angle, + calculate_theta, +) from tables import NoSuchNodeError from ..compat import COPY_IF_NEEDED @@ -72,16 +76,12 @@ def __init__( else: self.event_selection = EventSelection(parent=self) - def normalise_column_names(self, events: Table) -> QTable: - if self.irf_pre_processing and events["subarray_pointing_lat"].std() > 1e-3: - raise NotImplementedError( - "No support for making irfs from varying pointings yet" - ) - if any(events["subarray_pointing_frame"] != CoordinateFrameType.ALTAZ.value): - raise NotImplementedError( - "At the moment only pointing in altaz is supported." - ) - + def get_columns_keep_rename_scheme( + self, events: Table, already_derived: bool = False + ): + """ + Function to get the columns to keep, and the scheme for renaming columns + """ keep_columns = [ "obs_id", "event_id", @@ -101,56 +101,108 @@ def normalise_column_names(self, events: Table) -> QTable: "tels_with_trigger", ] - rename_from = [ - f"{self.energy_reconstructor}_energy", - f"{self.geometry_reconstructor}_az", - f"{self.geometry_reconstructor}_alt", - f"{self.gammaness_classifier}_prediction", - "subarray_pointing_lat", - "subarray_pointing_lon", - ] - rename_to = [ - "reco_energy", - "reco_az", - "reco_alt", - "gh_score", - "pointing_alt", - "pointing_az", - ] - if self.optional_dl3_columns: - rename_from_optional = [ - f"{self.energy_reconstructor}_energy_uncert", - f"{self.geometry_reconstructor}_ang_distance_uncert", - f"{self.geometry_reconstructor}_core_x", - f"{self.geometry_reconstructor}_core_y", - f"{self.geometry_reconstructor}_core_uncert_x", - f"{self.geometry_reconstructor}_core_uncert_y", - f"{self.geometry_reconstructor}_h_max", - f"{self.geometry_reconstructor}_h_max_uncert", + if already_derived: + rename_from = [] + rename_to = [] + keep_columns += ["reco_energy", "reco_az", "reco_alt", "gh_score"] + + if self.irf_pre_processing or self.optional_dl3_columns: + keep_columns += ["reco_fov_lat", "reco_fov_lon"] + + if self.irf_pre_processing: + keep_columns += [ + "pointing_az", + "pointing_alt", + "theta", + "true_source_fov_offset", + "reco_source_fov_offset", + "weight", + ] + else: + keep_columns += ["reco_ra", "reco_dec"] + if self.optional_dl3_columns: + keep_columns += [ + "multiplicity", + "reco_glon", + "reco_glat", + "reco_source_fov_offset", + "reco_source_fov_position_angle", + "reco_energy_uncert", + "reco_dir_uncert", + "reco_core_x", + "reco_core_y", + "reco_core_uncert_x", + "reco_core_uncert_y", + "reco_h_max", + "reco_h_max_uncert", + ] + + else: + rename_from = [ + f"{self.energy_reconstructor}_energy", + f"{self.geometry_reconstructor}_az", + f"{self.geometry_reconstructor}_alt", + f"{self.gammaness_classifier}_prediction", + "subarray_pointing_lat", + "subarray_pointing_lon", ] - rename_to_optional = [ - "reco_energy_uncert", - "reco_dir_uncert", - "reco_core_x", - "reco_core_y", - "reco_core_uncert_x", - "reco_core_uncert_y", - "reco_h_max", - "reco_h_max_uncert", + rename_to = [ + "reco_energy", + "reco_az", + "reco_alt", + "gh_score", + "pointing_alt", + "pointing_az", ] - if not self.raise_error_for_optional: - for i, c in enumerate(rename_from_optional): - if c not in events.colnames: - self.log.warning( - f"Optional column {c} is missing from the DL2 file." - ) - else: - rename_from.append(rename_from_optional[i]) - rename_to.append(rename_to_optional[i]) - else: - rename_from += rename_from_optional - rename_to += rename_to_optional + if self.optional_dl3_columns: + rename_from_optional = [ + f"{self.energy_reconstructor}_energy_uncert", + f"{self.geometry_reconstructor}_ang_distance_uncert", + f"{self.geometry_reconstructor}_core_x", + f"{self.geometry_reconstructor}_core_y", + f"{self.geometry_reconstructor}_core_uncert_x", + f"{self.geometry_reconstructor}_core_uncert_y", + f"{self.geometry_reconstructor}_h_max", + f"{self.geometry_reconstructor}_h_max_uncert", + ] + rename_to_optional = [ + "reco_energy_uncert", + "reco_dir_uncert", + "reco_core_x", + "reco_core_y", + "reco_core_uncert_x", + "reco_core_uncert_y", + "reco_h_max", + "reco_h_max_uncert", + ] + if not self.raise_error_for_optional: + for i, c in enumerate(rename_from_optional): + if c not in events.colnames: + self.log.warning( + f"Optional column {c} is missing from the DL2 file." + ) + else: + rename_from.append(rename_from_optional[i]) + rename_to.append(rename_to_optional[i]) + else: + rename_from += rename_from_optional + rename_to += rename_to_optional + + return keep_columns, rename_from, rename_to + def normalise_column_names(self, events: Table) -> QTable: + if self.irf_pre_processing and events["subarray_pointing_lat"].std() > 1e-3: + raise NotImplementedError( + "No support for making irfs from varying pointings yet" + ) + if any(events["subarray_pointing_frame"] != CoordinateFrameType.ALTAZ.value): + raise NotImplementedError( + "At the moment only pointing in altaz is supported." + ) + + keep_columns, rename_from, rename_to = self.get_columns_keep_rename_scheme( + events + ) keep_columns.extend(rename_from) for c in keep_columns: if c not in events.colnames: @@ -163,26 +215,14 @@ def normalise_column_names(self, events: Table) -> QTable: renamed_events.rename_columns(rename_from, rename_to) return renamed_events - def keep_nescessary_columns_only(self, events: Table) -> QTable: - keep_columns = [ - "obs_id", - "event_id", - "true_energy", - "true_az", - "true_alt", - "reco_energy", - "reco_az", - "reco_alt", - "reco_fov_lat", - "reco_fov_lon", - "pointing_az", - "pointing_alt", - "theta", - "true_source_fov_offset", - "reco_source_fov_offset", - "gh_score", - "weight", - ] + def keep_necessary_columns_only(self, events: Table) -> QTable: + keep_columns, _, _ = self.get_columns_keep_rename_scheme(events) + for c in keep_columns: + if c not in events.colnames: + raise ValueError( + f"Required column {c} is missing from the events table." + ) + return QTable(events[keep_columns], copy=COPY_IF_NEEDED) def make_derived_columns( @@ -216,6 +256,13 @@ def make_derived_columns( reco_nominal = reco.transform_to(nominal) events["reco_fov_lon"] = u.Quantity(-reco_nominal.fov_lon) # minus for GADF events["reco_fov_lat"] = u.Quantity(reco_nominal.fov_lat) + events["reco_source_fov_offset"] = calculate_source_fov_offset( + events, prefix="reco" + ) + if not self.irf_pre_processing: + events[ + "reco_source_fov_position_angle" + ] = calculate_source_fov_position_angle(events, prefix="reco") if self.irf_pre_processing: events["weight"] = ( @@ -229,10 +276,6 @@ def make_derived_columns( events["true_source_fov_offset"] = calculate_source_fov_offset( events, prefix="true" ) - events["reco_source_fov_offset"] = calculate_source_fov_offset( - events, prefix="reco" - ) - else: reco_icrs = reco.transform_to(ICRS()) events["reco_ra"] = reco_icrs.ra @@ -366,7 +409,7 @@ def load_preselected_events(self, chunk_size: int) -> tuple[QTable, int, dict]: ) events = self.epp.event_selection.calculate_selection(events) selected = events[events["selected"]] - selected = self.epp.keep_nescessary_columns_only(selected) + selected = self.epp.keep_necessary_columns_only(selected) bits.append(selected) bits.append(header) # Putting it last ensures the correct metadata is used From 8dcfea5bd18f1c7db0fb6d33504db1592f81788c Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 6 Mar 2025 15:21:06 +0100 Subject: [PATCH 24/58] Update table creation --- src/ctapipe/conftest.py | 3 +- src/ctapipe/irf/preprocessing.py | 100 +++++++++++++++++++++++++++++-- 2 files changed, 98 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/conftest.py b/src/ctapipe/conftest.py index f4555044f4d..7af605dc5e6 100644 --- a/src/ctapipe/conftest.py +++ b/src/ctapipe/conftest.py @@ -816,7 +816,8 @@ def irf_events_table(): N2 = 100 N = N1 + N2 epp = EventPreprocessor() - tab = epp.make_empty_table() + keep_columns, _, _ = epp.get_columns_keep_rename_scheme(None, True) + tab = epp.make_empty_table(keep_columns) ids, bulk, unitless = tab.colnames[:2], tab.colnames[2:-2], tab.colnames[-2:] diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 6020959f8cd..1c8a4af08a1 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -77,7 +77,7 @@ def __init__( self.event_selection = EventSelection(parent=self) def get_columns_keep_rename_scheme( - self, events: Table, already_derived: bool = False + self, events: Table = None, already_derived: bool = False ): """ Function to get the columns to keep, and the scheme for renaming columns @@ -176,6 +176,10 @@ def get_columns_keep_rename_scheme( "reco_h_max_uncert", ] if not self.raise_error_for_optional: + if events is None: + raise ValueError( + "Require events table to assess existence of optional columns" + ) for i, c in enumerate(rename_from_optional): if c not in events.colnames: self.log.warning( @@ -289,7 +293,7 @@ def make_derived_columns( return events - def make_empty_table(self) -> QTable: + def make_empty_table(self, columns_to_use: list[str]) -> QTable: """ This function defines the columns later functions expect to be present in the event table. @@ -317,6 +321,11 @@ def make_empty_table(self) -> QTable: unit=u.TeV, description="Reconstructed energy", ), + Column( + name="reco_energy_uncert", + unit=u.TeV, + description="Uncertainty on the reconstructed energy", + ), Column( name="reco_az", unit=u.deg, @@ -327,6 +336,31 @@ def make_empty_table(self) -> QTable: unit=u.deg, description="Reconstructed altitude", ), + Column( + name="reco_dir_uncert", + unit=u.deg, + description="Uncertainty on the reconstructed direction", + ), + Column( + name="reco_ra", + unit=u.deg, + description="Reconstructed direction, Right Ascension (ICRS)", + ), + Column( + name="reco_dec", + unit=u.deg, + description="Reconstructed direction, Declination (ICRS)", + ), + Column( + name="reco_glon", + unit=u.deg, + description="Reconstructed direction, galactic longitude", + ), + Column( + name="reco_glat", + unit=u.deg, + description="Reconstructed direction, galactic latitude", + ), Column( name="reco_fov_lat", unit=u.deg, @@ -349,11 +383,21 @@ def make_empty_table(self) -> QTable: unit=u.deg, description="Simulated angular offset from pointing direction", ), + Column( + name="true_source_fov_position_angle", + unit=u.deg, + description="Simulated angular position angle from pointing direction", + ), Column( name="reco_source_fov_offset", unit=u.deg, description="Reconstructed angular offset from pointing direction", ), + Column( + name="reco_source_fov_position_angle", + unit=u.deg, + description="Reconstructed angular position angle from pointing direction", + ), Column( name="gh_score", unit=u.dimensionless_unscaled, @@ -366,9 +410,56 @@ def make_empty_table(self) -> QTable: unit=u.dimensionless_unscaled, description="Event weight", ), + Column( + name="multiplicity", + description="Number of telescopes used for the reconstruction", + ), + Column( + name="reco_core_x", + unit=u.m, + description="Reconstructed position of the core of the shower, x coordinate", + ), + Column( + name="reco_core_y", + unit=u.m, + description="Reconstructed position of the core of the shower, y coordinate", + ), + Column( + name="reco_core_uncert_x", + unit=u.m, + description="Uncertainty on the reconstructed position of the core of the shower, x coordinate", + ), + Column( + name="reco_core_uncert_y", + unit=u.m, + description="Uncertainty on the reconstructed position of the core of the shower, y coordinate", + ), + Column( + name="reco_h_max", + unit=u.m, + description="Reconstructed altitude of the maximum of the shower", + ), + Column( + name="reco_h_max_uncert", + unit=u.m, + description="Uncertainty on the reconstructed altitude of the maximum of the shower", + ), ] - return QTable(columns) + # Rearrange in a dict, easier for searching after + columns_dict = {} + for i in range(len(columns)): + columns_dict[columns[i].name] = columns[i] + + # Select only the necessary columns + columns_for_keep = [] + for c in columns_to_use: + if c in columns_dict.keys(): + columns_for_keep.append(columns_dict[c]) + else: + raise ValueError(f"Missing columns definition for {c}") + + return QTable(columns_for_keep) class EventLoader(Component): @@ -400,7 +491,8 @@ def __init__( def load_preselected_events(self, chunk_size: int) -> tuple[QTable, int, dict]: opts = dict(dl2=True, simulated=True, observation_info=True) with TableLoader(self.file, parent=self, **opts) as load: - header = self.epp.make_empty_table() + keep_columns, _, _ = self.epp.get_columns_keep_rename_scheme(None, True) + header = self.epp.make_empty_table(keep_columns) bits = [header] for _, _, events in load.read_subarray_events_chunked(chunk_size, **opts): events = self.epp.normalise_column_names(events) From 3fcfe225e00c09d3e33ed3b27f4536edd7a57891 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 6 Mar 2025 15:25:01 +0100 Subject: [PATCH 25/58] Correct handling of optional columns for DL3 --- src/ctapipe/irf/preprocessing.py | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 1c8a4af08a1..e20a4d505d0 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -106,14 +106,12 @@ def get_columns_keep_rename_scheme( rename_to = [] keep_columns += ["reco_energy", "reco_az", "reco_alt", "gh_score"] - if self.irf_pre_processing or self.optional_dl3_columns: - keep_columns += ["reco_fov_lat", "reco_fov_lon"] - if self.irf_pre_processing: keep_columns += [ "pointing_az", "pointing_alt", - "theta", + "reco_fov_lat", + "reco_fov_lon" "theta", "true_source_fov_offset", "reco_source_fov_offset", "weight", @@ -121,11 +119,12 @@ def get_columns_keep_rename_scheme( else: keep_columns += ["reco_ra", "reco_dec"] if self.optional_dl3_columns: - keep_columns += [ + keep_columns_optional = [ "multiplicity", "reco_glon", "reco_glat", - "reco_source_fov_offset", + "reco_fov_lat", + "reco_fov_lon" "reco_source_fov_offset", "reco_source_fov_position_angle", "reco_energy_uncert", "reco_dir_uncert", @@ -137,6 +136,21 @@ def get_columns_keep_rename_scheme( "reco_h_max_uncert", ] + if not self.raise_error_for_optional: + if events is None: + raise ValueError( + "Require events table to assess existence of optional columns" + ) + for c in keep_columns_optional: + if c not in events.colnames: + self.log.warning( + f"Optional column {c} is missing from the events table." + ) + else: + keep_columns.append(c) + else: + keep_columns += keep_columns_optional + else: rename_from = [ f"{self.energy_reconstructor}_energy", From 0f795b4abb186736965c99a276cdf6e6428c0b3a Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 6 Mar 2025 15:27:14 +0100 Subject: [PATCH 26/58] Fix on optional columns --- src/ctapipe/irf/preprocessing.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index e20a4d505d0..9f17040274c 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -96,10 +96,6 @@ def get_columns_keep_rename_scheme( keep_columns += [ "time", ] - if self.optional_dl3_columns: - keep_columns += [ - "tels_with_trigger", - ] if already_derived: rename_from = [] @@ -111,7 +107,8 @@ def get_columns_keep_rename_scheme( "pointing_az", "pointing_alt", "reco_fov_lat", - "reco_fov_lon" "theta", + "reco_fov_lon", + "theta", "true_source_fov_offset", "reco_source_fov_offset", "weight", @@ -152,6 +149,11 @@ def get_columns_keep_rename_scheme( keep_columns += keep_columns_optional else: + if self.optional_dl3_columns: + keep_columns += [ + "tels_with_trigger", + ] + rename_from = [ f"{self.energy_reconstructor}_energy", f"{self.geometry_reconstructor}_az", From f3ba31484a85237592a854aa8c94f2340f7ddb18 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 6 Mar 2025 15:51:21 +0100 Subject: [PATCH 27/58] Format column for DL3 format --- src/ctapipe/irf/dl3.py | 62 +++++++++++++++++++++++++++++++++ src/ctapipe/tools/create_dl3.py | 3 +- 2 files changed, 64 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 0e684bd8c50..6d0f0955237 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -1,2 +1,64 @@ +from astropy.table import QTable + +from ctapipe.compat import COPY_IF_NEEDED + + def get_hdu_header_events(): return {"HDUCLASS": "GADF", "HDUCLAS1": "EVENTS"} + + +def transform_events_columns_for_gadf_format(events): + rename_from = ["event_id", "time", "reco_ra", "reco_dec", "reco_energy"] + rename_to = ["EVENT_ID", "TIME", "RA", "DEC", "ENERGY"] + + rename_from_optional = [ + "multiplicity", + "reco_glon", + "reco_glat", + "reco_alt", + "reco_az", + "reco_fov_lon", + "reco_fov_lat", + "reco_source_fov_offset", + "reco_source_fov_position_angle", + "gh_score", + "reco_dir_uncert", + "reco_energy_uncert", + "reco_core_x", + "reco_core_y", + "reco_core_uncert", + "reco_x_max", + "reco_x_max_uncert", + ] + rename_to_optional = [ + "MULTIP", + "GLON", + "GLAT", + "ALT", + "AZ", + "DETX", + "DETY", + "THETA", + "PHI", + "GAMANESS", + "DIR_ERR", + "ENERGY_ERR", + "COREX", + "COREY", + "CORE_ERR", + "XMAX", + "XMAX_ERR", + ] + + for i, c in enumerate(rename_from_optional): + if c not in events.colnames: + pass + # self.log.warning(f"Optional column {c} is missing from the events_table.") + else: + rename_from.append(rename_from_optional[i]) + rename_to.append(rename_to_optional[i]) + + renamed_events = QTable(events, copy=COPY_IF_NEEDED) + renamed_events.rename_columns(rename_from, rename_to) + renamed_events = renamed_events[rename_to] + return renamed_events diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index a44f2e69ddb..466ac44d7c5 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -10,7 +10,7 @@ __all__ = ["DL3Tool"] from ..irf.cuts import EventSelection -from ..irf.dl3 import get_hdu_header_events +from ..irf.dl3 import get_hdu_header_events, transform_events_columns_for_gadf_format class DL3Tool(Tool): @@ -95,6 +95,7 @@ def start(self): parent=self, file=self.dl2_file, quality_selection_only=False ) events = self.event_loader.load_preselected_events(self.chunk_size) + events = transform_events_columns_for_gadf_format(events) hdu_dl3 = fits.HDUList( [ From aedf6d8aa57255bfb0a58be8f90fa2f1a1b7afbf Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 6 Mar 2025 17:39:00 +0100 Subject: [PATCH 28/58] Change x max for h max --- src/ctapipe/irf/dl3.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 6d0f0955237..c580608e4d5 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -27,8 +27,8 @@ def transform_events_columns_for_gadf_format(events): "reco_core_x", "reco_core_y", "reco_core_uncert", - "reco_x_max", - "reco_x_max_uncert", + "reco_h_max", + "reco_h_max_uncert", ] rename_to_optional = [ "MULTIP", @@ -46,8 +46,8 @@ def transform_events_columns_for_gadf_format(events): "COREX", "COREY", "CORE_ERR", - "XMAX", - "XMAX_ERR", + "HMAX", + "HMAX_ERR", ] for i, c in enumerate(rename_from_optional): From e2a119b286b95ebfde5e767b46a35d52e8a5e32d Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 6 Mar 2025 17:51:24 +0100 Subject: [PATCH 29/58] Set base object for DL3 formating --- src/ctapipe/irf/dl3.py | 123 ++++++++++++++++++------------- src/ctapipe/irf/preprocessing.py | 4 +- src/ctapipe/tools/create_dl3.py | 33 +++++++-- 3 files changed, 102 insertions(+), 58 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index c580608e4d5..352c278e7b0 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -1,62 +1,85 @@ from astropy.table import QTable from ctapipe.compat import COPY_IF_NEEDED +from ctapipe.core import Component +from ctapipe.core.traits import Bool -def get_hdu_header_events(): - return {"HDUCLASS": "GADF", "HDUCLAS1": "EVENTS"} +class DL3_GADF(Component): + optional_dl3_columns = Bool( + default_value=False, help="If true add optional columns to produce file" + ).tag(config=False) + raise_error_for_optional = Bool( + default_value=True, + help="If true will raise error in the case optional column are missing", + ).tag(config=False) -def transform_events_columns_for_gadf_format(events): - rename_from = ["event_id", "time", "reco_ra", "reco_dec", "reco_energy"] - rename_to = ["EVENT_ID", "TIME", "RA", "DEC", "ENERGY"] + def __init__(self, **kwargs): + super().__init__(**kwargs) - rename_from_optional = [ - "multiplicity", - "reco_glon", - "reco_glat", - "reco_alt", - "reco_az", - "reco_fov_lon", - "reco_fov_lat", - "reco_source_fov_offset", - "reco_source_fov_position_angle", - "gh_score", - "reco_dir_uncert", - "reco_energy_uncert", - "reco_core_x", - "reco_core_y", - "reco_core_uncert", - "reco_h_max", - "reco_h_max_uncert", - ] - rename_to_optional = [ - "MULTIP", - "GLON", - "GLAT", - "ALT", - "AZ", - "DETX", - "DETY", - "THETA", - "PHI", - "GAMANESS", - "DIR_ERR", - "ENERGY_ERR", - "COREX", - "COREY", - "CORE_ERR", - "HMAX", - "HMAX_ERR", - ] + def get_hdu_header_events(self): + return {"HDUCLASS": "GADF", "HDUCLAS1": "EVENTS"} - for i, c in enumerate(rename_from_optional): - if c not in events.colnames: - pass - # self.log.warning(f"Optional column {c} is missing from the events_table.") - else: - rename_from.append(rename_from_optional[i]) - rename_to.append(rename_to_optional[i]) + def transform_events_columns_for_gadf_format(self, events): + rename_from = ["event_id", "time", "reco_ra", "reco_dec", "reco_energy"] + rename_to = ["EVENT_ID", "TIME", "RA", "DEC", "ENERGY"] + + rename_from_optional = [ + "multiplicity", + "reco_glon", + "reco_glat", + "reco_alt", + "reco_az", + "reco_fov_lon", + "reco_fov_lat", + "reco_source_fov_offset", + "reco_source_fov_position_angle", + "gh_score", + "reco_dir_uncert", + "reco_energy_uncert", + "reco_core_x", + "reco_core_y", + "reco_core_uncert", + "reco_h_max", + "reco_h_max_uncert", + ] + rename_to_optional = [ + "MULTIP", + "GLON", + "GLAT", + "ALT", + "AZ", + "DETX", + "DETY", + "THETA", + "PHI", + "GAMANESS", + "DIR_ERR", + "ENERGY_ERR", + "COREX", + "COREY", + "CORE_ERR", + "HMAX", + "HMAX_ERR", + ] + + if not self.raise_error_for_optional: + for i, c in enumerate(rename_from_optional): + if c not in events.colnames: + self.log.warning( + f"Optional column {c} is missing from the events_table." + ) + else: + rename_from.append(rename_from_optional[i]) + rename_to.append(rename_to_optional[i]) + + for c in rename_from: + if c not in events.colnames: + raise ValueError( + "Input files must conform to the ctapipe DL2 data model. " + f"Required column {c} is missing." + ) renamed_events = QTable(events, copy=COPY_IF_NEEDED) renamed_events.rename_columns(rename_from, rename_to) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 9f17040274c..23e1b0391d8 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -60,12 +60,12 @@ class EventPreprocessor(Component): optional_dl3_columns = Bool( default_value=False, help="If true add optional columns to produce file" - ).tag(config=True) + ).tag(config=False) raise_error_for_optional = Bool( default_value=True, help="If true will raise error in the case optional column are missing", - ).tag(config=True) + ).tag(config=False) def __init__( self, quality_selection_only: bool = True, config=None, parent=None, **kwargs diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index 466ac44d7c5..c028d31da37 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -10,7 +10,7 @@ __all__ = ["DL3Tool"] from ..irf.cuts import EventSelection -from ..irf.dl3 import get_hdu_header_events, transform_events_columns_for_gadf_format +from ..irf.dl3 import DL3_GADF class DL3Tool(Tool): @@ -48,6 +48,15 @@ class DL3Tool(Tool): help="If true, allow to overwrite already existing output file", ).tag(config=True) + optional_dl3_columns = Bool( + default_value=False, help="If true add optional columns to produce file" + ).tag(config=True) + + raise_error_for_optional = Bool( + default_value=True, + help="If true will raise error in the case optional column are missing", + ).tag(config=True) + # Which classes are registered for configuration classes = ( [ @@ -69,13 +78,13 @@ class DL3Tool(Tool): flags = { **flag( "optional-column", - "EventPreprocessor.optional_dl3_columns", + "DL3Tool.optional_dl3_columns", "Add optional columns for events in the DL3 file", "Do not add optional column for events in the DL3 file", ), **flag( "raise-error-for-optional", - "EventPreprocessor.raise_error_for_optional", + "DL3Tool.raise_error_for_optional", "Raise an error if an optional column is missing", "Only display a warning if an optional column is missing, it will lead to optional columns missing in the DL3 file", ), @@ -86,8 +95,16 @@ def setup(self): Initialize components from config and load g/h (and theta) cuts. """ - # Force the preprocessing for DL3 + # Setting preprocessing for DL3 EventPreprocessor.irf_pre_processing = False + EventPreprocessor.optional_dl3_columns = self.optional_dl3_columns + EventPreprocessor.raise_error_for_optional = self.raise_error_for_optional + + # Setting the GADF format object + DL3_GADF.optional_dl3_columns = self.optional_dl3_columns + DL3_GADF.raise_error_for_optional = self.raise_error_for_optional + + self.dl3_format = DL3_GADF() def start(self): self.log.info("Loading events from DL2") @@ -95,7 +112,7 @@ def start(self): parent=self, file=self.dl2_file, quality_selection_only=False ) events = self.event_loader.load_preselected_events(self.chunk_size) - events = transform_events_columns_for_gadf_format(events) + events = self.dl3_format.transform_events_columns_for_gadf_format(events) hdu_dl3 = fits.HDUList( [ @@ -105,7 +122,11 @@ def start(self): ] ) hdu_dl3.append( - fits.BinTableHDU(data=events, name="EVENTS", header=get_hdu_header_events()) + fits.BinTableHDU( + data=events, + name="EVENTS", + header=self.dl3_format.get_hdu_header_events(), + ) ) self.log.info("Loading IRFs") From 475a4bdf5ed1e13d29407d8343376a656b6b2701 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 6 Mar 2025 17:52:07 +0100 Subject: [PATCH 30/58] Update message --- src/ctapipe/irf/dl3.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 352c278e7b0..93e72738e39 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -68,7 +68,7 @@ def transform_events_columns_for_gadf_format(self, events): for i, c in enumerate(rename_from_optional): if c not in events.colnames: self.log.warning( - f"Optional column {c} is missing from the events_table." + f"Optional column {c} is missing from the events table." ) else: rename_from.append(rename_from_optional[i]) @@ -77,8 +77,7 @@ def transform_events_columns_for_gadf_format(self, events): for c in rename_from: if c not in events.colnames: raise ValueError( - "Input files must conform to the ctapipe DL2 data model. " - f"Required column {c} is missing." + f"Required column {c} is missing from the events table." ) renamed_events = QTable(events, copy=COPY_IF_NEEDED) From fa25592dd034f89d64f0f065dd766a2afb61d911 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 6 Mar 2025 17:56:15 +0100 Subject: [PATCH 31/58] Compute average uncertainty on core position --- src/ctapipe/irf/preprocessing.py | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 23e1b0391d8..65ad7145fe8 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -127,8 +127,7 @@ def get_columns_keep_rename_scheme( "reco_dir_uncert", "reco_core_x", "reco_core_y", - "reco_core_uncert_x", - "reco_core_uncert_y", + "reco_core_uncert", "reco_h_max", "reco_h_max_uncert", ] @@ -307,6 +306,11 @@ def make_derived_columns( events["multiplicity"] = np.sum(events["tels_with_trigger"], axis=1) + events["reco_core_uncert"] = np.sqrt( + events["reco_core_uncert_x"] ** 2 + + events["reco_core_uncert_y"] ** 2 + ) + return events def make_empty_table(self, columns_to_use: list[str]) -> QTable: @@ -441,14 +445,9 @@ def make_empty_table(self, columns_to_use: list[str]) -> QTable: description="Reconstructed position of the core of the shower, y coordinate", ), Column( - name="reco_core_uncert_x", - unit=u.m, - description="Uncertainty on the reconstructed position of the core of the shower, x coordinate", - ), - Column( - name="reco_core_uncert_y", + name="reco_core_uncert", unit=u.m, - description="Uncertainty on the reconstructed position of the core of the shower, y coordinate", + description="Uncertainty on the reconstructed position of the core of the shower", ), Column( name="reco_h_max", From 9b0c9f9d1b56b06e96552bf923384999e66b8b4a Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 6 Mar 2025 20:22:22 +0100 Subject: [PATCH 32/58] Switch to fully use the object for DL3 format --- src/ctapipe/io/hdf5eventsource.py | 24 ++++---- src/ctapipe/irf/dl3.py | 77 +++++++++++++++++++++++++ src/ctapipe/tools/create_dl3.py | 33 +++++------ src/ctapipe/tools/tests/test_merge.py | 6 +- src/ctapipe/tools/tests/test_process.py | 12 ++-- 5 files changed, 111 insertions(+), 41 deletions(-) diff --git a/src/ctapipe/io/hdf5eventsource.py b/src/ctapipe/io/hdf5eventsource.py index 8e8cb3ad6b8..779e1e6d6ba 100644 --- a/src/ctapipe/io/hdf5eventsource.py +++ b/src/ctapipe/io/hdf5eventsource.py @@ -340,7 +340,7 @@ def has_simulated_dl1(self): True for files with telescope-wise event information in the simulation group """ if self.is_simulation: - if "telescope" in self.file_.root.simulation.event: + if "telescope" in self.file_.root.simulation.events: return True return False @@ -384,7 +384,7 @@ def simulation_config(self) -> dict[int, SimulationConfigContainer]: return self._simulation_configs def __len__(self): - n_events = len(self.file_.root.dl1.event.subarray.trigger) + n_events = len(self.file_.root.dl1.events.subarray.trigger) if self.max_events is not None: return min(n_events, self.max_events) return n_events @@ -434,7 +434,7 @@ def _parse_sb_and_ob_configs(self): return scheduling_blocks, observation_blocks def _is_hillas_in_camera_frame(self): - parameters_group = self.file_.root.dl1.event.telescope.parameters + parameters_group = self.file_.root.dl1.events.telescope.parameters telescope_tables = parameters_group._v_children.values() # in case of no parameters, it doesn't matter, we just return False @@ -456,7 +456,7 @@ def _generator(self): table.name: self.reader.read( f"/r1/event/telescope/{table.name}", R1CameraContainer ) - for table in self.file_.root.r1.event.telescope + for table in self.file_.root.r1.events.telescope } if DataLevel.DL1_IMAGES in self.datalevels: @@ -472,14 +472,14 @@ def _generator(self): DL1CameraContainer, ignore_columns=ignore_columns, ) - for table in self.file_.root.dl1.event.telescope.images + for table in self.file_.root.dl1.events.telescope.images } if self.has_simulated_dl1: simulated_image_iterators = { - table.name: self.file_.root.simulation.event.telescope.images[ + table.name: self.file_.root.simulation.events.telescope.images[ table.name ].iterrows() - for table in self.file_.root.simulation.event.telescope.images + for table in self.file_.root.simulation.events.telescope.images } if DataLevel.DL1_PARAMETERS in self.datalevels: @@ -516,7 +516,7 @@ def _generator(self): "peak_time", ], ) - for table in self.file_.root.dl1.event.telescope.parameters + for table in self.file_.root.dl1.events.telescope.parameters } if self.has_simulated_dl1: simulated_param_readers = { @@ -537,7 +537,7 @@ def _generator(self): "true_intensity", ], ) - for table in self.file_.root.dl1.event.telescope.parameters + for table in self.file_.root.dl1.events.telescope.parameters } if self.has_muon_parameters: @@ -550,7 +550,7 @@ def _generator(self): MuonEfficiencyContainer, ], ) - for table in self.file_.root.dl1.event.telescope.muon + for table in self.file_.root.dl1.events.telescope.muon } dl2_readers = {} @@ -602,14 +602,14 @@ def _generator(self): SimulatedShowerContainer, prefixes="true", ) - if "impact" in self.file_.root.simulation.event.telescope: + if "impact" in self.file_.root.simulation.events.telescope: true_impact_readers = { table.name: self.reader.read( f"/simulation/event/telescope/impact/{table.name}", containers=TelescopeImpactParameterContainer, prefixes=["true_impact"], ) - for table in self.file_.root.simulation.event.telescope.impact + for table in self.file_.root.simulation.events.telescope.impact } # Setup iterators for the array events diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 93e72738e39..23e3837155a 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -1,3 +1,7 @@ +from datetime import datetime + +from astropy.io import fits +from astropy.io.fits.hdu.base import ExtensionHDU from astropy.table import QTable from ctapipe.compat import COPY_IF_NEEDED @@ -15,8 +19,81 @@ class DL3_GADF(Component): help="If true will raise error in the case optional column are missing", ).tag(config=False) + raise_error_for_missing_hdu = Bool( + default_value=True, + help="If true will raise error if HDU are missing from the final DL3 file", + ).tag(config=False) + + overwrite = Bool( + default_value=False, + help="If true, allow to overwrite already existing output file", + ).tag(config=False) + def __init__(self, **kwargs): super().__init__(**kwargs) + self._events = None + self._pointing = None + self._gti = None + self._aeff = None + self._psf = None + self._edisp = None + self._bkg = None + + def write_file(self, path): + events = self.transform_events_columns_for_gadf_format(self.events) + + hdu_dl3 = fits.HDUList( + [ + fits.PrimaryHDU( + header={"CREATED": datetime.now(tz=datetime.UTC).isoformat()} + ) + ] + ) + hdu_dl3.append( + fits.BinTableHDU( + data=events, + name="EVENTS", + header=self.get_hdu_header_events(), + ) + ) + hdu_dl3.append(self.aeff) + hdu_dl3.append(self.psf) + hdu_dl3.append(self.edisp) + hdu_dl3.append(self.bkg) + + hdu_dl3.writeto(path, checksum=True, overwrite=self.overwrite) + + @property + def events(self): + return self._events + + @events.setter + def events(self, events: QTable): + self._events = events + + @property + def aeff(self): + return self._aeff + + @aeff.setter + def aeff(self, aeff: ExtensionHDU): + self._aeff = aeff + + @property + def psf(self): + return self._psf + + @psf.setter + def psf(self, psf: ExtensionHDU): + self._psf = psf + + @property + def bkg(self): + return self._psf + + @psf.setter + def psf(self, bkg: ExtensionHDU): + self._bkg = bkg def get_hdu_header_events(self): return {"HDUCLASS": "GADF", "HDUCLAS1": "EVENTS"} diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index c028d31da37..1495e1a9b1b 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -1,5 +1,3 @@ -from datetime import datetime - from astropy.io import fits from ctapipe.core import Tool, traits @@ -103,6 +101,7 @@ def setup(self): # Setting the GADF format object DL3_GADF.optional_dl3_columns = self.optional_dl3_columns DL3_GADF.raise_error_for_optional = self.raise_error_for_optional + DL3_GADF.overwrite = self.overwrite self.dl3_format = DL3_GADF() @@ -111,31 +110,25 @@ def start(self): self.event_loader = EventLoader( parent=self, file=self.dl2_file, quality_selection_only=False ) - events = self.event_loader.load_preselected_events(self.chunk_size) - events = self.dl3_format.transform_events_columns_for_gadf_format(events) - - hdu_dl3 = fits.HDUList( - [ - fits.PrimaryHDU( - header={"CREATED": datetime.now(tz=datetime.UTC).isoformat()} - ) - ] - ) - hdu_dl3.append( - fits.BinTableHDU( - data=events, - name="EVENTS", - header=self.dl3_format.get_hdu_header_events(), - ) + self.dl3_format.events = self.event_loader.load_preselected_events( + self.chunk_size ) self.log.info("Loading IRFs") hdu_irfs = fits.open(self.irfs_file, checksum=True) for i in range(1, len(hdu_irfs)): - hdu_dl3.append(hdu_irfs[i]) + if "HDUCLAS2" in hdu_irfs[i].header.keys(): + if hdu_irfs[i].header["HDUCLAS2"] == "EFF_AREA": + self.dl3_format.aeff = hdu_irfs[i] + elif hdu_irfs[i].header["HDUCLAS2"] == "EDISP": + self.dl3_format.edisp = hdu_irfs[i] + elif hdu_irfs[i].header["HDUCLAS2"] == "RPSF": + self.dl3_format.psf = hdu_irfs[i] + elif hdu_irfs[i].header["HDUCLAS2"] == "BKG": + self.dl3_format.bkg = hdu_irfs[i] self.log.info("Writing DL3 File") - hdu_dl3.writeto(self.output_file, checksum=True, overwrite=self.overwrite) + self.dl3_format.write_file(self.output_file) def finish(self): pass diff --git a/src/ctapipe/tools/tests/test_merge.py b/src/ctapipe/tools/tests/test_merge.py index 0879d27f1df..2875ddb88ed 100644 --- a/src/ctapipe/tools/tests/test_merge.py +++ b/src/ctapipe/tools/tests/test_merge.py @@ -103,9 +103,9 @@ def test_skip_images(tmp_path, dl1_file, dl1_proton_file): ) with tables.open_file(output, "r") as f: - assert "images" not in f.root.dl1.event.telescope - assert "images" in f.root.simulation.event.telescope - assert "parameters" in f.root.dl1.event.telescope + assert "images" not in f.root.dl1.events.telescope + assert "images" in f.root.simulation.events.telescope + assert "parameters" in f.root.dl1.events.telescope t = read_table(output, "/simulation/event/telescope/images/tel_001") assert "true_image" not in t.colnames diff --git a/src/ctapipe/tools/tests/test_process.py b/src/ctapipe/tools/tests/test_process.py index 424e10d1b9d..e1fc0eb6588 100644 --- a/src/ctapipe/tools/tests/test_process.py +++ b/src/ctapipe/tools/tests/test_process.py @@ -88,8 +88,8 @@ def test_stage_1_dl1(tmp_path, dl1_image_file, dl1_parameters_file): # check tables were written with tables.open_file(dl1b_from_dl1a_file, mode="r") as testfile: assert testfile.root.dl1 - assert testfile.root.dl1.event.telescope - assert testfile.root.dl1.event.subarray + assert testfile.root.dl1.events.telescope + assert testfile.root.dl1.events.subarray assert testfile.root.configuration.instrument.subarray.layout assert testfile.root.configuration.instrument.telescope.optics assert testfile.root.configuration.instrument.telescope.camera.geometry_0 @@ -226,7 +226,7 @@ def test_stage_2_from_dl1_images(tmp_path, dl1_image_file): # check tables were written with tables.open_file(output, mode="r") as testfile: - assert testfile.root.dl2.event.subarray.geometry.HillasReconstructor + assert testfile.root.dl2.events.subarray.geometry.HillasReconstructor def test_stage_2_from_dl1_params(tmp_path, dl1_parameters_file): @@ -249,7 +249,7 @@ def test_stage_2_from_dl1_params(tmp_path, dl1_parameters_file): # check tables were written with tables.open_file(output, mode="r") as testfile: - assert testfile.root.dl2.event.subarray.geometry.HillasReconstructor + assert testfile.root.dl2.events.subarray.geometry.HillasReconstructor def test_ml_preprocessing_from_simtel(tmp_path): @@ -274,8 +274,8 @@ def test_ml_preprocessing_from_simtel(tmp_path): # check tables were written with tables.open_file(output, mode="r") as testfile: - assert testfile.root.dl1.event.telescope.parameters.tel_002 - assert testfile.root.dl2.event.subarray.geometry.HillasReconstructor + assert testfile.root.dl1.events.telescope.parameters.tel_002 + assert testfile.root.dl2.events.subarray.geometry.HillasReconstructor def test_image_modifications(tmp_path, dl1_image_file): From 8e451237f13891de9cbd847c5ad222e112b6fc8c Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 6 Mar 2025 21:31:26 +0100 Subject: [PATCH 33/58] Fix setter --- src/ctapipe/irf/dl3.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 23e3837155a..a55b849c5e1 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -91,8 +91,8 @@ def psf(self, psf: ExtensionHDU): def bkg(self): return self._psf - @psf.setter - def psf(self, bkg: ExtensionHDU): + @bkg.setter + def bkg(self, bkg: ExtensionHDU): self._bkg = bkg def get_hdu_header_events(self): From 5c06050478cb415f026b04acb26de0f8b7b90718 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 6 Mar 2025 21:40:45 +0100 Subject: [PATCH 34/58] Fix getter --- src/ctapipe/irf/dl3.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index a55b849c5e1..2d0d5574f28 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -89,7 +89,7 @@ def psf(self, psf: ExtensionHDU): @property def bkg(self): - return self._psf + return self._bkg @bkg.setter def bkg(self, bkg: ExtensionHDU): From e06054fa7261e8a7c0e3333cd195d164de8fed05 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 6 Mar 2025 21:43:55 +0100 Subject: [PATCH 35/58] Add edisp getter and setter --- src/ctapipe/irf/dl3.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 2d0d5574f28..04217462f6b 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -87,6 +87,14 @@ def psf(self): def psf(self, psf: ExtensionHDU): self._psf = psf + @property + def edisp(self): + return self._edisp + + @edisp.setter + def edisp(self, edisp: ExtensionHDU): + self._edisp = edisp + @property def bkg(self): return self._bkg From 0a1d94b4b73826d6c646ad43a888df81e8407aaa Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Mon, 10 Mar 2025 16:31:55 +0100 Subject: [PATCH 36/58] Direct access for overwrite --- src/ctapipe/irf/dl3.py | 2 +- src/ctapipe/tools/create_dl3.py | 12 ++++++------ 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 04217462f6b..613fdcaf4b8 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -27,7 +27,7 @@ class DL3_GADF(Component): overwrite = Bool( default_value=False, help="If true, allow to overwrite already existing output file", - ).tag(config=False) + ).tag(config=True) def __init__(self, **kwargs): super().__init__(**kwargs) diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index 1495e1a9b1b..51c14600f26 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -41,11 +41,6 @@ class DL3Tool(Tool): help="How many subarray events to load at once while selecting.", ).tag(config=True) - overwrite = Bool( - default_value=False, - help="If true, allow to overwrite already existing output file", - ).tag(config=True) - optional_dl3_columns = Bool( default_value=False, help="If true add optional columns to produce file" ).tag(config=True) @@ -70,10 +65,15 @@ class DL3Tool(Tool): "irfs-file": "DL3Tool.irfs_file", "output": "DL3Tool.output_file", "chunk-size": "DL3Tool.chunk_size", - "overwrite": "DL3Tool.overwrite", } flags = { + **flag( + "overwrite", + "DL3_GADF.overwrite", + "Will allow to overwrite existing DL3 file", + "Prevent overwriting existing DL3 file", + ), **flag( "optional-column", "DL3Tool.optional_dl3_columns", From 60f14aa9d6962a6b4062db46011b508cb7ea4290 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Mon, 10 Mar 2025 16:42:24 +0100 Subject: [PATCH 37/58] Add a base abstract class for generic DL3 format --- src/ctapipe/irf/dl3.py | 98 +++++++++++++++++++++++---------- src/ctapipe/tools/create_dl3.py | 10 ++-- 2 files changed, 74 insertions(+), 34 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 613fdcaf4b8..3372c7b53a0 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -1,5 +1,7 @@ +from abc import abstractmethod from datetime import datetime +from astropy.coordinates import EarthLocation from astropy.io import fits from astropy.io.fits.hdu.base import ExtensionHDU from astropy.table import QTable @@ -9,7 +11,12 @@ from ctapipe.core.traits import Bool -class DL3_GADF(Component): +class DL3_Format(Component): + overwrite = Bool( + default_value=False, + help="If true, allow to overwrite already existing output file", + ).tag(config=True) + optional_dl3_columns = Bool( default_value=False, help="If true add optional columns to produce file" ).tag(config=False) @@ -24,11 +31,6 @@ class DL3_GADF(Component): help="If true will raise error if HDU are missing from the final DL3 file", ).tag(config=False) - overwrite = Bool( - default_value=False, - help="If true, allow to overwrite already existing output file", - ).tag(config=True) - def __init__(self, **kwargs): super().__init__(**kwargs) self._events = None @@ -38,30 +40,11 @@ def __init__(self, **kwargs): self._psf = None self._edisp = None self._bkg = None + self._location = None + @abstractmethod def write_file(self, path): - events = self.transform_events_columns_for_gadf_format(self.events) - - hdu_dl3 = fits.HDUList( - [ - fits.PrimaryHDU( - header={"CREATED": datetime.now(tz=datetime.UTC).isoformat()} - ) - ] - ) - hdu_dl3.append( - fits.BinTableHDU( - data=events, - name="EVENTS", - header=self.get_hdu_header_events(), - ) - ) - hdu_dl3.append(self.aeff) - hdu_dl3.append(self.psf) - hdu_dl3.append(self.edisp) - hdu_dl3.append(self.bkg) - - hdu_dl3.writeto(path, checksum=True, overwrite=self.overwrite) + pass @property def events(self): @@ -69,6 +52,10 @@ def events(self): @events.setter def events(self, events: QTable): + if self._events is not None: + self.log.warning( + "Events table for DL3 file was already set, replacing current event table" + ) self._events = events @property @@ -77,6 +64,10 @@ def aeff(self): @aeff.setter def aeff(self, aeff: ExtensionHDU): + if self._aeff is not None: + self.log.warning( + "Effective area for DL3 file was already set, replacing current effective area" + ) self._aeff = aeff @property @@ -85,6 +76,8 @@ def psf(self): @psf.setter def psf(self, psf: ExtensionHDU): + if self._psf is not None: + self.log.warning("PSF for DL3 file was already set, replacing current PSF") self._psf = psf @property @@ -93,6 +86,10 @@ def edisp(self): @edisp.setter def edisp(self, edisp: ExtensionHDU): + if self._edisp is not None: + self.log.warning( + "EDISP for DL3 file was already set, replacing current EDISP" + ) self._edisp = edisp @property @@ -101,8 +98,53 @@ def bkg(self): @bkg.setter def bkg(self, bkg: ExtensionHDU): + if self._bkg is not None: + self.log.warning( + "Background for DL3 file was already set, replacing current background" + ) self._bkg = bkg + @property + def location(self): + return self._location + + @location.setter + def location(self, location: EarthLocation): + if self._location is not None: + self.log.warning( + "Telescope location table for DL3 file was already set, replacing current location" + ) + self._location = location + + +class DL3_GADF(DL3_Format): + def __init__(self, **kwargs): + super().__init__(**kwargs) + + def write_file(self, path): + events = self.transform_events_columns_for_gadf_format(self.events) + + hdu_dl3 = fits.HDUList( + [ + fits.PrimaryHDU( + header={"CREATED": datetime.now(tz=datetime.UTC).isoformat()} + ) + ] + ) + hdu_dl3.append( + fits.BinTableHDU( + data=events, + name="EVENTS", + header=self.get_hdu_header_events(), + ) + ) + hdu_dl3.append(self.aeff) + hdu_dl3.append(self.psf) + hdu_dl3.append(self.edisp) + hdu_dl3.append(self.bkg) + + hdu_dl3.writeto(path, checksum=True, overwrite=self.overwrite) + def get_hdu_header_events(self): return {"HDUCLASS": "GADF", "HDUCLAS1": "EVENTS"} diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index 51c14600f26..550e7d2e50d 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -4,12 +4,11 @@ from ctapipe.core.traits import Bool, Integer, classes_with_traits, flag from ..irf import EventLoader, EventPreprocessor +from ..irf.cuts import EventSelection +from ..irf.dl3 import DL3_GADF, DL3_Format __all__ = ["DL3Tool"] -from ..irf.cuts import EventSelection -from ..irf.dl3 import DL3_GADF - class DL3Tool(Tool): name = "ctapipe-create-dl3" @@ -99,9 +98,8 @@ def setup(self): EventPreprocessor.raise_error_for_optional = self.raise_error_for_optional # Setting the GADF format object - DL3_GADF.optional_dl3_columns = self.optional_dl3_columns - DL3_GADF.raise_error_for_optional = self.raise_error_for_optional - DL3_GADF.overwrite = self.overwrite + DL3_Format.optional_dl3_columns = self.optional_dl3_columns + DL3_Format.raise_error_for_optional = self.raise_error_for_optional self.dl3_format = DL3_GADF() From d2d87ff706d2f3d1189392f1f3d878bfec8c85de Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Mon, 10 Mar 2025 17:48:40 +0100 Subject: [PATCH 38/58] Add loading of location and gti --- src/ctapipe/irf/dl3.py | 26 +++++++++++++++------ src/ctapipe/irf/preprocessing.py | 39 ++++++++++++++++++++++++++++---- src/ctapipe/tools/create_dl3.py | 3 +++ 3 files changed, 56 insertions(+), 12 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 3372c7b53a0..f543d8bcfb3 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -1,10 +1,12 @@ from abc import abstractmethod from datetime import datetime +from typing import List, Tuple from astropy.coordinates import EarthLocation from astropy.io import fits from astropy.io.fits.hdu.base import ExtensionHDU from astropy.table import QTable +from astropy.time import Time from ctapipe.compat import COPY_IF_NEEDED from ctapipe.core import Component @@ -47,7 +49,7 @@ def write_file(self, path): pass @property - def events(self): + def events(self) -> QTable: return self._events @events.setter @@ -59,7 +61,7 @@ def events(self, events: QTable): self._events = events @property - def aeff(self): + def aeff(self) -> ExtensionHDU: return self._aeff @aeff.setter @@ -71,7 +73,7 @@ def aeff(self, aeff: ExtensionHDU): self._aeff = aeff @property - def psf(self): + def psf(self) -> ExtensionHDU: return self._psf @psf.setter @@ -81,7 +83,7 @@ def psf(self, psf: ExtensionHDU): self._psf = psf @property - def edisp(self): + def edisp(self) -> ExtensionHDU: return self._edisp @edisp.setter @@ -93,7 +95,7 @@ def edisp(self, edisp: ExtensionHDU): self._edisp = edisp @property - def bkg(self): + def bkg(self) -> ExtensionHDU: return self._bkg @bkg.setter @@ -105,17 +107,27 @@ def bkg(self, bkg: ExtensionHDU): self._bkg = bkg @property - def location(self): + def location(self) -> EarthLocation: return self._location @location.setter def location(self, location: EarthLocation): if self._location is not None: self.log.warning( - "Telescope location table for DL3 file was already set, replacing current location" + "Telescope location for DL3 file was already set, replacing current location" ) self._location = location + @property + def gti(self) -> List[Tuple[Time, Time]]: + return self._gti + + @gti.setter + def gti(self, gti: List[Tuple[Time, Time]]): + if self._gti is not None: + self.log.warning("GTi for DL3 file was already set, replacing current gti") + self._gti = gti + class DL3_GADF(DL3_Format): def __init__(self, **kwargs): diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 65ad7145fe8..cf82401febc 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -1,6 +1,7 @@ """Module containing classes related to event loading and preprocessing""" from pathlib import Path +from typing import List, Tuple import astropy.units as u import numpy as np @@ -502,14 +503,16 @@ def __init__( else: self.target_spectrum = None self.file = file + self.opts_loader = dict(dl2=True, simulated=True, observation_info=True) def load_preselected_events(self, chunk_size: int) -> tuple[QTable, int, dict]: - opts = dict(dl2=True, simulated=True, observation_info=True) - with TableLoader(self.file, parent=self, **opts) as load: + with TableLoader(self.file, parent=self, **self.opts_loader) as load: keep_columns, _, _ = self.epp.get_columns_keep_rename_scheme(None, True) header = self.epp.make_empty_table(keep_columns) bits = [header] - for _, _, events in load.read_subarray_events_chunked(chunk_size, **opts): + for _, _, events in load.read_subarray_events_chunked( + chunk_size, **self.opts_loader + ): events = self.epp.normalise_column_names(events) events = self.epp.make_derived_columns( events, load.subarray.reference_location @@ -523,11 +526,37 @@ def load_preselected_events(self, chunk_size: int) -> tuple[QTable, int, dict]: table = vstack(bits, join_type="exact", metadata_conflicts="silent") return table + def get_observation_information( + self, + ) -> Tuple[EarthLocation, List[Tuple[Time, Time]]]: + with TableLoader(self.file, parent=self, **self.opts_loader) as load: + array_location = load.subarray.reference_location + + # TODO : proper retrieval of gti information + start_time = None + stop_time = None + for _, _, events in load.read_subarray_events_chunked( + 10000, **self.opts_loader + ): + start_time = ( + Time(np.min(events["time"])) + if start_time is None + else min(start_time, Time(np.min(events["time"]))) + ) + stop_time = ( + Time(np.max(events["time"])) + if start_time is None + else max(stop_time, Time(np.max(events["time"]))) + ) + + return array_location, [ + (start_time, stop_time), + ] + def get_simulation_information( self, obs_time: u.Quantity ) -> tuple[SimulatedEventsInfo, PowerLaw]: - opts = dict(dl2=True, simulated=True, observation_info=True) - with TableLoader(self.file, parent=self, **opts) as load: + with TableLoader(self.file, parent=self, **self.opts_loader) as load: sim = load.read_simulation_configuration() try: show = load.read_shower_distribution() diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index 550e7d2e50d..defd0def2b8 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -111,6 +111,9 @@ def start(self): self.dl3_format.events = self.event_loader.load_preselected_events( self.chunk_size ) + array_location, gti = self.event_loader.get_observation_information() + self.dl3_format.location = array_location + self.dl3_format.gti = gti self.log.info("Loading IRFs") hdu_irfs = fits.open(self.irfs_file, checksum=True) From 663faf403016d9178cc12c76c981bfa7c5597730 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Mon, 10 Mar 2025 18:02:14 +0100 Subject: [PATCH 39/58] Fix --- src/ctapipe/irf/preprocessing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index cf82401febc..2162cbe3822 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -545,7 +545,7 @@ def get_observation_information( ) stop_time = ( Time(np.max(events["time"])) - if start_time is None + if stop_time is None else max(stop_time, Time(np.max(events["time"]))) ) From 0f712c6dc04009142d40ecdaadd951712604959a Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Tue, 11 Mar 2025 10:50:59 +0100 Subject: [PATCH 40/58] Update headers of the DL3 file --- src/ctapipe/irf/dl3.py | 73 +++++++++++++++++++++++++++++++++++++----- 1 file changed, 65 insertions(+), 8 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index f543d8bcfb3..c818628f257 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -2,15 +2,18 @@ from datetime import datetime from typing import List, Tuple +import astropy.units as u from astropy.coordinates import EarthLocation from astropy.io import fits +from astropy.io.fits import Header from astropy.io.fits.hdu.base import ExtensionHDU from astropy.table import QTable -from astropy.time import Time +from astropy.time import Time, TimeDelta from ctapipe.compat import COPY_IF_NEEDED from ctapipe.core import Component from ctapipe.core.traits import Bool +from ctapipe.version import version as ctapipe_version class DL3_Format(Component): @@ -132,22 +135,21 @@ def gti(self, gti: List[Tuple[Time, Time]]): class DL3_GADF(DL3_Format): def __init__(self, **kwargs): super().__init__(**kwargs) + self.file_creation_time = datetime.now(tz=datetime.UTC).isoformat() def write_file(self, path): + self.file_creation_time = datetime.now(tz=datetime.UTC).isoformat() + events = self.transform_events_columns_for_gadf_format(self.events) hdu_dl3 = fits.HDUList( - [ - fits.PrimaryHDU( - header={"CREATED": datetime.now(tz=datetime.UTC).isoformat()} - ) - ] + [fits.PrimaryHDU(header=Header(self.get_hdu_header_base_format()))] ) hdu_dl3.append( fits.BinTableHDU( data=events, name="EVENTS", - header=self.get_hdu_header_events(), + header=Header(self.get_hdu_header_events()), ) ) hdu_dl3.append(self.aeff) @@ -157,8 +159,63 @@ def write_file(self, path): hdu_dl3.writeto(path, checksum=True, overwrite=self.overwrite) + def get_hdu_header_base_format(self): + return { + "HDUCLASS": "GADF", + "HDUVERS": "v0.3", + "HDUDOC": "https://gamma-astro-data-formats.readthedocs.io/en/v0.3/index.html", + "CREATOR": "ctapipe " + ctapipe_version, + "CREATED": self.file_creation_time, + } + + def get_hdu_header_base_time(self): + if self._gti is None: + raise ValueError("No available time information for the DL3 file") + start_time = None + stop_time = None + ontime = TimeDelta(0.0 * u.s) + for gti_interval in self._gti: + ontime += gti_interval[1] - gti_interval[0] + start_time = ( + gti_interval[0] + if start_time is None + else min(start_time, gti_interval[0]) + ) + stop_time = ( + gti_interval[1] + if stop_time is None + else max(stop_time, gti_interval[1]) + ) + + reference_time = Time(datetime.fromisoformat("1970-01-01T00:00:00+00:00")) + + return { + "MJDREFI": int(reference_time.mjd), + "MJDREFF": reference_time.mjd % 1, + "TIMEUNIT": "s", + "TIMEREF": "GEOCENTER", + "TIMESYS": "UTC", + "TSTART": start_time, + "TSTOP": stop_time, + "ONTIME": ontime.to_value(u.s), + "TELAPSE": (stop_time - start_time).to_value(u.s), + "DATE-OBS": start_time.fits, + "DATE-BEG": start_time.fits, + "DATE-AVG": (start_time + (stop_time - start_time) / 2.0).fits, + "DATE-END": stop_time.fits, + } + def get_hdu_header_events(self): - return {"HDUCLASS": "GADF", "HDUCLAS1": "EVENTS"} + header = self.get_hdu_header_base_format() + header.update({"HDUCLAS1": "EVENTS"}) + header.update(self.get_hdu_header_base_time()) + return header + + def get_hdu_header_gti(self): + header = self.get_hdu_header_base_format() + header.update({"HDUCLAS1": "GTI"}) + header.update(self.get_hdu_header_base_time()) + return header def transform_events_columns_for_gadf_format(self, events): rename_from = ["event_id", "time", "reco_ra", "reco_dec", "reco_energy"] From a9ca20cfb48c8770483fc7de63428fe9196ba362 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Tue, 11 Mar 2025 13:29:23 +0100 Subject: [PATCH 41/58] Update way of reading metainformation of the runs --- src/ctapipe/irf/preprocessing.py | 47 ++++++++++++++++---------------- src/ctapipe/tools/create_dl3.py | 6 ++-- 2 files changed, 26 insertions(+), 27 deletions(-) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 2162cbe3822..e0f44b18c7b 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -1,7 +1,7 @@ """Module containing classes related to event loading and preprocessing""" from pathlib import Path -from typing import List, Tuple +from typing import Any, Dict import astropy.units as u import numpy as np @@ -526,32 +526,31 @@ def load_preselected_events(self, chunk_size: int) -> tuple[QTable, int, dict]: table = vstack(bits, join_type="exact", metadata_conflicts="silent") return table - def get_observation_information( - self, - ) -> Tuple[EarthLocation, List[Tuple[Time, Time]]]: + def get_observation_information(self) -> Dict[str, Any]: with TableLoader(self.file, parent=self, **self.opts_loader) as load: - array_location = load.subarray.reference_location - - # TODO : proper retrieval of gti information - start_time = None - stop_time = None - for _, _, events in load.read_subarray_events_chunked( - 10000, **self.opts_loader - ): - start_time = ( - Time(np.min(events["time"])) - if start_time is None - else min(start_time, Time(np.min(events["time"]))) - ) - stop_time = ( - Time(np.max(events["time"])) - if stop_time is None - else max(stop_time, Time(np.max(events["time"]))) + # Extract telescope location + meta = {} + meta["location"] = load.subarray.reference_location + + obs_info_table = load.read_observation_information() + if len(obs_info_table) == 0: + self.log.error("No observation information available in the DL2 file") + raise ValueError("Missing observation information in the DL2 file") + + if len(np.unique(obs_info_table["obs_id"])): + self.log.warning( + "Multiple observations included in the DL2 file, only the id of the first observation will be reported in the DL3 file, data from all observations will be included" ) + meta["obs_id"] = np.unique(obs_info_table["obs_id"])[0] - return array_location, [ - (start_time, stop_time), - ] + # Extract GTI + list_gti = [] + for i in range(len(obs_info_table)): + start_time = obs_info_table["actual_start_time"][i] + stop_time = start_time + obs_info_table["actual_duration"][i] + list_gti.append((start_time, stop_time)) + + return meta def get_simulation_information( self, obs_time: u.Quantity diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index defd0def2b8..ffa33a4f692 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -111,9 +111,9 @@ def start(self): self.dl3_format.events = self.event_loader.load_preselected_events( self.chunk_size ) - array_location, gti = self.event_loader.get_observation_information() - self.dl3_format.location = array_location - self.dl3_format.gti = gti + meta = self.event_loader.get_observation_information() + self.dl3_format.location = meta["location"] + self.dl3_format.gti = meta["gti"] self.log.info("Loading IRFs") hdu_irfs = fits.open(self.irfs_file, checksum=True) From 74cca7f57544c82e269e1cc5de75908b7dd7cd09 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Wed, 12 Mar 2025 12:08:03 +0100 Subject: [PATCH 42/58] Add reading of additional observation metadata from dl2 file --- src/ctapipe/irf/preprocessing.py | 120 ++++++++++++++++++++++++++++--- 1 file changed, 112 insertions(+), 8 deletions(-) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index e0f44b18c7b..30f7c079da6 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -7,7 +7,7 @@ import numpy as np from astropy.coordinates import ICRS, AltAz, EarthLocation, Galactic, SkyCoord from astropy.table import Column, QTable, Table, vstack -from astropy.time import Time +from astropy.time import Time, TimeDelta from pyirf.simulations import SimulatedEventsInfo from pyirf.spectral import ( DIFFUSE_FLUX_UNIT, @@ -23,16 +23,19 @@ from tables import NoSuchNodeError from ..compat import COPY_IF_NEEDED -from ..containers import CoordinateFrameType +from ..containers import CoordinateFrameType, PointingMode from ..coordinates import NominalFrame from ..core import Component from ..core.traits import Bool, Unicode from ..io import TableLoader +from ..version import version as ctapipe_version from .cuts import EventQualitySelection, EventSelection from .spectra import SPECTRA, Spectra __all__ = ["EventLoader", "EventPreprocessor"] +from ..io.astropy_helpers import join_allow_empty + class EventPreprocessor(Component): """Defines pre-selection cuts and the necessary renaming of columns.""" @@ -532,23 +535,124 @@ def get_observation_information(self) -> Dict[str, Any]: meta = {} meta["location"] = load.subarray.reference_location + # Read observation information obs_info_table = load.read_observation_information() - if len(obs_info_table) == 0: + schedule_info_table = load.read_scheduling_blocks() + obs_all_info_table = join_allow_empty( + obs_info_table, schedule_info_table, "sb_id", "inner" + ) + if len(obs_all_info_table) == 0: self.log.error("No observation information available in the DL2 file") raise ValueError("Missing observation information in the DL2 file") + elif len(obs_all_info_table) != len(obs_info_table): + self.log.error( + "Scheduling blocks are missing for some observation block" + ) + raise ValueError( + "Scheduling blocks are missing for some observation block" + ) - if len(np.unique(obs_info_table["obs_id"])): + # Check for obs ids + if len(np.unique(obs_all_info_table["obs_id"])): self.log.warning( "Multiple observations included in the DL2 file, only the id of the first observation will be reported in the DL3 file, data from all observations will be included" ) - meta["obs_id"] = np.unique(obs_info_table["obs_id"])[0] + meta["obs_id"] = np.unique(obs_all_info_table["obs_id"])[0] # Extract GTI list_gti = [] - for i in range(len(obs_info_table)): - start_time = obs_info_table["actual_start_time"][i] - stop_time = start_time + obs_info_table["actual_duration"][i] + for i in range(len(obs_all_info_table)): + start_time = Time(obs_all_info_table["actual_start_time"][i]) + stop_time = start_time + TimeDelta( + obs_all_info_table["actual_duration"][i] + ) list_gti.append((start_time, stop_time)) + meta["gti"] = list_gti + + # Extract pointing information + if len(np.unique(obs_all_info_table["pointing_mode"])) != 1: + self.log.error("Multiple pointing mode are not supported") + raise ValueError("Multiple pointing mode are not supported") + meta["pointing"] = {} + meta["pointing"]["pointing_mode"] = PointingMode( + obs_all_info_table["pointing_mode"][0] + ).name + list_pointing = [] + for i in range(len(obs_all_info_table)): + if ( + CoordinateFrameType( + obs_all_info_table["subarray_pointing_frame"][i] + ).name + == "ALTAZ" + ): + pointing_start = AltAz( + alt=obs_all_info_table["subarray_pointing_lat"][i], + az=obs_all_info_table["subarray_pointing_lon"][i], + location=meta["location"], + obstime=meta["gti"][i][0], + ) + pointing_stop = AltAz( + alt=obs_all_info_table["subarray_pointing_lat"][i], + az=obs_all_info_table["subarray_pointing_lon"][i], + location=meta["location"], + obstime=meta["gti"][i][1], + ) + elif ( + CoordinateFrameType( + obs_all_info_table["subarray_pointing_frame"][i] + ).name + == "ICRS" + ): + pointing_start = ICRS( + dec=obs_all_info_table["subarray_pointing_lat"][i], + ra=obs_all_info_table["subarray_pointing_lon"][i], + ) + pointing_stop = pointing_start + elif ( + CoordinateFrameType( + obs_all_info_table["subarray_pointing_frame"][i] + ).name + == "GALACTIC" + ): + pointing_start = Galactic( + b=obs_all_info_table["subarray_pointing_lat"][i], + l=obs_all_info_table["subarray_pointing_lon"][i], + ) + pointing_stop = pointing_start + list_pointing.append((meta["gti"][i][0], pointing_start)) + list_pointing.append((meta["gti"][i][1], pointing_stop)) + meta["pointing"]["pointing_list"] = list_pointing + + # Telescope information + self.log.warning( + "Name of organisation, array and subarray are not read from the file and are default value" + ) + meta["telescope_information"] = { + "organisation": "UNKNOWN", + "array": "UNKNOWN", + "subarray": "UNKNOWN", + "telescope_list": np.array( + load.subarray.get_tel_ids(load.subarray.tel) + ), + } + + # Target information + self.log.warning( + "Name of the target, coordinates and observer are not read from the file and are default value" + ) + meta["target"] = { + "observer": "UNKNOWN", + "object": "UNKNOWN", + "object_coordinate": ICRS(np.nan * u.deg, np.nan * u.deg), + } + + # Software information + self.log.warning("Software version are unknown") + meta["software version"] = { + "analysis_version": "ctapipe " + ctapipe_version, + "calibration_version": "UNKNOWN", + "dst_version": "UNKNOWN", + } return meta From 681ce64395e0d5354e6d9b53ff499876372092d4 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Wed, 12 Mar 2025 13:34:02 +0100 Subject: [PATCH 43/58] Add new metadata to dl3 class --- src/ctapipe/irf/dl3.py | 69 +++++++++++++++++++++++++++----- src/ctapipe/irf/preprocessing.py | 2 +- src/ctapipe/tools/create_dl3.py | 7 +++- 3 files changed, 67 insertions(+), 11 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index c818628f257..766b61b99b8 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -1,9 +1,9 @@ from abc import abstractmethod from datetime import datetime -from typing import List, Tuple +from typing import Any, Dict, List, Tuple import astropy.units as u -from astropy.coordinates import EarthLocation +from astropy.coordinates import EarthLocation, SkyCoord from astropy.io import fits from astropy.io.fits import Header from astropy.io.fits.hdu.base import ExtensionHDU @@ -46,6 +46,9 @@ def __init__(self, **kwargs): self._edisp = None self._bkg = None self._location = None + self._telescope_information = None + self._target_information = None + self._software_information = None @abstractmethod def write_file(self, path): @@ -63,6 +66,28 @@ def events(self, events: QTable): ) self._events = events + @property + def pointing(self) -> List[Tuple[Time, SkyCoord]]: + return self._pointing + + @pointing.setter + def pointing(self, pointing: List[Tuple[Time, SkyCoord]]): + if self._pointing is not None: + self.log.warning( + "Pointing for DL3 file was already set, replacing current pointing" + ) + self._pointing = pointing + + @property + def gti(self) -> List[Tuple[Time, Time]]: + return self._gti + + @gti.setter + def gti(self, gti: List[Tuple[Time, Time]]): + if self._gti is not None: + self.log.warning("GTI for DL3 file was already set, replacing current gti") + self._gti = gti + @property def aeff(self) -> ExtensionHDU: return self._aeff @@ -122,14 +147,40 @@ def location(self, location: EarthLocation): self._location = location @property - def gti(self) -> List[Tuple[Time, Time]]: - return self._gti + def telescope_information(self) -> Dict[str, Any]: + return self._telescope_information - @gti.setter - def gti(self, gti: List[Tuple[Time, Time]]): - if self._gti is not None: - self.log.warning("GTi for DL3 file was already set, replacing current gti") - self._gti = gti + @telescope_information.setter + def telescope_information(self, telescope_information: Dict[str, Any]): + if self._telescope_information is not None: + self.log.warning( + "Telescope information for DL3 file was already set, replacing current information" + ) + self._telescope_information = telescope_information + + @property + def target_information(self) -> Dict[str, Any]: + return self._target_information + + @target_information.setter + def target_information(self, target_information: Dict[str, Any]): + if self._target_information is not None: + self.log.warning( + "Target information for DL3 file was already set, replacing current target information" + ) + self._target_information = target_information + + @property + def software_information(self) -> Dict[str, Any]: + return self._software_information + + @software_information.setter + def software_information(self, software_information: Dict[str, Any]): + if self._software_information is not None: + self.log.warning( + "Software information for DL3 file was already set, replacing current software information" + ) + self._software_information = software_information class DL3_GADF(DL3_Format): diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 30f7c079da6..5075e0d04b4 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -648,7 +648,7 @@ def get_observation_information(self) -> Dict[str, Any]: # Software information self.log.warning("Software version are unknown") - meta["software version"] = { + meta["software_version"] = { "analysis_version": "ctapipe " + ctapipe_version, "calibration_version": "UNKNOWN", "dst_version": "UNKNOWN", diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index ffa33a4f692..bcf6b583213 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -112,9 +112,14 @@ def start(self): self.chunk_size ) meta = self.event_loader.get_observation_information() - self.dl3_format.location = meta["location"] + self.dl3_format.pointing = meta["pointing"] self.dl3_format.gti = meta["gti"] + self.dl3_format.location = meta["location"] + self.dl3_format.telescope_information = meta["telescope_information"] + self.dl3_format.target_information = meta["target"] + self.dl3_format.software_information = meta["software_version"] + self.log.info("Loading IRFs") hdu_irfs = fits.open(self.irfs_file, checksum=True) for i in range(1, len(hdu_irfs)): From 1ce26f5d6ae73b0586b4dce4526ff1dd234d1baf Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Wed, 12 Mar 2025 18:27:16 +0100 Subject: [PATCH 44/58] Add dead time fraction as metadata --- src/ctapipe/irf/dl3.py | 13 +++++++++++++ src/ctapipe/irf/preprocessing.py | 6 ++++++ src/ctapipe/tools/create_dl3.py | 1 + 3 files changed, 20 insertions(+) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 766b61b99b8..451d1bf09ba 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -45,6 +45,7 @@ def __init__(self, **kwargs): self._psf = None self._edisp = None self._bkg = None + self._dead_time_fraction = None self._location = None self._telescope_information = None self._target_information = None @@ -146,6 +147,18 @@ def location(self, location: EarthLocation): ) self._location = location + @property + def dead_time_fraction(self) -> float: + return self._dead_time_fraction + + @dead_time_fraction.setter + def dead_time_fraction(self, dead_time_fraction: float): + if self.dead_time_fraction is not None: + self.log.warning( + "Dead time fraction for DL3 file was already set, replacing current dead time fraction" + ) + self._dead_time_fraction = dead_time_fraction + @property def telescope_information(self) -> Dict[str, Any]: return self._telescope_information diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 5075e0d04b4..034499adfbb 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -569,6 +569,12 @@ def get_observation_information(self) -> Dict[str, Any]: list_gti.append((start_time, stop_time)) meta["gti"] = list_gti + # Dead time fraction + self.log.warning( + "Dead time fraction is not read from the file, a default value of 1 is used" + ) + meta["dead_time_fraction"] = 1.0 + # Extract pointing information if len(np.unique(obs_all_info_table["pointing_mode"])) != 1: self.log.error("Multiple pointing mode are not supported") diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index bcf6b583213..2459262024d 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -114,6 +114,7 @@ def start(self): meta = self.event_loader.get_observation_information() self.dl3_format.pointing = meta["pointing"] self.dl3_format.gti = meta["gti"] + self.dl3_format.dead_time_fraction = meta["dead_time_fraction"] self.dl3_format.location = meta["location"] self.dl3_format.telescope_information = meta["telescope_information"] From 25f8f8bb44dfa7ec28c9a515fe28ccc254f3b836 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Wed, 12 Mar 2025 18:28:17 +0100 Subject: [PATCH 45/58] Proper selection of mandatory only columns --- src/ctapipe/irf/dl3.py | 97 +++++++++++++++++++++--------------------- 1 file changed, 49 insertions(+), 48 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 451d1bf09ba..6b2102b346d 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -285,54 +285,55 @@ def transform_events_columns_for_gadf_format(self, events): rename_from = ["event_id", "time", "reco_ra", "reco_dec", "reco_energy"] rename_to = ["EVENT_ID", "TIME", "RA", "DEC", "ENERGY"] - rename_from_optional = [ - "multiplicity", - "reco_glon", - "reco_glat", - "reco_alt", - "reco_az", - "reco_fov_lon", - "reco_fov_lat", - "reco_source_fov_offset", - "reco_source_fov_position_angle", - "gh_score", - "reco_dir_uncert", - "reco_energy_uncert", - "reco_core_x", - "reco_core_y", - "reco_core_uncert", - "reco_h_max", - "reco_h_max_uncert", - ] - rename_to_optional = [ - "MULTIP", - "GLON", - "GLAT", - "ALT", - "AZ", - "DETX", - "DETY", - "THETA", - "PHI", - "GAMANESS", - "DIR_ERR", - "ENERGY_ERR", - "COREX", - "COREY", - "CORE_ERR", - "HMAX", - "HMAX_ERR", - ] - - if not self.raise_error_for_optional: - for i, c in enumerate(rename_from_optional): - if c not in events.colnames: - self.log.warning( - f"Optional column {c} is missing from the events table." - ) - else: - rename_from.append(rename_from_optional[i]) - rename_to.append(rename_to_optional[i]) + if self.optional_dl3_columns: + rename_from_optional = [ + "multiplicity", + "reco_glon", + "reco_glat", + "reco_alt", + "reco_az", + "reco_fov_lon", + "reco_fov_lat", + "reco_source_fov_offset", + "reco_source_fov_position_angle", + "gh_score", + "reco_dir_uncert", + "reco_energy_uncert", + "reco_core_x", + "reco_core_y", + "reco_core_uncert", + "reco_h_max", + "reco_h_max_uncert", + ] + rename_to_optional = [ + "MULTIP", + "GLON", + "GLAT", + "ALT", + "AZ", + "DETX", + "DETY", + "THETA", + "PHI", + "GAMANESS", + "DIR_ERR", + "ENERGY_ERR", + "COREX", + "COREY", + "CORE_ERR", + "HMAX", + "HMAX_ERR", + ] + + if not self.raise_error_for_optional: + for i, c in enumerate(rename_from_optional): + if c not in events.colnames: + self.log.warning( + f"Optional column {c} is missing from the events table." + ) + else: + rename_from.append(rename_from_optional[i]) + rename_to.append(rename_to_optional[i]) for c in rename_from: if c not in events.colnames: From 8e03bdbf58f84be6ea45120024ec6190a75102c8 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Wed, 12 Mar 2025 19:26:36 +0100 Subject: [PATCH 46/58] Add metadata header to the DL3 file --- .codespell-ignores | 1 + src/ctapipe/irf/dl3.py | 64 +++++++++++++++++++++++++++++--- src/ctapipe/irf/preprocessing.py | 2 +- src/ctapipe/tools/create_dl3.py | 1 + 4 files changed, 61 insertions(+), 7 deletions(-) diff --git a/.codespell-ignores b/.codespell-ignores index 446610f6aa9..c1bccebd1bb 100644 --- a/.codespell-ignores +++ b/.codespell-ignores @@ -3,3 +3,4 @@ usera nd studi referenc +livetime diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 6b2102b346d..6e964d248ad 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -3,7 +3,8 @@ from typing import Any, Dict, List, Tuple import astropy.units as u -from astropy.coordinates import EarthLocation, SkyCoord +import numpy as np +from astropy.coordinates import ICRS, EarthLocation, SkyCoord from astropy.io import fits from astropy.io.fits import Header from astropy.io.fits.hdu.base import ExtensionHDU @@ -38,6 +39,7 @@ class DL3_Format(Component): def __init__(self, **kwargs): super().__init__(**kwargs) + self._obs_id = None self._events = None self._pointing = None self._gti = None @@ -55,6 +57,18 @@ def __init__(self, **kwargs): def write_file(self, path): pass + @property + def obs_id(self) -> int: + return self._obs_id + + @obs_id.setter + def obs_id(self, obs_id: int): + if self._obs_id is not None: + self.log.warning( + "Obs id for DL3 file was already set, replacing current obs id" + ) + self._obs_id = obs_id + @property def events(self) -> QTable: return self._events @@ -204,14 +218,12 @@ def __init__(self, **kwargs): def write_file(self, path): self.file_creation_time = datetime.now(tz=datetime.UTC).isoformat() - events = self.transform_events_columns_for_gadf_format(self.events) - hdu_dl3 = fits.HDUList( [fits.PrimaryHDU(header=Header(self.get_hdu_header_base_format()))] ) hdu_dl3.append( fits.BinTableHDU( - data=events, + data=self.transform_events_columns_for_gadf_format(self.events), name="EVENTS", header=Header(self.get_hdu_header_events()), ) @@ -233,12 +245,14 @@ def get_hdu_header_base_format(self): } def get_hdu_header_base_time(self): - if self._gti is None: + if self.gti is None: raise ValueError("No available time information for the DL3 file") + if self.dead_time_fraction is None: + raise ValueError("No available dead time fraction for the DL3 file") start_time = None stop_time = None ontime = TimeDelta(0.0 * u.s) - for gti_interval in self._gti: + for gti_interval in self.gti: ontime += gti_interval[1] - gti_interval[0] start_time = ( gti_interval[0] @@ -262,6 +276,8 @@ def get_hdu_header_base_time(self): "TSTART": start_time, "TSTOP": stop_time, "ONTIME": ontime.to_value(u.s), + "LIVETIME": ontime.to_value(u.s) * self.dead_time_fraction, + "DEADC": self.dead_time_fraction, "TELAPSE": (stop_time - start_time).to_value(u.s), "DATE-OBS": start_time.fits, "DATE-BEG": start_time.fits, @@ -269,10 +285,46 @@ def get_hdu_header_base_time(self): "DATE-END": stop_time.fits, } + def get_hdu_header_observation_information(self): + if self.obs_id is None: + raise ValueError("Observation ID is missing.") + header = {"OBS_ID": self.obs_id} + if self.target_information is not None: + header["OBSERVER"] = self.target_information["observer"] + header["OBJECT"] = self.target_information["object_name"] + object_coordinate = self.target_information[ + "object_coordinate" + ].transform_to(ICRS()) + header["RA_OBJ"] = object_coordinate.ra.to_value(u.deg) + header["DEC_OBJ"] = object_coordinate.dec.to_value(u.deg) + + def get_hdu_header_subarray_information(self): + if self.telescope_information is None: + raise ValueError("Telescope information are missing.") + header = { + "ORIGIN": self.telescope_information["organisation"], + "TELESCOP": self.telescope_information["array"], + "INSTRUME": self.telescope_information["subarray"], + "TELLIST": self.telescope_information["telescope_list"], + "N_TELS": np.sum(self.telescope_information["telescope_list"]), + } + return header + + def get_hdu_header_software_information(self): + header = {} + if self.software_information is not None: + header["DST_VER"] = self.software_information["dst_version"] + header["ANA_VER"] = self.software_information["analysis_version"] + header["CAL_VER"] = self.software_information["calibration_version"] + return header + def get_hdu_header_events(self): header = self.get_hdu_header_base_format() header.update({"HDUCLAS1": "EVENTS"}) header.update(self.get_hdu_header_base_time()) + header.update(self.get_hdu_header_observation_information()) + header.update(self.get_hdu_header_subarray_information()) + header.update(self.get_hdu_header_software_information()) return header def get_hdu_header_gti(self): diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 034499adfbb..4e753b1318c 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -648,7 +648,7 @@ def get_observation_information(self) -> Dict[str, Any]: ) meta["target"] = { "observer": "UNKNOWN", - "object": "UNKNOWN", + "object_name": "UNKNOWN", "object_coordinate": ICRS(np.nan * u.deg, np.nan * u.deg), } diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index 2459262024d..f53410c4943 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -112,6 +112,7 @@ def start(self): self.chunk_size ) meta = self.event_loader.get_observation_information() + self.dl3_format.obs_id = meta["obs_id"] self.dl3_format.pointing = meta["pointing"] self.dl3_format.gti = meta["gti"] self.dl3_format.dead_time_fraction = meta["dead_time_fraction"] From df027d1f86979fa2024af14aac6c88f4bd6f774f Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Wed, 12 Mar 2025 19:43:23 +0100 Subject: [PATCH 47/58] Add GTI table in the DL3 file --- src/ctapipe/irf/dl3.py | 38 ++++++++++++++++++++++++++++++++------ 1 file changed, 32 insertions(+), 6 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 6e964d248ad..e0fa9b3516a 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -214,6 +214,7 @@ class DL3_GADF(DL3_Format): def __init__(self, **kwargs): super().__init__(**kwargs) self.file_creation_time = datetime.now(tz=datetime.UTC).isoformat() + self.reference_time = Time(datetime.fromisoformat("1970-01-01T00:00:00+00:00")) def write_file(self, path): self.file_creation_time = datetime.now(tz=datetime.UTC).isoformat() @@ -228,6 +229,13 @@ def write_file(self, path): header=Header(self.get_hdu_header_events()), ) ) + hdu_dl3.append( + fits.BinTableHDU( + data=self.create_gti_table(), + name="GTI", + header=Header(self.get_hdu_header_gti()), + ) + ) hdu_dl3.append(self.aeff) hdu_dl3.append(self.psf) hdu_dl3.append(self.edisp) @@ -265,11 +273,9 @@ def get_hdu_header_base_time(self): else max(stop_time, gti_interval[1]) ) - reference_time = Time(datetime.fromisoformat("1970-01-01T00:00:00+00:00")) - return { - "MJDREFI": int(reference_time.mjd), - "MJDREFF": reference_time.mjd % 1, + "MJDREFI": int(self.reference_time.mjd), + "MJDREFF": self.reference_time.mjd % 1, "TIMEUNIT": "s", "TIMEREF": "GEOCENTER", "TIMESYS": "UTC", @@ -285,11 +291,11 @@ def get_hdu_header_base_time(self): "DATE-END": stop_time.fits, } - def get_hdu_header_observation_information(self): + def get_hdu_header_observation_information(self, obs_id_only=False): if self.obs_id is None: raise ValueError("Observation ID is missing.") header = {"OBS_ID": self.obs_id} - if self.target_information is not None: + if self.target_information is not None and not obs_id_only: header["OBSERVER"] = self.target_information["observer"] header["OBJECT"] = self.target_information["object_name"] object_coordinate = self.target_information[ @@ -297,6 +303,7 @@ def get_hdu_header_observation_information(self): ].transform_to(ICRS()) header["RA_OBJ"] = object_coordinate.ra.to_value(u.deg) header["DEC_OBJ"] = object_coordinate.dec.to_value(u.deg) + return header def get_hdu_header_subarray_information(self): if self.telescope_information is None: @@ -331,6 +338,7 @@ def get_hdu_header_gti(self): header = self.get_hdu_header_base_format() header.update({"HDUCLAS1": "GTI"}) header.update(self.get_hdu_header_base_time()) + header.update(self.get_hdu_header_observation_information(obs_id_only=True)) return header def transform_events_columns_for_gadf_format(self, events): @@ -397,3 +405,21 @@ def transform_events_columns_for_gadf_format(self, events): renamed_events.rename_columns(rename_from, rename_to) renamed_events = renamed_events[rename_to] return renamed_events + + def create_gti_table(self) -> QTable: + table_structure = {"START": [], "STOP": []} + for gti_interval in self.gti: + table_structure["START"].append( + (gti_interval[0] - self.reference_time).to(u.s) + ) + table_structure["STOP"].append( + (gti_interval[1] - self.reference_time).to(u.s) + ) + + table = QTable(table_structure).sort("START") + for i in range(len(QTable) - 1): + if table_structure["STOP"][i] > table_structure["START"][i + 1]: + self.log.warning("Overlapping GTI intervals") + break + + return table From 32231555205bf1440987efba570d71ca042029fa Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Wed, 12 Mar 2025 21:16:36 +0100 Subject: [PATCH 48/58] Add pointing table and header to the DL3 file --- src/ctapipe/irf/dl3.py | 150 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 141 insertions(+), 9 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index e0fa9b3516a..b87403e3601 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -4,7 +4,7 @@ import astropy.units as u import numpy as np -from astropy.coordinates import ICRS, EarthLocation, SkyCoord +from astropy.coordinates import ICRS, AltAz, EarthLocation, SkyCoord from astropy.io import fits from astropy.io.fits import Header from astropy.io.fits.hdu.base import ExtensionHDU @@ -12,6 +12,7 @@ from astropy.time import Time, TimeDelta from ctapipe.compat import COPY_IF_NEEDED +from ctapipe.coordinates.tests.test_coordinates import location from ctapipe.core import Component from ctapipe.core.traits import Bool from ctapipe.version import version as ctapipe_version @@ -236,10 +237,29 @@ def write_file(self, path): header=Header(self.get_hdu_header_gti()), ) ) + hdu_dl3.append( + fits.BinTableHDU( + data=self.create_pointing_table(), + name="POINTING", + header=Header(self.get_hdu_header_pointing()), + ) + ) + + if self.aeff is None: + raise ValueError("Missing effective area IRF") hdu_dl3.append(self.aeff) + hdu_dl3[-1].header["OBS_ID"] = self.obs_id + if self.psf is None: + raise ValueError("Missing PSF IRF") hdu_dl3.append(self.psf) + hdu_dl3[-1].header["OBS_ID"] = self.obs_id + if self.edisp is None: + raise ValueError("Missing EDISP IRF") hdu_dl3.append(self.edisp) - hdu_dl3.append(self.bkg) + hdu_dl3[-1].header["OBS_ID"] = self.obs_id + if self.bkg is not None: + hdu_dl3.append(self.bkg) + hdu_dl3[-1].header["OBS_ID"] = self.obs_id hdu_dl3.writeto(path, checksum=True, overwrite=self.overwrite) @@ -291,7 +311,7 @@ def get_hdu_header_base_time(self): "DATE-END": stop_time.fits, } - def get_hdu_header_observation_information(self, obs_id_only=False): + def get_hdu_header_base_observation_information(self, obs_id_only=False): if self.obs_id is None: raise ValueError("Observation ID is missing.") header = {"OBS_ID": self.obs_id} @@ -305,7 +325,7 @@ def get_hdu_header_observation_information(self, obs_id_only=False): header["DEC_OBJ"] = object_coordinate.dec.to_value(u.deg) return header - def get_hdu_header_subarray_information(self): + def get_hdu_header_base_subarray_information(self): if self.telescope_information is None: raise ValueError("Telescope information are missing.") header = { @@ -317,7 +337,7 @@ def get_hdu_header_subarray_information(self): } return header - def get_hdu_header_software_information(self): + def get_hdu_header_base_software_information(self): header = {} if self.software_information is not None: header["DST_VER"] = self.software_information["dst_version"] @@ -325,20 +345,101 @@ def get_hdu_header_software_information(self): header["CAL_VER"] = self.software_information["calibration_version"] return header + def get_hdu_header_base_pointing(self): + if self.pointing is None: + raise ValueError("Pointing information are missing") + if self.location is None: + raise ValueError("Telescope location information are missing") + + gti_table = self.create_gti_table() + time_evaluation = [] + for i in range(len(gti_table)): + time_evaluation += list( + np.linspace(gti_table["START"][i], gti_table["STOP"][i], 100) + ) + time_evaluation = self.reference_time + TimeDelta(time_evaluation) + + pointing_table = self.create_pointing_table() + if self.pointing["pointing_mode"] == "TRACK": + obs_mode = "POINTING" + icrs_coordinate = SkyCoord( + ra=np.interp( + time_evaluation, + xp=pointing_table["TIME"], + yp=pointing_table["RA_PNT"], + ), + dec=np.interp( + time_evaluation, + xp=pointing_table["TIME"], + yp=pointing_table["DEC_PNT"], + ), + ) + altaz_coordinate = icrs_coordinate.transform_to( + AltAz(location=self.location, obstime=time_evaluation) + ) + elif self.pointing["pointing_mode"] == "DRIFT": + obs_mode = "DRIFT" + altaz_coordinate = AltAz( + alt=np.interp( + time_evaluation, + xp=pointing_table["TIME"], + yp=pointing_table["ALT_PNT"], + ), + az=np.interp( + time_evaluation, + xp=pointing_table["TIME"], + yp=pointing_table["AZ_PNT"], + ), + location=self.location, + obstime=time_evaluation, + ) + icrs_coordinate = altaz_coordinate.transform_to(ICRS()) + else: + raise ValueError("Unknown pointing mode") + + header = { + "RADECSYS": "ICRS", + "EQUINOX": 2000.0, + "OBS_MODE": obs_mode, + "RA_PNT": np.mean(icrs_coordinate.ra.to_value(u.deg)), + "DEC_PNT": np.mean(icrs_coordinate.dec.to_value(u.deg)), + "ALT_PNT": np.mean(altaz_coordinate.alt.to_value(u.deg)), + "AZ_PNT": np.mean(altaz_coordinate.az.to_value(u.deg)), + "GEOLON": self.location.lon.to_value(u.deg), + "GEOLAT": self.location.lat.to_value(u.deg), + "ALTITUDE": self.location.altitude.to_value(u.m), + "OBSGEO-X": self.location.x.to_value(u.m), + "OBSGEO-Y": self.location.y.to_value(u.m), + "OBSGEO-Z": self.location.z.to_value(u.m), + } + return header + def get_hdu_header_events(self): header = self.get_hdu_header_base_format() header.update({"HDUCLAS1": "EVENTS"}) header.update(self.get_hdu_header_base_time()) - header.update(self.get_hdu_header_observation_information()) - header.update(self.get_hdu_header_subarray_information()) - header.update(self.get_hdu_header_software_information()) + header.update(self.get_hdu_header_base_pointing()) + header.update(self.get_hdu_header_base_observation_information()) + header.update(self.get_hdu_header_base_subarray_information()) + header.update(self.get_hdu_header_base_software_information()) return header def get_hdu_header_gti(self): header = self.get_hdu_header_base_format() header.update({"HDUCLAS1": "GTI"}) header.update(self.get_hdu_header_base_time()) - header.update(self.get_hdu_header_observation_information(obs_id_only=True)) + header.update( + self.get_hdu_header_base_observation_information(obs_id_only=True) + ) + return header + + def get_hdu_header_pointing(self): + header = self.get_hdu_header_base_format() + header.update({"HDUCLAS1": "POINTING"}) + header.update(self.get_hdu_header_base_pointing()) + header.update( + self.get_hdu_header_base_observation_information(obs_id_only=True) + ) return header def transform_events_columns_for_gadf_format(self, events): @@ -423,3 +524,34 @@ def create_gti_table(self) -> QTable: break return table + + def create_pointing_table(self) -> QTable: + if self.pointing is None: + raise ValueError("Pointing information are missing") + if self.location is None: + raise ValueError("Telescope location information are missing") + + table_structure = { + "TIME": [], + "RA_PNT": [], + "DEC_PNT": [], + "ALT_PNT": [], + "AZ_PNT": [], + } + + for pointing in self.pointing["pointing_list"]: + time = pointing[0] + pointing_icrs = pointing[1].transform_to(ICRS()) + pointing_altaz = pointing[1].transform_to( + AltAz(location=location, obstime=time) + ) + table_structure = { + "TIME": (time - self.reference_time).to(u.s), + "RA_PNT": pointing_icrs.ra.to(u.deg), + "DEC_PNT": pointing_icrs.dec.to(u.deg), + "ALT_PNT": pointing_altaz.alt.to(u.deg), + "AZ_PNT": pointing_altaz.az.to(u.deg), + } + + table = QTable(table_structure).sort("TIME") + return table From e3c52457205c64df7dc09f71df945017a0636f7c Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Wed, 12 Mar 2025 21:17:20 +0100 Subject: [PATCH 49/58] Fix interpolation --- src/ctapipe/irf/dl3.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index b87403e3601..88c03bdcb78 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -366,12 +366,12 @@ def get_hdu_header_base_pointing(self): ra=np.interp( time_evaluation, xp=pointing_table["TIME"], - yp=pointing_table["RA_PNT"], + fp=pointing_table["RA_PNT"], ), dec=np.interp( time_evaluation, xp=pointing_table["TIME"], - yp=pointing_table["DEC_PNT"], + fp=pointing_table["DEC_PNT"], ), ) altaz_coordinate = icrs_coordinate.transform_to( @@ -383,12 +383,12 @@ def get_hdu_header_base_pointing(self): alt=np.interp( time_evaluation, xp=pointing_table["TIME"], - yp=pointing_table["ALT_PNT"], + fp=pointing_table["ALT_PNT"], ), az=np.interp( time_evaluation, xp=pointing_table["TIME"], - yp=pointing_table["AZ_PNT"], + fp=pointing_table["AZ_PNT"], ), location=self.location, obstime=time_evaluation, From ae70fd42044f0bd782d8f8cead8487f53538d418 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Thu, 13 Mar 2025 12:08:08 +0100 Subject: [PATCH 50/58] Add docstring to the dl3.py file --- src/ctapipe/irf/dl3.py | 196 +++++++++++++++++++++++++++++--- src/ctapipe/tools/create_dl3.py | 3 +- 2 files changed, 180 insertions(+), 19 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 88c03bdcb78..db74356b45a 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -19,6 +19,10 @@ class DL3_Format(Component): + """ + Base class for writing a DL3 file + """ + overwrite = Bool( default_value=False, help="If true, allow to overwrite already existing output file", @@ -33,16 +37,12 @@ class DL3_Format(Component): help="If true will raise error in the case optional column are missing", ).tag(config=False) - raise_error_for_missing_hdu = Bool( - default_value=True, - help="If true will raise error if HDU are missing from the final DL3 file", - ).tag(config=False) - def __init__(self, **kwargs): super().__init__(**kwargs) self._obs_id = None self._events = None self._pointing = None + self._pointing_mode = None self._gti = None self._aeff = None self._psf = None @@ -64,6 +64,12 @@ def obs_id(self) -> int: @obs_id.setter def obs_id(self, obs_id: int): + """ + Parameters + ---------- + obs_id : int + Observation ID + """ if self._obs_id is not None: self.log.warning( "Obs id for DL3 file was already set, replacing current obs id" @@ -76,6 +82,12 @@ def events(self) -> QTable: @events.setter def events(self, events: QTable): + """ + Parameters + ---------- + events : QTable + A table with a line for each event + """ if self._events is not None: self.log.warning( "Events table for DL3 file was already set, replacing current event table" @@ -88,18 +100,48 @@ def pointing(self) -> List[Tuple[Time, SkyCoord]]: @pointing.setter def pointing(self, pointing: List[Tuple[Time, SkyCoord]]): + """ + Parameters + ---------- + pointing : QTable + A list with for each entry containing the time at which the coordinate where evaluated and the associated coordinates + """ if self._pointing is not None: self.log.warning( "Pointing for DL3 file was already set, replacing current pointing" ) self._pointing = pointing + @property + def pointing_mode(self) -> str: + return self._pointing_mode + + @pointing_mode.setter + def pointing_mode(self, pointing_mode: str): + """ + Parameters + ---------- + pointing_mode : str + The name of the pointing mode used for the observation + """ + if self._pointing_mode is not None: + self.log.warning( + "Pointing for DL3 file was already set, replacing current pointing" + ) + self._pointing_mode = pointing_mode + @property def gti(self) -> List[Tuple[Time, Time]]: return self._gti @gti.setter def gti(self, gti: List[Tuple[Time, Time]]): + """ + Parameters + ---------- + gti : QTable + A list with for each entry containing the time at which the coordinate where evaluated and + """ if self._gti is not None: self.log.warning("GTI for DL3 file was already set, replacing current gti") self._gti = gti @@ -110,6 +152,12 @@ def aeff(self) -> ExtensionHDU: @aeff.setter def aeff(self, aeff: ExtensionHDU): + """ + Parameters + ---------- + aeff : ExtensionHDU + The effective area HDU read from the fits file containing IRFs + """ if self._aeff is not None: self.log.warning( "Effective area for DL3 file was already set, replacing current effective area" @@ -122,6 +170,12 @@ def psf(self) -> ExtensionHDU: @psf.setter def psf(self, psf: ExtensionHDU): + """ + Parameters + ---------- + psf : ExtensionHDU + The PSF HDU read from the fits file containing IRFs + """ if self._psf is not None: self.log.warning("PSF for DL3 file was already set, replacing current PSF") self._psf = psf @@ -132,6 +186,12 @@ def edisp(self) -> ExtensionHDU: @edisp.setter def edisp(self, edisp: ExtensionHDU): + """ + Parameters + ---------- + edisp : ExtensionHDU + The EDISP HDU read from the fits file containing IRFs + """ if self._edisp is not None: self.log.warning( "EDISP for DL3 file was already set, replacing current EDISP" @@ -144,6 +204,12 @@ def bkg(self) -> ExtensionHDU: @bkg.setter def bkg(self, bkg: ExtensionHDU): + """ + Parameters + ---------- + bkg : ExtensionHDU + The background HDU read from the fits file containing IRFs + """ if self._bkg is not None: self.log.warning( "Background for DL3 file was already set, replacing current background" @@ -156,6 +222,12 @@ def location(self) -> EarthLocation: @location.setter def location(self, location: EarthLocation): + """ + Parameters + ---------- + location : EarthLocation + The location of the telescope + """ if self._location is not None: self.log.warning( "Telescope location for DL3 file was already set, replacing current location" @@ -168,6 +240,12 @@ def dead_time_fraction(self) -> float: @dead_time_fraction.setter def dead_time_fraction(self, dead_time_fraction: float): + """ + Parameters + ---------- + dead_time_fraction : float + The location of the telescope + """ if self.dead_time_fraction is not None: self.log.warning( "Dead time fraction for DL3 file was already set, replacing current dead time fraction" @@ -180,6 +258,12 @@ def telescope_information(self) -> Dict[str, Any]: @telescope_information.setter def telescope_information(self, telescope_information: Dict[str, Any]): + """ + Parameters + ---------- + telescope_information : dict[str, any] + A dictionary containing general information about telescope with as key : organisation, array, subarray, telescope_list + """ if self._telescope_information is not None: self.log.warning( "Telescope information for DL3 file was already set, replacing current information" @@ -192,6 +276,12 @@ def target_information(self) -> Dict[str, Any]: @target_information.setter def target_information(self, target_information: Dict[str, Any]): + """ + Parameters + ---------- + target_information : dict[str, any] + A dictionary containing general information about the targeted source with as key : observer, object_name, object_coordinate + """ if self._target_information is not None: self.log.warning( "Target information for DL3 file was already set, replacing current target information" @@ -204,6 +294,12 @@ def software_information(self) -> Dict[str, Any]: @software_information.setter def software_information(self, software_information: Dict[str, Any]): + """ + Parameters + ---------- + software_information : dict[str, any] + A dictionary containing general information about the software used to produce the file with as key : analysis_version, calibration_version, dst_version + """ if self._software_information is not None: self.log.warning( "Software information for DL3 file was already set, replacing current software information" @@ -212,12 +308,25 @@ def software_information(self, software_information: Dict[str, Any]): class DL3_GADF(DL3_Format): + """ + Class to write DL3 in GADF format, subclass of DL3_Format + """ + def __init__(self, **kwargs): super().__init__(**kwargs) self.file_creation_time = datetime.now(tz=datetime.UTC).isoformat() self.reference_time = Time(datetime.fromisoformat("1970-01-01T00:00:00+00:00")) def write_file(self, path): + """ + This function will write the new DL3 file + All the content associated with the file should have been specified previously, otherwise error will be raised + + Parameters + ---------- + path : str + The full path and filename of the new file to write + """ self.file_creation_time = datetime.now(tz=datetime.UTC).isoformat() hdu_dl3 = fits.HDUList( @@ -263,7 +372,10 @@ def write_file(self, path): hdu_dl3.writeto(path, checksum=True, overwrite=self.overwrite) - def get_hdu_header_base_format(self): + def get_hdu_header_base_format(self) -> Dict[str, Any]: + """ + Return the base information that should be included in all HDU of the final fits file + """ return { "HDUCLASS": "GADF", "HDUVERS": "v0.3", @@ -272,7 +384,10 @@ def get_hdu_header_base_format(self): "CREATED": self.file_creation_time, } - def get_hdu_header_base_time(self): + def get_hdu_header_base_time(self) -> Dict[str, Any]: + """ + Return the information about time parameters used in several HDU + """ if self.gti is None: raise ValueError("No available time information for the DL3 file") if self.dead_time_fraction is None: @@ -311,7 +426,17 @@ def get_hdu_header_base_time(self): "DATE-END": stop_time.fits, } - def get_hdu_header_base_observation_information(self, obs_id_only=False): + def get_hdu_header_base_observation_information( + self, obs_id_only: bool = False + ) -> Dict[str, Any]: + """ + Return generic information on the observation setting (id, target, ...) + + Parameters + ---------- + obs_id_only : bool + If true, will return a dict with as only information the obs_id + """ if self.obs_id is None: raise ValueError("Observation ID is missing.") header = {"OBS_ID": self.obs_id} @@ -325,7 +450,10 @@ def get_hdu_header_base_observation_information(self, obs_id_only=False): header["DEC_OBJ"] = object_coordinate.dec.to_value(u.deg) return header - def get_hdu_header_base_subarray_information(self): + def get_hdu_header_base_subarray_information(self) -> Dict[str, Any]: + """ + Return generic information on the array used for observations + """ if self.telescope_information is None: raise ValueError("Telescope information are missing.") header = { @@ -337,7 +465,10 @@ def get_hdu_header_base_subarray_information(self): } return header - def get_hdu_header_base_software_information(self): + def get_hdu_header_base_software_information(self) -> Dict[str, Any]: + """ + Return information about the software versions used to process the observation + """ header = {} if self.software_information is not None: header["DST_VER"] = self.software_information["dst_version"] @@ -345,9 +476,14 @@ def get_hdu_header_base_software_information(self): header["CAL_VER"] = self.software_information["calibration_version"] return header - def get_hdu_header_base_pointing(self): + def get_hdu_header_base_pointing(self) -> Dict[str, Any]: + """ + Return information on the pointing during the observation + """ if self.pointing is None: raise ValueError("Pointing information are missing") + if self.pointing_mode is None: + raise ValueError("Pointing mode is missing") if self.location is None: raise ValueError("Telescope location information are missing") @@ -360,7 +496,7 @@ def get_hdu_header_base_pointing(self): time_evaluation = self.reference_time + TimeDelta(time_evaluation) pointing_table = self.create_pointing_table() - if self.pointing["pointing_mode"] == "TRACK": + if self.pointing_mode == "TRACK": obs_mode = "POINTING" icrs_coordinate = SkyCoord( ra=np.interp( @@ -377,7 +513,7 @@ def get_hdu_header_base_pointing(self): altaz_coordinate = icrs_coordinate.transform_to( AltAz(location=self.location, obstime=time_evaluation) ) - elif self.pointing["pointing_mode"] == "DRIFT": + elif self.pointing_mode == "DRIFT": obs_mode = "DRIFT" altaz_coordinate = AltAz( alt=np.interp( @@ -414,7 +550,10 @@ def get_hdu_header_base_pointing(self): } return header - def get_hdu_header_events(self): + def get_hdu_header_events(self) -> Dict[str, Any]: + """ + The output dictionary contain all the necessary information that should be added to the header of the events HDU + """ header = self.get_hdu_header_base_format() header.update({"HDUCLAS1": "EVENTS"}) header.update(self.get_hdu_header_base_time()) @@ -424,7 +563,10 @@ def get_hdu_header_events(self): header.update(self.get_hdu_header_base_software_information()) return header - def get_hdu_header_gti(self): + def get_hdu_header_gti(self) -> Dict[str, Any]: + """ + The output dictionary contain all the necessary information that should be added to the header of the GTI HDU + """ header = self.get_hdu_header_base_format() header.update({"HDUCLAS1": "GTI"}) header.update(self.get_hdu_header_base_time()) @@ -433,7 +575,10 @@ def get_hdu_header_gti(self): ) return header - def get_hdu_header_pointing(self): + def get_hdu_header_pointing(self) -> Dict[str, Any]: + """ + The output dictionary contain all the necessary information that should be added to the header of the pointing HDU + """ header = self.get_hdu_header_base_format() header.update({"HDUCLAS1": "POINTING"}) header.update(self.get_hdu_header_base_pointing()) @@ -442,7 +587,16 @@ def get_hdu_header_pointing(self): ) return header - def transform_events_columns_for_gadf_format(self, events): + def transform_events_columns_for_gadf_format(self, events: QTable) -> QTable: + """ + Return an event table containing only the columns that should be added to the EVENTS HDU + It also rename all the columns to match the name expected in the GADF format + + Parameters + ---------- + events : QTable + The base events table to process + """ rename_from = ["event_id", "time", "reco_ra", "reco_dec", "reco_energy"] rename_to = ["EVENT_ID", "TIME", "RA", "DEC", "ENERGY"] @@ -508,6 +662,9 @@ def transform_events_columns_for_gadf_format(self, events): return renamed_events def create_gti_table(self) -> QTable: + """ + Build a table that contains GTI information with the GADF names and format, to be concerted directly as a TableHDU + """ table_structure = {"START": [], "STOP": []} for gti_interval in self.gti: table_structure["START"].append( @@ -526,6 +683,9 @@ def create_gti_table(self) -> QTable: return table def create_pointing_table(self) -> QTable: + """ + Build a table that contains pointing information with the GADF names and format, to be concerted directly as a TableHDU + """ if self.pointing is None: raise ValueError("Pointing information are missing") if self.location is None: @@ -539,7 +699,7 @@ def create_pointing_table(self) -> QTable: "AZ_PNT": [], } - for pointing in self.pointing["pointing_list"]: + for pointing in self.pointing: time = pointing[0] pointing_icrs = pointing[1].transform_to(ICRS()) pointing_altaz = pointing[1].transform_to( diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index f53410c4943..ff1e7d3b639 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -113,7 +113,8 @@ def start(self): ) meta = self.event_loader.get_observation_information() self.dl3_format.obs_id = meta["obs_id"] - self.dl3_format.pointing = meta["pointing"] + self.dl3_format.pointing = meta["pointing"]["pointing_list"] + self.dl3_format.pointing_mode = meta["pointing"]["pointing_mode"] self.dl3_format.gti = meta["gti"] self.dl3_format.dead_time_fraction = meta["dead_time_fraction"] From 917240fb9ad92eb9e37e24c86da1cda4b6e89fcc Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Mon, 17 Mar 2025 15:47:14 +0100 Subject: [PATCH 51/58] Fix errors encountered in automatic tests --- src/ctapipe/conftest.py | 9 +- src/ctapipe/io/hdf5eventsource.py | 24 +- src/ctapipe/irf/preprocessing.py | 2 +- src/ctapipe/tools/optimize_event_selection.py | 6 +- .../tools/tests/test_calculate_pixel_stats.py | 129 -------- src/ctapipe/tools/tests/test_compute_irf.py | 284 ------------------ src/ctapipe/tools/tests/test_merge.py | 6 +- .../tests/test_optimize_event_selection.py | 116 ------- src/ctapipe/tools/tests/test_process.py | 12 +- 9 files changed, 35 insertions(+), 553 deletions(-) delete mode 100644 src/ctapipe/tools/tests/test_calculate_pixel_stats.py delete mode 100644 src/ctapipe/tools/tests/test_compute_irf.py delete mode 100644 src/ctapipe/tools/tests/test_optimize_event_selection.py diff --git a/src/ctapipe/conftest.py b/src/ctapipe/conftest.py index 7af605dc5e6..b681004f1f2 100644 --- a/src/ctapipe/conftest.py +++ b/src/ctapipe/conftest.py @@ -819,7 +819,14 @@ def irf_events_table(): keep_columns, _, _ = epp.get_columns_keep_rename_scheme(None, True) tab = epp.make_empty_table(keep_columns) - ids, bulk, unitless = tab.colnames[:2], tab.colnames[2:-2], tab.colnames[-2:] + ids, bulk, unitless = [], [], [] + for c in tab.columns: + if "id" in c: + ids.append(c) + elif tab[c].unit == u.dimensionless_unscaled: + unitless.append(c) + else: + bulk.append(c) id_tab = QTable( data=np.zeros((N, len(ids)), dtype=np.uint64), diff --git a/src/ctapipe/io/hdf5eventsource.py b/src/ctapipe/io/hdf5eventsource.py index 779e1e6d6ba..8e8cb3ad6b8 100644 --- a/src/ctapipe/io/hdf5eventsource.py +++ b/src/ctapipe/io/hdf5eventsource.py @@ -340,7 +340,7 @@ def has_simulated_dl1(self): True for files with telescope-wise event information in the simulation group """ if self.is_simulation: - if "telescope" in self.file_.root.simulation.events: + if "telescope" in self.file_.root.simulation.event: return True return False @@ -384,7 +384,7 @@ def simulation_config(self) -> dict[int, SimulationConfigContainer]: return self._simulation_configs def __len__(self): - n_events = len(self.file_.root.dl1.events.subarray.trigger) + n_events = len(self.file_.root.dl1.event.subarray.trigger) if self.max_events is not None: return min(n_events, self.max_events) return n_events @@ -434,7 +434,7 @@ def _parse_sb_and_ob_configs(self): return scheduling_blocks, observation_blocks def _is_hillas_in_camera_frame(self): - parameters_group = self.file_.root.dl1.events.telescope.parameters + parameters_group = self.file_.root.dl1.event.telescope.parameters telescope_tables = parameters_group._v_children.values() # in case of no parameters, it doesn't matter, we just return False @@ -456,7 +456,7 @@ def _generator(self): table.name: self.reader.read( f"/r1/event/telescope/{table.name}", R1CameraContainer ) - for table in self.file_.root.r1.events.telescope + for table in self.file_.root.r1.event.telescope } if DataLevel.DL1_IMAGES in self.datalevels: @@ -472,14 +472,14 @@ def _generator(self): DL1CameraContainer, ignore_columns=ignore_columns, ) - for table in self.file_.root.dl1.events.telescope.images + for table in self.file_.root.dl1.event.telescope.images } if self.has_simulated_dl1: simulated_image_iterators = { - table.name: self.file_.root.simulation.events.telescope.images[ + table.name: self.file_.root.simulation.event.telescope.images[ table.name ].iterrows() - for table in self.file_.root.simulation.events.telescope.images + for table in self.file_.root.simulation.event.telescope.images } if DataLevel.DL1_PARAMETERS in self.datalevels: @@ -516,7 +516,7 @@ def _generator(self): "peak_time", ], ) - for table in self.file_.root.dl1.events.telescope.parameters + for table in self.file_.root.dl1.event.telescope.parameters } if self.has_simulated_dl1: simulated_param_readers = { @@ -537,7 +537,7 @@ def _generator(self): "true_intensity", ], ) - for table in self.file_.root.dl1.events.telescope.parameters + for table in self.file_.root.dl1.event.telescope.parameters } if self.has_muon_parameters: @@ -550,7 +550,7 @@ def _generator(self): MuonEfficiencyContainer, ], ) - for table in self.file_.root.dl1.events.telescope.muon + for table in self.file_.root.dl1.event.telescope.muon } dl2_readers = {} @@ -602,14 +602,14 @@ def _generator(self): SimulatedShowerContainer, prefixes="true", ) - if "impact" in self.file_.root.simulation.events.telescope: + if "impact" in self.file_.root.simulation.event.telescope: true_impact_readers = { table.name: self.reader.read( f"/simulation/event/telescope/impact/{table.name}", containers=TelescopeImpactParameterContainer, prefixes=["true_impact"], ) - for table in self.file_.root.simulation.events.telescope.impact + for table in self.file_.root.simulation.event.telescope.impact } # Setup iterators for the array events diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 4e753b1318c..b92b6289102 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -239,7 +239,7 @@ def normalise_column_names(self, events: Table) -> QTable: return renamed_events def keep_necessary_columns_only(self, events: Table) -> QTable: - keep_columns, _, _ = self.get_columns_keep_rename_scheme(events) + keep_columns, _, _ = self.get_columns_keep_rename_scheme(events, True) for c in keep_columns: if c not in events.colnames: raise ValueError( diff --git a/src/ctapipe/tools/optimize_event_selection.py b/src/ctapipe/tools/optimize_event_selection.py index 445c9cf3ac9..71cd5e08476 100644 --- a/src/ctapipe/tools/optimize_event_selection.py +++ b/src/ctapipe/tools/optimize_event_selection.py @@ -7,7 +7,7 @@ from ..core import Provenance, Tool, traits from ..core.traits import AstroQuantity, Integer, classes_with_traits -from ..irf import EventLoader, Spectra +from ..irf import EventLoader, EventPreprocessor, Spectra __all__ = ["EventSelectionOptimizer"] @@ -109,6 +109,10 @@ def setup(self): """ Initialize components from config. """ + + # Force the preprocessing for IRF + EventPreprocessor.irf_pre_processing = True + self.optimizer = CutOptimizerBase.from_name( self.optimization_algorithm, parent=self ) diff --git a/src/ctapipe/tools/tests/test_calculate_pixel_stats.py b/src/ctapipe/tools/tests/test_calculate_pixel_stats.py deleted file mode 100644 index 9e2e9674664..00000000000 --- a/src/ctapipe/tools/tests/test_calculate_pixel_stats.py +++ /dev/null @@ -1,129 +0,0 @@ -#!/usr/bin/env python3 -""" -Test ctapipe-calculate-pixel-statistics tool -""" - -import pytest -from traitlets.config.loader import Config - -from ctapipe.core import run_tool -from ctapipe.core.tool import ToolConfigurationError -from ctapipe.instrument import SubarrayDescription -from ctapipe.io import read_table -from ctapipe.tools.calculate_pixel_stats import PixelStatisticsCalculatorTool - - -def test_calculate_pixel_stats_tool(tmp_path, dl1_image_file): - """check statistics calculation from pixel-wise image data files""" - - # Create a configuration suitable for the test - tel_id = 3 - config = Config( - { - "PixelStatisticsCalculatorTool": { - "allowed_tels": [tel_id], - "input_column_name": "image", - "output_table_name": "statistics", - }, - "PixelStatisticsCalculator": { - "stats_aggregator_type": [ - ("id", tel_id, "PlainAggregator"), - ], - }, - "PlainAggregator": { - "chunk_size": 1, - }, - } - ) - # Set the output file path - monitoring_file = tmp_path / "monitoring.dl1.h5" - # Run the tool with the configuration and the input file - run_tool( - PixelStatisticsCalculatorTool(config=config), - argv=[ - f"--input_url={dl1_image_file}", - f"--output_path={monitoring_file}", - "--overwrite", - ], - cwd=tmp_path, - raises=True, - ) - # Check that the output file has been created - assert monitoring_file.exists() - # Check that the output file is not empty - assert ( - read_table( - monitoring_file, - path=f"/dl1/monitoring/telescope/statistics/tel_{tel_id:03d}", - )["mean"] - is not None - ) - # Read subarray description from the created monitoring file - subarray = SubarrayDescription.from_hdf(monitoring_file) - # Check for the selected telescope - assert subarray.tel_ids[0] == tel_id - - -def test_tool_config_error(tmp_path, dl1_image_file): - """check tool configuration error""" - - # Run the tool with the configuration and the input file - config = Config( - { - "PixelStatisticsCalculatorTool": { - "allowed_tels": [3], - "input_column_name": "image_charges", - "output_table_name": "statistics", - } - } - ) - # Set the output file path - monitoring_file = tmp_path / "monitoring.dl1.h5" - # Check if ToolConfigurationError is raised - # when the column name of the pixel-wise image data is not correct - with pytest.raises( - ToolConfigurationError, match="Column 'image_charges' not found" - ): - run_tool( - PixelStatisticsCalculatorTool(config=config), - argv=[ - f"--input_url={dl1_image_file}", - f"--output_path={monitoring_file}", - "--StatisticsAggregator.chunk_size=1", - "--overwrite", - ], - cwd=tmp_path, - raises=True, - ) - # Check if ToolConfigurationError is raised - # when the input and output files are the same - with pytest.raises( - ToolConfigurationError, match="Input and output files are same." - ): - run_tool( - PixelStatisticsCalculatorTool(), - argv=[ - f"--input_url={dl1_image_file}", - f"--output_path={dl1_image_file}", - "--overwrite", - ], - cwd=tmp_path, - raises=True, - ) - # Check if ToolConfigurationError is raised - # when the chunk size is larger than the number of events in the input file - with pytest.raises( - ToolConfigurationError, match="Change --StatisticsAggregator.chunk_size" - ): - run_tool( - PixelStatisticsCalculatorTool(), - argv=[ - f"--input_url={dl1_image_file}", - f"--output_path={monitoring_file}", - "--PixelStatisticsCalculatorTool.allowed_tels=3", - "--StatisticsAggregator.chunk_size=2500", - "--overwrite", - ], - cwd=tmp_path, - raises=True, - ) diff --git a/src/ctapipe/tools/tests/test_compute_irf.py b/src/ctapipe/tools/tests/test_compute_irf.py deleted file mode 100644 index acd35a50aaf..00000000000 --- a/src/ctapipe/tools/tests/test_compute_irf.py +++ /dev/null @@ -1,284 +0,0 @@ -import json -import logging - -import pytest -from astropy.io import fits - -from ctapipe.core import ToolConfigurationError, run_tool - -pytest.importorskip("pyirf") - - -@pytest.fixture(scope="module") -def dummy_cuts_file( - gamma_diffuse_full_reco_file, - proton_full_reco_file, - event_loader_config_path, - irf_tmp_path, -): - from ctapipe.tools.optimize_event_selection import EventSelectionOptimizer - - output_path = irf_tmp_path / "dummy_cuts.fits" - run_tool( - EventSelectionOptimizer(), - argv=[ - f"--gamma-file={gamma_diffuse_full_reco_file}", - f"--proton-file={proton_full_reco_file}", - # Use diffuse gammas weighted to electron spectrum as stand-in - f"--electron-file={gamma_diffuse_full_reco_file}", - f"--output={output_path}", - f"--config={event_loader_config_path}", - ], - ) - return output_path - - -@pytest.mark.parametrize("include_background", (False, True)) -@pytest.mark.parametrize("spatial_selection_applied", (True, False)) -def test_irf_tool( - gamma_diffuse_full_reco_file, - proton_full_reco_file, - event_loader_config_path, - dummy_cuts_file, - tmp_path, - include_background, - spatial_selection_applied, -): - from ctapipe.tools.compute_irf import IrfTool - - output_path = tmp_path / "irf.fits.gz" - output_benchmarks_path = tmp_path / "benchmarks.fits.gz" - - argv = [ - f"--gamma-file={gamma_diffuse_full_reco_file}", - f"--cuts={dummy_cuts_file}", - f"--output={output_path}", - f"--config={event_loader_config_path}", - ] - if spatial_selection_applied: - argv.append("--spatial-selection-applied") - - if include_background: - argv.append(f"--proton-file={proton_full_reco_file}") - # Use diffuse gammas weighted to electron spectrum as stand-in - argv.append(f"--electron-file={gamma_diffuse_full_reco_file}") - else: - argv.append("--no-do-background") - - ret = run_tool(IrfTool(), argv=argv) - assert ret == 0 - - assert output_path.exists() - assert not output_benchmarks_path.exists() - # Readability by gammapy is tested by pyirf tests, so not repeated here - with fits.open(output_path) as hdul: - assert isinstance(hdul["ENERGY DISPERSION"], fits.BinTableHDU) - assert isinstance(hdul["EFFECTIVE AREA"], fits.BinTableHDU) - assert isinstance(hdul["PSF"], fits.BinTableHDU) - if spatial_selection_applied: - assert isinstance(hdul["RAD_MAX"], fits.BinTableHDU) - - if include_background: - assert isinstance(hdul["BACKGROUND"], fits.BinTableHDU) - - output_path.unlink() # Delete output file - - # Include benchmarks - argv.append(f"--benchmark-output={output_benchmarks_path}") - ret = run_tool(IrfTool(), argv=argv) - assert ret == 0 - - assert output_path.exists() - assert output_benchmarks_path.exists() - with fits.open(output_benchmarks_path) as hdul: - assert isinstance(hdul["ENERGY BIAS RESOLUTION"], fits.BinTableHDU) - assert isinstance(hdul["ANGULAR RESOLUTION"], fits.BinTableHDU) - if include_background: - assert isinstance(hdul["SENSITIVITY"], fits.BinTableHDU) - - -def test_irf_tool_no_electrons( - gamma_diffuse_full_reco_file, - proton_full_reco_file, - event_loader_config_path, - dummy_cuts_file, - tmp_path, -): - from ctapipe.tools.compute_irf import IrfTool - - output_path = tmp_path / "irf.fits.gz" - output_benchmarks_path = tmp_path / "benchmarks.fits.gz" - logpath = tmp_path / "test_irf_tool_no_electrons.log" - logger = logging.getLogger("ctapipe.tools.compute_irf") - logger.addHandler(logging.FileHandler(logpath)) - - ret = run_tool( - IrfTool(), - argv=[ - f"--gamma-file={gamma_diffuse_full_reco_file}", - f"--proton-file={proton_full_reco_file}", - f"--cuts={dummy_cuts_file}", - f"--output={output_path}", - f"--benchmark-output={output_benchmarks_path}", - f"--config={event_loader_config_path}", - f"--log-file={logpath}", - ], - ) - assert ret == 0 - assert output_path.exists() - assert output_benchmarks_path.exists() - assert "Estimating background without electron file." in logpath.read_text() - - -def test_irf_tool_only_gammas( - gamma_diffuse_full_reco_file, event_loader_config_path, dummy_cuts_file, tmp_path -): - from ctapipe.tools.compute_irf import IrfTool - - output_path = tmp_path / "irf.fits.gz" - output_benchmarks_path = tmp_path / "benchmarks.fits.gz" - - with pytest.raises( - ValueError, - match="At least a proton file required when specifying `do_background`.", - ): - run_tool( - IrfTool(), - argv=[ - f"--gamma-file={gamma_diffuse_full_reco_file}", - f"--cuts={dummy_cuts_file}", - f"--output={output_path}", - f"--benchmark-output={output_benchmarks_path}", - f"--config={event_loader_config_path}", - ], - raises=True, - ) - - ret = run_tool( - IrfTool(), - argv=[ - f"--gamma-file={gamma_diffuse_full_reco_file}", - f"--cuts={dummy_cuts_file}", - f"--output={output_path}", - f"--benchmark-output={output_benchmarks_path}", - f"--config={event_loader_config_path}", - "--no-do-background", - ], - ) - assert ret == 0 - assert output_path.exists() - assert output_benchmarks_path.exists() - - -# TODO: Add test using point-like gammas - - -def test_point_like_irf_no_theta_cut( - gamma_diffuse_full_reco_file, - proton_full_reco_file, - event_loader_config_path, - dummy_cuts_file, - tmp_path, -): - from ctapipe.irf import OptimizationResult - from ctapipe.tools.compute_irf import IrfTool - - gh_cuts_path = tmp_path / "gh_cuts.fits" - cuts = OptimizationResult.read(dummy_cuts_file) - cuts.spatial_selection_table = None - cuts.write(gh_cuts_path) - assert gh_cuts_path.exists() - - output_path = tmp_path / "irf.fits.gz" - output_benchmarks_path = tmp_path / "benchmarks.fits.gz" - - with pytest.raises( - ToolConfigurationError, - match=rf"{gh_cuts_path} does not contain any direction cut", - ): - run_tool( - IrfTool(), - argv=[ - f"--gamma-file={gamma_diffuse_full_reco_file}", - f"--proton-file={proton_full_reco_file}", - # Use diffuse gammas weighted to electron spectrum as stand-in - f"--electron-file={gamma_diffuse_full_reco_file}", - f"--cuts={gh_cuts_path}", - f"--output={output_path}", - f"--benchmark-output={output_benchmarks_path}", - f"--config={event_loader_config_path}", - "--spatial-selection-applied", - ], - raises=True, - ) - - -def test_irf_tool_wrong_cuts( - gamma_diffuse_full_reco_file, proton_full_reco_file, dummy_cuts_file, tmp_path -): - from ctapipe.tools.compute_irf import IrfTool - - output_path = tmp_path / "irf.fits.gz" - output_benchmarks_path = tmp_path / "benchmarks.fits.gz" - - with pytest.raises(RuntimeError): - run_tool( - IrfTool(), - argv=[ - f"--gamma-file={gamma_diffuse_full_reco_file}", - f"--proton-file={proton_full_reco_file}", - # Use diffuse gammas weighted to electron spectrum as stand-in - f"--electron-file={gamma_diffuse_full_reco_file}", - f"--cuts={dummy_cuts_file}", - f"--output={output_path}", - f"--benchmark-output={output_benchmarks_path}", - ], - raises=True, - ) - - config_path = tmp_path / "config.json" - with config_path.open("w") as f: - json.dump( - { - "EventPreprocessor": { - "energy_reconstructor": "ExtraTreesRegressor", - "geometry_reconstructor": "HillasReconstructor", - "gammaness_classifier": "ExtraTreesClassifier", - "EventQualityQuery": { - "quality_criteria": [ - # No criteria for minimum event multiplicity - ("valid classifier", "ExtraTreesClassifier_is_valid"), - ("valid geom reco", "HillasReconstructor_is_valid"), - ("valid energy reco", "ExtraTreesRegressor_is_valid"), - ], - }, - } - }, - f, - ) - - logpath = tmp_path / "test_irf_tool_wrong_cuts.log" - logger = logging.getLogger("ctapipe.tools.compute_irf") - logger.addHandler(logging.FileHandler(logpath)) - - ret = run_tool( - IrfTool(), - argv=[ - f"--gamma-file={gamma_diffuse_full_reco_file}", - f"--proton-file={proton_full_reco_file}", - # Use diffuse gammas weighted to electron spectrum as stand-in - f"--electron-file={gamma_diffuse_full_reco_file}", - f"--cuts={dummy_cuts_file}", - f"--output={output_path}", - f"--benchmark-output={output_benchmarks_path}", - f"--config={config_path}", - f"--log-file={logpath}", - ], - ) - assert ret == 0 - assert output_path.exists() - assert output_benchmarks_path.exists() - assert ( - "Quality criteria are different from quality criteria " - "used for calculating g/h / theta cuts." in logpath.read_text() - ) diff --git a/src/ctapipe/tools/tests/test_merge.py b/src/ctapipe/tools/tests/test_merge.py index 2875ddb88ed..0879d27f1df 100644 --- a/src/ctapipe/tools/tests/test_merge.py +++ b/src/ctapipe/tools/tests/test_merge.py @@ -103,9 +103,9 @@ def test_skip_images(tmp_path, dl1_file, dl1_proton_file): ) with tables.open_file(output, "r") as f: - assert "images" not in f.root.dl1.events.telescope - assert "images" in f.root.simulation.events.telescope - assert "parameters" in f.root.dl1.events.telescope + assert "images" not in f.root.dl1.event.telescope + assert "images" in f.root.simulation.event.telescope + assert "parameters" in f.root.dl1.event.telescope t = read_table(output, "/simulation/event/telescope/images/tel_001") assert "true_image" not in t.colnames diff --git a/src/ctapipe/tools/tests/test_optimize_event_selection.py b/src/ctapipe/tools/tests/test_optimize_event_selection.py deleted file mode 100644 index cb787653a0a..00000000000 --- a/src/ctapipe/tools/tests/test_optimize_event_selection.py +++ /dev/null @@ -1,116 +0,0 @@ -import logging - -import astropy.units as u -import pytest -from astropy.table import QTable - -from ctapipe.core import QualityQuery, run_tool - -pytest.importorskip("pyirf") - - -def test_cuts_optimization( - gamma_diffuse_full_reco_file, - proton_full_reco_file, - event_loader_config_path, - tmp_path, -): - from ctapipe.irf import ( - OptimizationResult, - ResultValidRange, - ) - from ctapipe.tools.optimize_event_selection import EventSelectionOptimizer - - output_path = tmp_path / "cuts.fits" - - argv = [ - f"--gamma-file={gamma_diffuse_full_reco_file}", - f"--proton-file={proton_full_reco_file}", - # Use diffuse gammas weighted to electron spectrum as stand-in - f"--electron-file={gamma_diffuse_full_reco_file}", - f"--output={output_path}", - f"--config={event_loader_config_path}", - ] - ret = run_tool(EventSelectionOptimizer(), argv=argv) - assert ret == 0 - - result = OptimizationResult.read(output_path) - assert isinstance(result, OptimizationResult) - assert isinstance(result.quality_query, QualityQuery) - assert isinstance(result.valid_energy, ResultValidRange) - assert isinstance(result.valid_offset, ResultValidRange) - assert isinstance(result.gh_cuts, QTable) - assert result.clf_prefix == "ExtraTreesClassifier" - assert "cut" in result.gh_cuts.colnames - assert isinstance(result.spatial_selection_table, QTable) - assert "cut" in result.spatial_selection_table.colnames - - for c in ["low", "center", "high"]: - assert c in result.gh_cuts.colnames - assert result.gh_cuts[c].unit == u.TeV - assert c in result.spatial_selection_table.colnames - assert result.spatial_selection_table[c].unit == u.TeV - - -def test_cuts_opt_no_electrons( - gamma_diffuse_full_reco_file, - proton_full_reco_file, - event_loader_config_path, - tmp_path, -): - from ctapipe.tools.optimize_event_selection import EventSelectionOptimizer - - output_path = tmp_path / "cuts.fits" - logpath = tmp_path / "test_cuts_opt_no_electrons.log" - logger = logging.getLogger("ctapipe.tools.optimize_event_selection") - logger.addHandler(logging.FileHandler(logpath)) - - ret = run_tool( - EventSelectionOptimizer(), - argv=[ - f"--gamma-file={gamma_diffuse_full_reco_file}", - f"--proton-file={proton_full_reco_file}", - f"--output={output_path}", - f"--config={event_loader_config_path}", - f"--log-file={logpath}", - ], - ) - assert ret == 0 - assert "Optimizing cuts without electron file." in logpath.read_text() - - -def test_cuts_opt_only_gammas( - gamma_diffuse_full_reco_file, event_loader_config_path, tmp_path -): - from ctapipe.tools.optimize_event_selection import EventSelectionOptimizer - - output_path = tmp_path / "cuts.fits" - - with pytest.raises( - ValueError, - match=( - "Need a proton file for cut optimization using " - "PointSourceSensitivityOptimizer" - ), - ): - run_tool( - EventSelectionOptimizer(), - argv=[ - f"--gamma-file={gamma_diffuse_full_reco_file}", - f"--output={output_path}", - f"--config={event_loader_config_path}", - ], - raises=True, - ) - - ret = run_tool( - EventSelectionOptimizer(), - argv=[ - f"--gamma-file={gamma_diffuse_full_reco_file}", - f"--output={output_path}", - f"--config={event_loader_config_path}", - "--EventSelectionOptimizer.optimization_algorithm=PercentileCuts", - ], - ) - assert ret == 0 - assert output_path.exists() diff --git a/src/ctapipe/tools/tests/test_process.py b/src/ctapipe/tools/tests/test_process.py index e1fc0eb6588..424e10d1b9d 100644 --- a/src/ctapipe/tools/tests/test_process.py +++ b/src/ctapipe/tools/tests/test_process.py @@ -88,8 +88,8 @@ def test_stage_1_dl1(tmp_path, dl1_image_file, dl1_parameters_file): # check tables were written with tables.open_file(dl1b_from_dl1a_file, mode="r") as testfile: assert testfile.root.dl1 - assert testfile.root.dl1.events.telescope - assert testfile.root.dl1.events.subarray + assert testfile.root.dl1.event.telescope + assert testfile.root.dl1.event.subarray assert testfile.root.configuration.instrument.subarray.layout assert testfile.root.configuration.instrument.telescope.optics assert testfile.root.configuration.instrument.telescope.camera.geometry_0 @@ -226,7 +226,7 @@ def test_stage_2_from_dl1_images(tmp_path, dl1_image_file): # check tables were written with tables.open_file(output, mode="r") as testfile: - assert testfile.root.dl2.events.subarray.geometry.HillasReconstructor + assert testfile.root.dl2.event.subarray.geometry.HillasReconstructor def test_stage_2_from_dl1_params(tmp_path, dl1_parameters_file): @@ -249,7 +249,7 @@ def test_stage_2_from_dl1_params(tmp_path, dl1_parameters_file): # check tables were written with tables.open_file(output, mode="r") as testfile: - assert testfile.root.dl2.events.subarray.geometry.HillasReconstructor + assert testfile.root.dl2.event.subarray.geometry.HillasReconstructor def test_ml_preprocessing_from_simtel(tmp_path): @@ -274,8 +274,8 @@ def test_ml_preprocessing_from_simtel(tmp_path): # check tables were written with tables.open_file(output, mode="r") as testfile: - assert testfile.root.dl1.events.telescope.parameters.tel_002 - assert testfile.root.dl2.events.subarray.geometry.HillasReconstructor + assert testfile.root.dl1.event.telescope.parameters.tel_002 + assert testfile.root.dl2.event.subarray.geometry.HillasReconstructor def test_image_modifications(tmp_path, dl1_image_file): From 58ecd8237ef2deef5f615a2956dcaf6fc76775d6 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Fri, 21 Mar 2025 14:26:28 +0100 Subject: [PATCH 52/58] Multiple fix --- src/ctapipe/irf/cuts.py | 3 -- src/ctapipe/irf/dl3.py | 72 ++++++++++++++++++-------------- src/ctapipe/irf/preprocessing.py | 67 ++++++++++++++++++++--------- src/ctapipe/tools/create_dl3.py | 19 +++++---- 4 files changed, 98 insertions(+), 63 deletions(-) diff --git a/src/ctapipe/irf/cuts.py b/src/ctapipe/irf/cuts.py index 8e0425995da..64e9866b059 100644 --- a/src/ctapipe/irf/cuts.py +++ b/src/ctapipe/irf/cuts.py @@ -22,9 +22,6 @@ class EventQualitySelection(QualityQuery): "multiplicity 4", "np.count_nonzero(HillasReconstructor_telescopes,axis=1) >= 4", ), - ("valid classifier", "RandomForestClassifier_is_valid"), - ("valid geom reco", "HillasReconstructor_is_valid"), - ("valid energy reco", "RandomForestRegressor_is_valid"), ], help=QualityQuery.quality_criteria.help, ).tag(config=True) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index db74356b45a..4d4e0219198 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -1,5 +1,5 @@ from abc import abstractmethod -from datetime import datetime +from datetime import UTC, datetime from typing import Any, Dict, List, Tuple import astropy.units as u @@ -30,12 +30,12 @@ class DL3_Format(Component): optional_dl3_columns = Bool( default_value=False, help="If true add optional columns to produce file" - ).tag(config=False) + ).tag(config=True) raise_error_for_optional = Bool( default_value=True, help="If true will raise error in the case optional column are missing", - ).tag(config=False) + ).tag(config=True) def __init__(self, **kwargs): super().__init__(**kwargs) @@ -314,7 +314,7 @@ class DL3_GADF(DL3_Format): def __init__(self, **kwargs): super().__init__(**kwargs) - self.file_creation_time = datetime.now(tz=datetime.UTC).isoformat() + self.file_creation_time = datetime.now(tz=UTC) self.reference_time = Time(datetime.fromisoformat("1970-01-01T00:00:00+00:00")) def write_file(self, path): @@ -327,7 +327,7 @@ def write_file(self, path): path : str The full path and filename of the new file to write """ - self.file_creation_time = datetime.now(tz=datetime.UTC).isoformat() + self.file_creation_time = datetime.now(tz=UTC) hdu_dl3 = fits.HDUList( [fits.PrimaryHDU(header=Header(self.get_hdu_header_base_format()))] @@ -381,7 +381,7 @@ def get_hdu_header_base_format(self) -> Dict[str, Any]: "HDUVERS": "v0.3", "HDUDOC": "https://gamma-astro-data-formats.readthedocs.io/en/v0.3/index.html", "CREATOR": "ctapipe " + ctapipe_version, - "CREATED": self.file_creation_time, + "CREATED": self.file_creation_time.isoformat(), } def get_hdu_header_base_time(self) -> Dict[str, Any]: @@ -414,8 +414,8 @@ def get_hdu_header_base_time(self) -> Dict[str, Any]: "TIMEUNIT": "s", "TIMEREF": "GEOCENTER", "TIMESYS": "UTC", - "TSTART": start_time, - "TSTOP": stop_time, + "TSTART": (start_time - self.reference_time).to_value(u.s), + "TSTOP": (stop_time - self.reference_time).to_value(u.s), "ONTIME": ontime.to_value(u.s), "LIVETIME": ontime.to_value(u.s) * self.dead_time_fraction, "DEADC": self.dead_time_fraction, @@ -446,8 +446,16 @@ def get_hdu_header_base_observation_information( object_coordinate = self.target_information[ "object_coordinate" ].transform_to(ICRS()) - header["RA_OBJ"] = object_coordinate.ra.to_value(u.deg) - header["DEC_OBJ"] = object_coordinate.dec.to_value(u.deg) + header["RA_OBJ"] = ( + "nan" + if np.isnan(object_coordinate.ra.to_value(u.deg)) + else object_coordinate.ra.to_value(u.deg) + ) + header["DEC_OBJ"] = ( + "nan" + if np.isnan(object_coordinate.dec.to_value(u.deg)) + else object_coordinate.dec.to_value(u.deg) + ) return header def get_hdu_header_base_subarray_information(self) -> Dict[str, Any]: @@ -460,8 +468,8 @@ def get_hdu_header_base_subarray_information(self) -> Dict[str, Any]: "ORIGIN": self.telescope_information["organisation"], "TELESCOP": self.telescope_information["array"], "INSTRUME": self.telescope_information["subarray"], - "TELLIST": self.telescope_information["telescope_list"], - "N_TELS": np.sum(self.telescope_information["telescope_list"]), + "TELLIST": str(self.telescope_information["telescope_list"]), + "N_TELS": len(self.telescope_information["telescope_list"]), } return header @@ -488,24 +496,25 @@ def get_hdu_header_base_pointing(self) -> Dict[str, Any]: raise ValueError("Telescope location information are missing") gti_table = self.create_gti_table() - time_evaluation = [] + delta_time_evaluation = [] for i in range(len(gti_table)): - time_evaluation += list( + delta_time_evaluation += list( np.linspace(gti_table["START"][i], gti_table["STOP"][i], 100) ) - time_evaluation = self.reference_time + TimeDelta(time_evaluation) + delta_time_evaluation = u.Quantity(delta_time_evaluation) + time_evaluation = self.reference_time + TimeDelta(delta_time_evaluation) pointing_table = self.create_pointing_table() if self.pointing_mode == "TRACK": obs_mode = "POINTING" icrs_coordinate = SkyCoord( ra=np.interp( - time_evaluation, + delta_time_evaluation, xp=pointing_table["TIME"], fp=pointing_table["RA_PNT"], ), dec=np.interp( - time_evaluation, + delta_time_evaluation, xp=pointing_table["TIME"], fp=pointing_table["DEC_PNT"], ), @@ -517,12 +526,12 @@ def get_hdu_header_base_pointing(self) -> Dict[str, Any]: obs_mode = "DRIFT" altaz_coordinate = AltAz( alt=np.interp( - time_evaluation, + delta_time_evaluation, xp=pointing_table["TIME"], fp=pointing_table["ALT_PNT"], ), az=np.interp( - time_evaluation, + delta_time_evaluation, xp=pointing_table["TIME"], fp=pointing_table["AZ_PNT"], ), @@ -543,7 +552,7 @@ def get_hdu_header_base_pointing(self) -> Dict[str, Any]: "AZ_PNT": np.mean(altaz_coordinate.az.to_value(u.deg)), "GEOLON": self.location.lon.to_value(u.deg), "GEOLAT": self.location.lat.to_value(u.deg), - "ALTITUDE": self.location.altitude.to_value(u.m), + "ALTITUDE": self.location.height.to_value(u.m), "OBSGEO-X": self.location.x.to_value(u.m), "OBSGEO-Y": self.location.y.to_value(u.m), "OBSGEO-Z": self.location.z.to_value(u.m), @@ -674,13 +683,13 @@ def create_gti_table(self) -> QTable: (gti_interval[1] - self.reference_time).to(u.s) ) - table = QTable(table_structure).sort("START") - for i in range(len(QTable) - 1): + QTable(table_structure).sort("START") + for i in range(len(table_structure) - 1): if table_structure["STOP"][i] > table_structure["START"][i + 1]: self.log.warning("Overlapping GTI intervals") break - return table + return QTable(table_structure) def create_pointing_table(self) -> QTable: """ @@ -705,13 +714,12 @@ def create_pointing_table(self) -> QTable: pointing_altaz = pointing[1].transform_to( AltAz(location=location, obstime=time) ) - table_structure = { - "TIME": (time - self.reference_time).to(u.s), - "RA_PNT": pointing_icrs.ra.to(u.deg), - "DEC_PNT": pointing_icrs.dec.to(u.deg), - "ALT_PNT": pointing_altaz.alt.to(u.deg), - "AZ_PNT": pointing_altaz.az.to(u.deg), - } - - table = QTable(table_structure).sort("TIME") + table_structure["TIME"].append((time - self.reference_time).to(u.s)) + table_structure["RA_PNT"].append(pointing_icrs.ra.to(u.deg)) + table_structure["DEC_PNT"].append(pointing_icrs.dec.to(u.deg)) + table_structure["ALT_PNT"].append(pointing_altaz.alt.to(u.deg)) + table_structure["AZ_PNT"].append(pointing_altaz.az.to(u.deg)) + + table = QTable(table_structure) + table.sort("TIME") return table diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index b92b6289102..e4d7f5bbc15 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -60,16 +60,16 @@ class EventPreprocessor(Component): irf_pre_processing = Bool( default_value=True, help="If true the pre processing assume the purpose is for IRF production, if false, for DL3 production", - ).tag(config=False) + ).tag(config=True) optional_dl3_columns = Bool( default_value=False, help="If true add optional columns to produce file" - ).tag(config=False) + ).tag(config=True) raise_error_for_optional = Bool( default_value=True, help="If true will raise error in the case optional column are missing", - ).tag(config=False) + ).tag(config=True) def __init__( self, quality_selection_only: bool = True, config=None, parent=None, **kwargs @@ -125,7 +125,8 @@ def get_columns_keep_rename_scheme( "reco_glon", "reco_glat", "reco_fov_lat", - "reco_fov_lon" "reco_source_fov_offset", + "reco_fov_lon", + "reco_source_fov_offset", "reco_source_fov_position_angle", "reco_energy_uncert", "reco_dir_uncert", @@ -301,12 +302,12 @@ def make_derived_columns( ) else: reco_icrs = reco.transform_to(ICRS()) - events["reco_ra"] = reco_icrs.ra - events["reco_dec"] = reco_icrs.dec + events["reco_ra"] = u.Quantity(reco_icrs.ra) + events["reco_dec"] = u.Quantity(reco_icrs.dec) if self.optional_dl3_columns: reco_gal = reco_icrs.transform_to(Galactic()) - events["reco_glon"] = reco_gal.l - events["reco_glat"] = reco_gal.b + events["reco_glon"] = u.Quantity(reco_gal.l) + events["reco_glat"] = u.Quantity(reco_gal.b) events["multiplicity"] = np.sum(events["tels_with_trigger"], axis=1) @@ -325,6 +326,11 @@ def make_empty_table(self, columns_to_use: list[str]) -> QTable: columns = [ Column(name="obs_id", dtype=np.uint64, description="Observation block ID"), Column(name="event_id", dtype=np.uint64, description="Array event ID"), + Time( + val=[], + scale="tai", + format="mjd", + ), Column( name="true_energy", unit=u.TeV, @@ -468,17 +474,27 @@ def make_empty_table(self, columns_to_use: list[str]) -> QTable: # Rearrange in a dict, easier for searching after columns_dict = {} for i in range(len(columns)): - columns_dict[columns[i].name] = columns[i] + if type(columns[i]) is Time: + columns_dict["time"] = columns[i] + else: + columns_dict[columns[i].name] = columns[i] # Select only the necessary columns columns_for_keep = [] - for c in columns_to_use: + index_time_column = -1 + for i, c in enumerate(columns_to_use): if c in columns_dict.keys(): columns_for_keep.append(columns_dict[c]) + if c == "time": + index_time_column = i else: raise ValueError(f"Missing columns definition for {c}") - return QTable(columns_for_keep) + empty_table = QTable(columns_for_keep) + if index_time_column >= 0: + empty_table.rename_column("col" + str(index_time_column), "time") + + return empty_table class EventLoader(Component): @@ -561,10 +577,15 @@ def get_observation_information(self) -> Dict[str, Any]: # Extract GTI list_gti = [] + mask_nan = np.isnan(obs_all_info_table["actual_duration"]) + if np.sum(mask_nan) > 0: + self.log.warning("Duration of the run is nan, replaced with zero") + obs_all_info_table["actual_duration"][mask_nan] = 0.0 * u.s for i in range(len(obs_all_info_table)): start_time = Time(obs_all_info_table["actual_start_time"][i]) stop_time = start_time + TimeDelta( obs_all_info_table["actual_duration"][i] + * obs_all_info_table["actual_duration"].unit ) list_gti.append((start_time, stop_time)) meta["gti"] = list_gti @@ -592,14 +613,18 @@ def get_observation_information(self) -> Dict[str, Any]: == "ALTAZ" ): pointing_start = AltAz( - alt=obs_all_info_table["subarray_pointing_lat"][i], - az=obs_all_info_table["subarray_pointing_lon"][i], + alt=obs_all_info_table["subarray_pointing_lat"][i] + * obs_all_info_table["subarray_pointing_lat"].unit, + az=obs_all_info_table["subarray_pointing_lon"][i] + * obs_all_info_table["subarray_pointing_lon"].unit, location=meta["location"], obstime=meta["gti"][i][0], ) pointing_stop = AltAz( - alt=obs_all_info_table["subarray_pointing_lat"][i], - az=obs_all_info_table["subarray_pointing_lon"][i], + alt=obs_all_info_table["subarray_pointing_lat"][i] + * obs_all_info_table["subarray_pointing_lat"].unit, + az=obs_all_info_table["subarray_pointing_lon"][i] + * obs_all_info_table["subarray_pointing_lon"].unit, location=meta["location"], obstime=meta["gti"][i][1], ) @@ -610,8 +635,10 @@ def get_observation_information(self) -> Dict[str, Any]: == "ICRS" ): pointing_start = ICRS( - dec=obs_all_info_table["subarray_pointing_lat"][i], - ra=obs_all_info_table["subarray_pointing_lon"][i], + dec=obs_all_info_table["subarray_pointing_lat"][i] + * obs_all_info_table["subarray_pointing_lat"].unit, + ra=obs_all_info_table["subarray_pointing_lon"][i] + * obs_all_info_table["subarray_pointing_lon"].unit, ) pointing_stop = pointing_start elif ( @@ -621,8 +648,10 @@ def get_observation_information(self) -> Dict[str, Any]: == "GALACTIC" ): pointing_start = Galactic( - b=obs_all_info_table["subarray_pointing_lat"][i], - l=obs_all_info_table["subarray_pointing_lon"][i], + b=obs_all_info_table["subarray_pointing_lat"][i] + * obs_all_info_table["subarray_pointing_lat"].unit, + l=obs_all_info_table["subarray_pointing_lon"][i] + * obs_all_info_table["subarray_pointing_lon"].unit, ) pointing_stop = pointing_start list_pointing.append((meta["gti"][i][0], pointing_start)) diff --git a/src/ctapipe/tools/create_dl3.py b/src/ctapipe/tools/create_dl3.py index ff1e7d3b639..e8933826d98 100644 --- a/src/ctapipe/tools/create_dl3.py +++ b/src/ctapipe/tools/create_dl3.py @@ -68,13 +68,7 @@ class DL3Tool(Tool): flags = { **flag( - "overwrite", - "DL3_GADF.overwrite", - "Will allow to overwrite existing DL3 file", - "Prevent overwriting existing DL3 file", - ), - **flag( - "optional-column", + "optional-columns", "DL3Tool.optional_dl3_columns", "Add optional columns for events in the DL3 file", "Do not add optional column for events in the DL3 file", @@ -100,6 +94,7 @@ def setup(self): # Setting the GADF format object DL3_Format.optional_dl3_columns = self.optional_dl3_columns DL3_Format.raise_error_for_optional = self.raise_error_for_optional + DL3_Format.overwrite = self.overwrite self.dl3_format = DL3_GADF() @@ -128,10 +123,16 @@ def start(self): for i in range(1, len(hdu_irfs)): if "HDUCLAS2" in hdu_irfs[i].header.keys(): if hdu_irfs[i].header["HDUCLAS2"] == "EFF_AREA": - self.dl3_format.aeff = hdu_irfs[i] + if self.dl3_format.aeff is None: + self.dl3_format.aeff = hdu_irfs[i] + elif "EXTNAME" in hdu_irfs[i].header and not ( + "PROTONS" in hdu_irfs[i].header["EXTNAME"] + or "ELECTRONS" in hdu_irfs[i].header["EXTNAME"] + ): + self.dl3_format.aeff = hdu_irfs[i] elif hdu_irfs[i].header["HDUCLAS2"] == "EDISP": self.dl3_format.edisp = hdu_irfs[i] - elif hdu_irfs[i].header["HDUCLAS2"] == "RPSF": + elif hdu_irfs[i].header["HDUCLAS2"] == "PSF": self.dl3_format.psf = hdu_irfs[i] elif hdu_irfs[i].header["HDUCLAS2"] == "BKG": self.dl3_format.bkg = hdu_irfs[i] From 0fde6452c540f6dda343f2ae83537ed87e0ffcb1 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Fri, 21 Mar 2025 15:21:17 +0100 Subject: [PATCH 53/58] Add docstring --- src/ctapipe/irf/preprocessing.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index e4d7f5bbc15..47b653911a6 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -546,6 +546,9 @@ def load_preselected_events(self, chunk_size: int) -> tuple[QTable, int, dict]: return table def get_observation_information(self) -> Dict[str, Any]: + """ + Retrieve all useful information on the observation for DL3 production + """ with TableLoader(self.file, parent=self, **self.opts_loader) as load: # Extract telescope location meta = {} From 3de36aa8bc4ee90f4d0cdc25148e328541490ea1 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Mon, 24 Mar 2025 15:47:29 +0100 Subject: [PATCH 54/58] Multiple fix --- src/ctapipe/irf/__init__.py | 3 + src/ctapipe/irf/cuts/__init__.py | 4 + src/ctapipe/irf/cuts/quality_cuts.py | 58 ++++ .../irf/{cuts.py => cuts/selection_cuts.py} | 59 +--- src/ctapipe/irf/optimize/algorithm.py | 6 +- src/ctapipe/irf/optimize/results.py | 38 ++- ...optimize.py => test_optimize_algorithm.py} | 38 +-- .../optimize/test/test_optimize_results.py | 40 +++ src/ctapipe/tools/compute_irf.py | 6 +- src/ctapipe/tools/process.py | 2 +- .../tools/tests/test_calculate_pixel_stats.py | 129 ++++++++ src/ctapipe/tools/tests/test_compute_irf.py | 284 ++++++++++++++++++ .../tests/test_optimize_event_selection.py | 117 ++++++++ 13 files changed, 669 insertions(+), 115 deletions(-) create mode 100644 src/ctapipe/irf/cuts/__init__.py create mode 100644 src/ctapipe/irf/cuts/quality_cuts.py rename src/ctapipe/irf/{cuts.py => cuts/selection_cuts.py} (60%) rename src/ctapipe/irf/optimize/test/{test_optimize.py => test_optimize_algorithm.py} (66%) create mode 100644 src/ctapipe/irf/optimize/test/test_optimize_results.py create mode 100644 src/ctapipe/tools/tests/test_calculate_pixel_stats.py create mode 100644 src/ctapipe/tools/tests/test_compute_irf.py create mode 100644 src/ctapipe/tools/tests/test_optimize_event_selection.py diff --git a/src/ctapipe/irf/__init__.py b/src/ctapipe/irf/__init__.py index dbbd812e1c1..8270690b35a 100644 --- a/src/ctapipe/irf/__init__.py +++ b/src/ctapipe/irf/__init__.py @@ -18,6 +18,7 @@ check_bins_in_range, make_bins_per_decade, ) +from .cuts import EventQualitySelection, EventSelection from .irfs import ( BackgroundRate2dMaker, EffectiveArea2dMaker, @@ -56,4 +57,6 @@ "FLUX_UNIT", "check_bins_in_range", "make_bins_per_decade", + "EventSelection", + "EventQualitySelection", ] diff --git a/src/ctapipe/irf/cuts/__init__.py b/src/ctapipe/irf/cuts/__init__.py new file mode 100644 index 00000000000..2ade3dc5524 --- /dev/null +++ b/src/ctapipe/irf/cuts/__init__.py @@ -0,0 +1,4 @@ +from .quality_cuts import EventQualitySelection +from .selection_cuts import EventSelection + +__all__ = ["EventSelection", "EventQualitySelection"] diff --git a/src/ctapipe/irf/cuts/quality_cuts.py b/src/ctapipe/irf/cuts/quality_cuts.py new file mode 100644 index 00000000000..70eb18bafcd --- /dev/null +++ b/src/ctapipe/irf/cuts/quality_cuts.py @@ -0,0 +1,58 @@ +from astropy.table import Table + +from ...core import QualityQuery +from ...core.traits import List, Tuple, Unicode + + +class EventQualitySelection(QualityQuery): + """ + Event pre-selection quality criteria for IRF and DL3 computation with different defaults. + """ + + quality_criteria = List( + Tuple(Unicode(), Unicode()), + default_value=[ + ( + "multiplicity 4", + "np.count_nonzero(HillasReconstructor_telescopes,axis=1) >= 4", + ), + ("valid classifier", "RandomForestClassifier_is_valid"), + ("valid geom reco", "HillasReconstructor_is_valid"), + ("valid energy reco", "RandomForestRegressor_is_valid"), + ], + help=QualityQuery.quality_criteria.help, + ).tag(config=True) + + def calculate_selection(self, events: Table): + """ + Add the selection columns to the events, will only compute quality selection + + Parameters + ---------- + events: Table + The table containing the events on which selection need to be applied + + Returns + ------- + Table + events with selection columns added. + """ + return self.calculate_quality_selection(events) + + def calculate_quality_selection(self, events: Table): + """ + Add the selection columns to the events, will only compute quality selection + + Parameters + ---------- + events: Table + The table containing the events on which selection need to be applied + + Returns + ------- + Table + events with selection columns added. + """ + events["selected_quality"] = self.get_table_mask(events) + events["selected"] = events["selected_quality"] + return events diff --git a/src/ctapipe/irf/cuts.py b/src/ctapipe/irf/cuts/selection_cuts.py similarity index 60% rename from src/ctapipe/irf/cuts.py rename to src/ctapipe/irf/cuts/selection_cuts.py index 64e9866b059..f37c3d1c3b9 100644 --- a/src/ctapipe/irf/cuts.py +++ b/src/ctapipe/irf/cuts/selection_cuts.py @@ -3,62 +3,11 @@ from astropy.table import Table from pyirf.cuts import evaluate_binned_cut -from ..core import QualityQuery -from ..core.traits import List, Path, Tuple, Unicode -from .optimize.results import OptimizationResult +from ...core.traits import Path +from ..optimize.results import OptimizationResult +from .quality_cuts import EventQualitySelection -__all__ = ["EventQualitySelection", "EventSelection"] - - -class EventQualitySelection(QualityQuery): - """ - Event pre-selection quality criteria for IRF and DL3 computation with different defaults. - """ - - quality_criteria = List( - Tuple(Unicode(), Unicode()), - default_value=[ - ( - "multiplicity 4", - "np.count_nonzero(HillasReconstructor_telescopes,axis=1) >= 4", - ), - ], - help=QualityQuery.quality_criteria.help, - ).tag(config=True) - - def calculate_selection(self, events: Table): - """ - Add the selection columns to the events, will only compute quality selection - - Parameters - ---------- - events: Table - The table containing the events on which selection need to be applied - - Returns - ------- - Table - events with selection columns added. - """ - return self.calculate_quality_selection(events) - - def calculate_quality_selection(self, events: Table): - """ - Add the selection columns to the events, will only compute quality selection - - Parameters - ---------- - events: Table - The table containing the events on which selection need to be applied - - Returns - ------- - Table - events with selection columns added. - """ - events["selected_quality"] = self.get_table_mask(events) - events["selected"] = events["selected_quality"] - return events +__all__ = ["EventSelection"] class EventSelection(EventQualitySelection): diff --git a/src/ctapipe/irf/optimize/algorithm.py b/src/ctapipe/irf/optimize/algorithm.py index c4fc6a7aa2a..a117d4aea94 100644 --- a/src/ctapipe/irf/optimize/algorithm.py +++ b/src/ctapipe/irf/optimize/algorithm.py @@ -62,7 +62,7 @@ def __call__( Dictionary containing tables of events used for calculating cuts. This has to include "signal" events and can include "background" events. quality_query: ctapipe.irf.EventPreprocessor - ``ctapipe.core.QualityQuery`` subclass containing preselection + ``ctapipe.irf.cuts.EventQualitySelection`` subclass containing preselection criteria for events. clf_prefix: str Prefix of the output from the G/H classifier for which the @@ -209,7 +209,7 @@ def __call__( ) result = OptimizationResult( - quality_query=quality_query, + quality_selection=quality_query, gh_cuts=gh_cuts, clf_prefix=clf_prefix, valid_energy_min=self.reco_energy_min, @@ -335,7 +335,7 @@ def __call__( ) result = OptimizationResult( - quality_query=quality_query, + quality_selection=quality_query, gh_cuts=gh_cuts, clf_prefix=clf_prefix, valid_energy_min=valid_energy[0], diff --git a/src/ctapipe/irf/optimize/results.py b/src/ctapipe/irf/optimize/results.py index 19de853e085..83b1968ac13 100644 --- a/src/ctapipe/irf/optimize/results.py +++ b/src/ctapipe/irf/optimize/results.py @@ -5,9 +5,9 @@ from astropy.io import fits from astropy.table import QTable, Table -from ...core import QualityQuery from ...core.traits import Path from ...irf import ResultValidRange +from ..cuts.quality_cuts import EventQualitySelection class OptimizationResult: @@ -22,22 +22,28 @@ def __init__( gh_cuts: QTable, clf_prefix: str, spatial_selection_table: QTable | None = None, - quality_query: QualityQuery | Sequence | None = None, + quality_selection: EventQualitySelection | Sequence | None = None, ) -> None: - if quality_query: - if isinstance(quality_query, QualityQuery): - if len(quality_query.quality_criteria) == 0: - quality_query.quality_criteria = [ + if quality_selection: + if isinstance(quality_selection, EventQualitySelection): + if len(quality_selection.quality_criteria) == 0: + quality_selection.quality_criteria = [ (" ", " ") ] # Ensures table serialises properly - self.quality_query = quality_query - elif isinstance(quality_query, list): - self.quality_query = QualityQuery(quality_criteria=quality_query) + self.quality_selection = quality_selection + elif isinstance(quality_selection, list): + self.quality_selection = EventQualitySelection( + quality_criteria=quality_selection + ) else: - self.quality_query = QualityQuery(quality_criteria=list(quality_query)) + self.quality_selection = EventQualitySelection( + quality_criteria=list(quality_selection) + ) else: - self.quality_query = QualityQuery(quality_criteria=[(" ", " ")]) + self.quality_selection = EventQualitySelection( + quality_criteria=[(" ", " ")] + ) self.valid_energy = ResultValidRange(min=valid_energy_min, max=valid_energy_max) self.valid_offset = ResultValidRange(min=valid_offset_min, max=valid_offset_max) @@ -52,21 +58,21 @@ def __repr__(self): f"and {len(self.spatial_selection_table)} theta bins valid " f"between {self.valid_offset.min} to {self.valid_offset.max} " f"and {self.valid_energy.min} to {self.valid_energy.max} " - f"with {len(self.quality_query.quality_criteria)} quality criteria>" + f"with {len(self.quality_selection.quality_criteria)} quality criteria>" ) else: return ( f"" + f"with {len(self.quality_selection.quality_criteria)} quality criteria>" ) def write(self, output_name: Path | str, overwrite: bool = False) -> None: """Write an ``OptimizationResult`` to a file in FITS format.""" cut_expr_tab = Table( - rows=self.quality_query.quality_criteria, + rows=self.quality_selection.quality_criteria, names=["name", "cut_expr"], dtype=[np.str_, np.str_], ) @@ -109,14 +115,14 @@ def read(cls, file_name): if (" ", " ") in cut_expr_lst: cut_expr_lst.remove((" ", " ")) - quality_query = QualityQuery(quality_criteria=cut_expr_lst) + quality_query = EventQualitySelection(quality_criteria=cut_expr_lst) gh_cuts = QTable.read(hdul[2]) valid_energy = QTable.read(hdul[3]) valid_offset = QTable.read(hdul[4]) spatial_selection_table = QTable.read(hdul[5]) if len(hdul) > 5 else None return cls( - quality_query=quality_query, + quality_selection=quality_query, valid_energy_min=valid_energy["energy_min"], valid_energy_max=valid_energy["energy_max"], valid_offset_min=valid_offset["offset_min"], diff --git a/src/ctapipe/irf/optimize/test/test_optimize.py b/src/ctapipe/irf/optimize/test/test_optimize_algorithm.py similarity index 66% rename from src/ctapipe/irf/optimize/test/test_optimize.py rename to src/ctapipe/irf/optimize/test/test_optimize_algorithm.py index 0cb032a12bf..e1126c17571 100644 --- a/src/ctapipe/irf/optimize/test/test_optimize.py +++ b/src/ctapipe/irf/optimize/test/test_optimize_algorithm.py @@ -1,47 +1,11 @@ import astropy.units as u import numpy as np import pytest -from astropy.table import QTable -from ctapipe.core import QualityQuery, non_abstract_children +from ctapipe.core import non_abstract_children from ctapipe.irf.optimize.algorithm import CutOptimizerBase -def test_optimization_result(tmp_path, irf_event_loader_test_config): - from ctapipe.irf import ( - EventPreprocessor, - OptimizationResult, - ResultValidRange, - ) - - result_path = tmp_path / "result.h5" - epp = EventPreprocessor(config=irf_event_loader_test_config) - gh_cuts = QTable( - data=[[0.2, 0.8, 1.5] * u.TeV, [0.8, 1.5, 10] * u.TeV, [0.82, 0.91, 0.88]], - names=["low", "high", "cut"], - ) - result = OptimizationResult( - quality_query=epp.event_selection, - gh_cuts=gh_cuts, - clf_prefix="ExtraTreesClassifier", - valid_energy_min=0.2 * u.TeV, - valid_energy_max=10 * u.TeV, - valid_offset_min=0 * u.deg, - valid_offset_max=np.inf * u.deg, - spatial_selection_table=None, - ) - result.write(result_path) - assert result_path.exists() - - loaded = OptimizationResult.read(result_path) - assert isinstance(loaded, OptimizationResult) - assert isinstance(loaded.quality_query, QualityQuery) - assert isinstance(loaded.valid_energy, ResultValidRange) - assert isinstance(loaded.valid_offset, ResultValidRange) - assert isinstance(loaded.gh_cuts, QTable) - assert loaded.clf_prefix == "ExtraTreesClassifier" - - def test_gh_percentile_cut_calculator(): from ctapipe.irf import GhPercentileCutCalculator diff --git a/src/ctapipe/irf/optimize/test/test_optimize_results.py b/src/ctapipe/irf/optimize/test/test_optimize_results.py new file mode 100644 index 00000000000..b9ec78bb9bf --- /dev/null +++ b/src/ctapipe/irf/optimize/test/test_optimize_results.py @@ -0,0 +1,40 @@ +import numpy as np +from astropy import units as u +from astropy.table import QTable + +from ctapipe.irf.cuts import EventQualitySelection + + +def test_optimization_result(tmp_path, irf_event_loader_test_config): + from ctapipe.irf import ( + EventPreprocessor, + OptimizationResult, + ResultValidRange, + ) + + result_path = tmp_path / "result.h5" + epp = EventPreprocessor(config=irf_event_loader_test_config) + gh_cuts = QTable( + data=[[0.2, 0.8, 1.5] * u.TeV, [0.8, 1.5, 10] * u.TeV, [0.82, 0.91, 0.88]], + names=["low", "high", "cut"], + ) + result = OptimizationResult( + quality_selection=epp.event_selection, + gh_cuts=gh_cuts, + clf_prefix="ExtraTreesClassifier", + valid_energy_min=0.2 * u.TeV, + valid_energy_max=10 * u.TeV, + valid_offset_min=0 * u.deg, + valid_offset_max=np.inf * u.deg, + spatial_selection_table=None, + ) + result.write(result_path) + assert result_path.exists() + + loaded = OptimizationResult.read(result_path) + assert isinstance(loaded, OptimizationResult) + assert isinstance(loaded.quality_selection, EventQualitySelection) + assert isinstance(loaded.valid_energy, ResultValidRange) + assert isinstance(loaded.valid_offset, ResultValidRange) + assert isinstance(loaded.gh_cuts, QTable) + assert loaded.clf_prefix == "ExtraTreesClassifier" diff --git a/src/ctapipe/tools/compute_irf.py b/src/ctapipe/tools/compute_irf.py index 15e2e713a02..740125b91a8 100644 --- a/src/ctapipe/tools/compute_irf.py +++ b/src/ctapipe/tools/compute_irf.py @@ -477,7 +477,7 @@ def start(self): if ( loader.epp.event_selection.quality_criteria - != self.opt_result.quality_query.quality_criteria + != self.opt_result.quality_selection.quality_criteria ): self.log.warning( "Quality criteria are different from quality criteria used for " @@ -487,14 +487,14 @@ def start(self): loader.epp.event_selection.to_table(functions=True)[ "criteria", "func" ], - self.opt_result.event_selection.to_table(functions=True)[ + self.opt_result.quality_selection.to_table(functions=True)[ "criteria", "func" ], ) ) loader.epp.event_selection = EventQualitySelection( parent=loader, - quality_criteria=self.opt_result.event_selection.quality_criteria, + quality_criteria=self.opt_result.quality_selection.quality_criteria, ) self.log.debug( diff --git a/src/ctapipe/tools/process.py b/src/ctapipe/tools/process.py index 7ae3d86e52e..fd5a034d7e7 100644 --- a/src/ctapipe/tools/process.py +++ b/src/ctapipe/tools/process.py @@ -265,7 +265,7 @@ def _write_processing_statistics(self): reconstructor_names, reconstructors ): write_table( - reconstructor.quality_query.to_table(functions=True), + reconstructor.quality_selection.to_table(functions=True), self.write.output_path, f"/dl2/service/tel_event_statistics/{reconstructor_name}", append=True, diff --git a/src/ctapipe/tools/tests/test_calculate_pixel_stats.py b/src/ctapipe/tools/tests/test_calculate_pixel_stats.py new file mode 100644 index 00000000000..9e2e9674664 --- /dev/null +++ b/src/ctapipe/tools/tests/test_calculate_pixel_stats.py @@ -0,0 +1,129 @@ +#!/usr/bin/env python3 +""" +Test ctapipe-calculate-pixel-statistics tool +""" + +import pytest +from traitlets.config.loader import Config + +from ctapipe.core import run_tool +from ctapipe.core.tool import ToolConfigurationError +from ctapipe.instrument import SubarrayDescription +from ctapipe.io import read_table +from ctapipe.tools.calculate_pixel_stats import PixelStatisticsCalculatorTool + + +def test_calculate_pixel_stats_tool(tmp_path, dl1_image_file): + """check statistics calculation from pixel-wise image data files""" + + # Create a configuration suitable for the test + tel_id = 3 + config = Config( + { + "PixelStatisticsCalculatorTool": { + "allowed_tels": [tel_id], + "input_column_name": "image", + "output_table_name": "statistics", + }, + "PixelStatisticsCalculator": { + "stats_aggregator_type": [ + ("id", tel_id, "PlainAggregator"), + ], + }, + "PlainAggregator": { + "chunk_size": 1, + }, + } + ) + # Set the output file path + monitoring_file = tmp_path / "monitoring.dl1.h5" + # Run the tool with the configuration and the input file + run_tool( + PixelStatisticsCalculatorTool(config=config), + argv=[ + f"--input_url={dl1_image_file}", + f"--output_path={monitoring_file}", + "--overwrite", + ], + cwd=tmp_path, + raises=True, + ) + # Check that the output file has been created + assert monitoring_file.exists() + # Check that the output file is not empty + assert ( + read_table( + monitoring_file, + path=f"/dl1/monitoring/telescope/statistics/tel_{tel_id:03d}", + )["mean"] + is not None + ) + # Read subarray description from the created monitoring file + subarray = SubarrayDescription.from_hdf(monitoring_file) + # Check for the selected telescope + assert subarray.tel_ids[0] == tel_id + + +def test_tool_config_error(tmp_path, dl1_image_file): + """check tool configuration error""" + + # Run the tool with the configuration and the input file + config = Config( + { + "PixelStatisticsCalculatorTool": { + "allowed_tels": [3], + "input_column_name": "image_charges", + "output_table_name": "statistics", + } + } + ) + # Set the output file path + monitoring_file = tmp_path / "monitoring.dl1.h5" + # Check if ToolConfigurationError is raised + # when the column name of the pixel-wise image data is not correct + with pytest.raises( + ToolConfigurationError, match="Column 'image_charges' not found" + ): + run_tool( + PixelStatisticsCalculatorTool(config=config), + argv=[ + f"--input_url={dl1_image_file}", + f"--output_path={monitoring_file}", + "--StatisticsAggregator.chunk_size=1", + "--overwrite", + ], + cwd=tmp_path, + raises=True, + ) + # Check if ToolConfigurationError is raised + # when the input and output files are the same + with pytest.raises( + ToolConfigurationError, match="Input and output files are same." + ): + run_tool( + PixelStatisticsCalculatorTool(), + argv=[ + f"--input_url={dl1_image_file}", + f"--output_path={dl1_image_file}", + "--overwrite", + ], + cwd=tmp_path, + raises=True, + ) + # Check if ToolConfigurationError is raised + # when the chunk size is larger than the number of events in the input file + with pytest.raises( + ToolConfigurationError, match="Change --StatisticsAggregator.chunk_size" + ): + run_tool( + PixelStatisticsCalculatorTool(), + argv=[ + f"--input_url={dl1_image_file}", + f"--output_path={monitoring_file}", + "--PixelStatisticsCalculatorTool.allowed_tels=3", + "--StatisticsAggregator.chunk_size=2500", + "--overwrite", + ], + cwd=tmp_path, + raises=True, + ) diff --git a/src/ctapipe/tools/tests/test_compute_irf.py b/src/ctapipe/tools/tests/test_compute_irf.py new file mode 100644 index 00000000000..acd35a50aaf --- /dev/null +++ b/src/ctapipe/tools/tests/test_compute_irf.py @@ -0,0 +1,284 @@ +import json +import logging + +import pytest +from astropy.io import fits + +from ctapipe.core import ToolConfigurationError, run_tool + +pytest.importorskip("pyirf") + + +@pytest.fixture(scope="module") +def dummy_cuts_file( + gamma_diffuse_full_reco_file, + proton_full_reco_file, + event_loader_config_path, + irf_tmp_path, +): + from ctapipe.tools.optimize_event_selection import EventSelectionOptimizer + + output_path = irf_tmp_path / "dummy_cuts.fits" + run_tool( + EventSelectionOptimizer(), + argv=[ + f"--gamma-file={gamma_diffuse_full_reco_file}", + f"--proton-file={proton_full_reco_file}", + # Use diffuse gammas weighted to electron spectrum as stand-in + f"--electron-file={gamma_diffuse_full_reco_file}", + f"--output={output_path}", + f"--config={event_loader_config_path}", + ], + ) + return output_path + + +@pytest.mark.parametrize("include_background", (False, True)) +@pytest.mark.parametrize("spatial_selection_applied", (True, False)) +def test_irf_tool( + gamma_diffuse_full_reco_file, + proton_full_reco_file, + event_loader_config_path, + dummy_cuts_file, + tmp_path, + include_background, + spatial_selection_applied, +): + from ctapipe.tools.compute_irf import IrfTool + + output_path = tmp_path / "irf.fits.gz" + output_benchmarks_path = tmp_path / "benchmarks.fits.gz" + + argv = [ + f"--gamma-file={gamma_diffuse_full_reco_file}", + f"--cuts={dummy_cuts_file}", + f"--output={output_path}", + f"--config={event_loader_config_path}", + ] + if spatial_selection_applied: + argv.append("--spatial-selection-applied") + + if include_background: + argv.append(f"--proton-file={proton_full_reco_file}") + # Use diffuse gammas weighted to electron spectrum as stand-in + argv.append(f"--electron-file={gamma_diffuse_full_reco_file}") + else: + argv.append("--no-do-background") + + ret = run_tool(IrfTool(), argv=argv) + assert ret == 0 + + assert output_path.exists() + assert not output_benchmarks_path.exists() + # Readability by gammapy is tested by pyirf tests, so not repeated here + with fits.open(output_path) as hdul: + assert isinstance(hdul["ENERGY DISPERSION"], fits.BinTableHDU) + assert isinstance(hdul["EFFECTIVE AREA"], fits.BinTableHDU) + assert isinstance(hdul["PSF"], fits.BinTableHDU) + if spatial_selection_applied: + assert isinstance(hdul["RAD_MAX"], fits.BinTableHDU) + + if include_background: + assert isinstance(hdul["BACKGROUND"], fits.BinTableHDU) + + output_path.unlink() # Delete output file + + # Include benchmarks + argv.append(f"--benchmark-output={output_benchmarks_path}") + ret = run_tool(IrfTool(), argv=argv) + assert ret == 0 + + assert output_path.exists() + assert output_benchmarks_path.exists() + with fits.open(output_benchmarks_path) as hdul: + assert isinstance(hdul["ENERGY BIAS RESOLUTION"], fits.BinTableHDU) + assert isinstance(hdul["ANGULAR RESOLUTION"], fits.BinTableHDU) + if include_background: + assert isinstance(hdul["SENSITIVITY"], fits.BinTableHDU) + + +def test_irf_tool_no_electrons( + gamma_diffuse_full_reco_file, + proton_full_reco_file, + event_loader_config_path, + dummy_cuts_file, + tmp_path, +): + from ctapipe.tools.compute_irf import IrfTool + + output_path = tmp_path / "irf.fits.gz" + output_benchmarks_path = tmp_path / "benchmarks.fits.gz" + logpath = tmp_path / "test_irf_tool_no_electrons.log" + logger = logging.getLogger("ctapipe.tools.compute_irf") + logger.addHandler(logging.FileHandler(logpath)) + + ret = run_tool( + IrfTool(), + argv=[ + f"--gamma-file={gamma_diffuse_full_reco_file}", + f"--proton-file={proton_full_reco_file}", + f"--cuts={dummy_cuts_file}", + f"--output={output_path}", + f"--benchmark-output={output_benchmarks_path}", + f"--config={event_loader_config_path}", + f"--log-file={logpath}", + ], + ) + assert ret == 0 + assert output_path.exists() + assert output_benchmarks_path.exists() + assert "Estimating background without electron file." in logpath.read_text() + + +def test_irf_tool_only_gammas( + gamma_diffuse_full_reco_file, event_loader_config_path, dummy_cuts_file, tmp_path +): + from ctapipe.tools.compute_irf import IrfTool + + output_path = tmp_path / "irf.fits.gz" + output_benchmarks_path = tmp_path / "benchmarks.fits.gz" + + with pytest.raises( + ValueError, + match="At least a proton file required when specifying `do_background`.", + ): + run_tool( + IrfTool(), + argv=[ + f"--gamma-file={gamma_diffuse_full_reco_file}", + f"--cuts={dummy_cuts_file}", + f"--output={output_path}", + f"--benchmark-output={output_benchmarks_path}", + f"--config={event_loader_config_path}", + ], + raises=True, + ) + + ret = run_tool( + IrfTool(), + argv=[ + f"--gamma-file={gamma_diffuse_full_reco_file}", + f"--cuts={dummy_cuts_file}", + f"--output={output_path}", + f"--benchmark-output={output_benchmarks_path}", + f"--config={event_loader_config_path}", + "--no-do-background", + ], + ) + assert ret == 0 + assert output_path.exists() + assert output_benchmarks_path.exists() + + +# TODO: Add test using point-like gammas + + +def test_point_like_irf_no_theta_cut( + gamma_diffuse_full_reco_file, + proton_full_reco_file, + event_loader_config_path, + dummy_cuts_file, + tmp_path, +): + from ctapipe.irf import OptimizationResult + from ctapipe.tools.compute_irf import IrfTool + + gh_cuts_path = tmp_path / "gh_cuts.fits" + cuts = OptimizationResult.read(dummy_cuts_file) + cuts.spatial_selection_table = None + cuts.write(gh_cuts_path) + assert gh_cuts_path.exists() + + output_path = tmp_path / "irf.fits.gz" + output_benchmarks_path = tmp_path / "benchmarks.fits.gz" + + with pytest.raises( + ToolConfigurationError, + match=rf"{gh_cuts_path} does not contain any direction cut", + ): + run_tool( + IrfTool(), + argv=[ + f"--gamma-file={gamma_diffuse_full_reco_file}", + f"--proton-file={proton_full_reco_file}", + # Use diffuse gammas weighted to electron spectrum as stand-in + f"--electron-file={gamma_diffuse_full_reco_file}", + f"--cuts={gh_cuts_path}", + f"--output={output_path}", + f"--benchmark-output={output_benchmarks_path}", + f"--config={event_loader_config_path}", + "--spatial-selection-applied", + ], + raises=True, + ) + + +def test_irf_tool_wrong_cuts( + gamma_diffuse_full_reco_file, proton_full_reco_file, dummy_cuts_file, tmp_path +): + from ctapipe.tools.compute_irf import IrfTool + + output_path = tmp_path / "irf.fits.gz" + output_benchmarks_path = tmp_path / "benchmarks.fits.gz" + + with pytest.raises(RuntimeError): + run_tool( + IrfTool(), + argv=[ + f"--gamma-file={gamma_diffuse_full_reco_file}", + f"--proton-file={proton_full_reco_file}", + # Use diffuse gammas weighted to electron spectrum as stand-in + f"--electron-file={gamma_diffuse_full_reco_file}", + f"--cuts={dummy_cuts_file}", + f"--output={output_path}", + f"--benchmark-output={output_benchmarks_path}", + ], + raises=True, + ) + + config_path = tmp_path / "config.json" + with config_path.open("w") as f: + json.dump( + { + "EventPreprocessor": { + "energy_reconstructor": "ExtraTreesRegressor", + "geometry_reconstructor": "HillasReconstructor", + "gammaness_classifier": "ExtraTreesClassifier", + "EventQualityQuery": { + "quality_criteria": [ + # No criteria for minimum event multiplicity + ("valid classifier", "ExtraTreesClassifier_is_valid"), + ("valid geom reco", "HillasReconstructor_is_valid"), + ("valid energy reco", "ExtraTreesRegressor_is_valid"), + ], + }, + } + }, + f, + ) + + logpath = tmp_path / "test_irf_tool_wrong_cuts.log" + logger = logging.getLogger("ctapipe.tools.compute_irf") + logger.addHandler(logging.FileHandler(logpath)) + + ret = run_tool( + IrfTool(), + argv=[ + f"--gamma-file={gamma_diffuse_full_reco_file}", + f"--proton-file={proton_full_reco_file}", + # Use diffuse gammas weighted to electron spectrum as stand-in + f"--electron-file={gamma_diffuse_full_reco_file}", + f"--cuts={dummy_cuts_file}", + f"--output={output_path}", + f"--benchmark-output={output_benchmarks_path}", + f"--config={config_path}", + f"--log-file={logpath}", + ], + ) + assert ret == 0 + assert output_path.exists() + assert output_benchmarks_path.exists() + assert ( + "Quality criteria are different from quality criteria " + "used for calculating g/h / theta cuts." in logpath.read_text() + ) diff --git a/src/ctapipe/tools/tests/test_optimize_event_selection.py b/src/ctapipe/tools/tests/test_optimize_event_selection.py new file mode 100644 index 00000000000..6a581cf5dee --- /dev/null +++ b/src/ctapipe/tools/tests/test_optimize_event_selection.py @@ -0,0 +1,117 @@ +import logging + +import astropy.units as u +import pytest +from astropy.table import QTable + +from ctapipe.core import run_tool +from ctapipe.irf.cuts import EventQualitySelection + +pytest.importorskip("pyirf") + + +def test_cuts_optimization( + gamma_diffuse_full_reco_file, + proton_full_reco_file, + event_loader_config_path, + tmp_path, +): + from ctapipe.irf import ( + OptimizationResult, + ResultValidRange, + ) + from ctapipe.tools.optimize_event_selection import EventSelectionOptimizer + + output_path = tmp_path / "cuts.fits" + + argv = [ + f"--gamma-file={gamma_diffuse_full_reco_file}", + f"--proton-file={proton_full_reco_file}", + # Use diffuse gammas weighted to electron spectrum as stand-in + f"--electron-file={gamma_diffuse_full_reco_file}", + f"--output={output_path}", + f"--config={event_loader_config_path}", + ] + ret = run_tool(EventSelectionOptimizer(), argv=argv) + assert ret == 0 + + result = OptimizationResult.read(output_path) + assert isinstance(result, OptimizationResult) + assert isinstance(result.quality_selection, EventQualitySelection) + assert isinstance(result.valid_energy, ResultValidRange) + assert isinstance(result.valid_offset, ResultValidRange) + assert isinstance(result.gh_cuts, QTable) + assert result.clf_prefix == "ExtraTreesClassifier" + assert "cut" in result.gh_cuts.colnames + assert isinstance(result.spatial_selection_table, QTable) + assert "cut" in result.spatial_selection_table.colnames + + for c in ["low", "center", "high"]: + assert c in result.gh_cuts.colnames + assert result.gh_cuts[c].unit == u.TeV + assert c in result.spatial_selection_table.colnames + assert result.spatial_selection_table[c].unit == u.TeV + + +def test_cuts_opt_no_electrons( + gamma_diffuse_full_reco_file, + proton_full_reco_file, + event_loader_config_path, + tmp_path, +): + from ctapipe.tools.optimize_event_selection import EventSelectionOptimizer + + output_path = tmp_path / "cuts.fits" + logpath = tmp_path / "test_cuts_opt_no_electrons.log" + logger = logging.getLogger("ctapipe.tools.optimize_event_selection") + logger.addHandler(logging.FileHandler(logpath)) + + ret = run_tool( + EventSelectionOptimizer(), + argv=[ + f"--gamma-file={gamma_diffuse_full_reco_file}", + f"--proton-file={proton_full_reco_file}", + f"--output={output_path}", + f"--config={event_loader_config_path}", + f"--log-file={logpath}", + ], + ) + assert ret == 0 + assert "Optimizing cuts without electron file." in logpath.read_text() + + +def test_cuts_opt_only_gammas( + gamma_diffuse_full_reco_file, event_loader_config_path, tmp_path +): + from ctapipe.tools.optimize_event_selection import EventSelectionOptimizer + + output_path = tmp_path / "cuts.fits" + + with pytest.raises( + ValueError, + match=( + "Need a proton file for cut optimization using " + "PointSourceSensitivityOptimizer" + ), + ): + run_tool( + EventSelectionOptimizer(), + argv=[ + f"--gamma-file={gamma_diffuse_full_reco_file}", + f"--output={output_path}", + f"--config={event_loader_config_path}", + ], + raises=True, + ) + + ret = run_tool( + EventSelectionOptimizer(), + argv=[ + f"--gamma-file={gamma_diffuse_full_reco_file}", + f"--output={output_path}", + f"--config={event_loader_config_path}", + "--EventSelectionOptimizer.optimization_algorithm=PercentileCuts", + ], + ) + assert ret == 0 + assert output_path.exists() From 7914648c22c020a11df4af29e5a332461248057d Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Tue, 25 Mar 2025 11:28:42 +0100 Subject: [PATCH 55/58] Fix configurations --- src/ctapipe/conftest.py | 2 +- src/ctapipe/resources/compute_irf.yaml | 2 +- src/ctapipe/resources/optimize_cuts.yaml | 2 +- src/ctapipe/tools/tests/test_compute_irf.py | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ctapipe/conftest.py b/src/ctapipe/conftest.py index b681004f1f2..c5f56335e01 100644 --- a/src/ctapipe/conftest.py +++ b/src/ctapipe/conftest.py @@ -783,7 +783,7 @@ def irf_event_loader_test_config(): "energy_reconstructor": "ExtraTreesRegressor", "geometry_reconstructor": "HillasReconstructor", "gammaness_classifier": "ExtraTreesClassifier", - "EventQualityQuery": { + "EventQualitySelection": { "quality_criteria": [ ( "multiplicity 4", diff --git a/src/ctapipe/resources/compute_irf.yaml b/src/ctapipe/resources/compute_irf.yaml index 806f73d160b..2fed8528371 100644 --- a/src/ctapipe/resources/compute_irf.yaml +++ b/src/ctapipe/resources/compute_irf.yaml @@ -23,7 +23,7 @@ EventPreprocessor: geometry_reconstructor: "HillasReconstructor" gammaness_classifier: "RandomForestClassifier" - EventQualityQuery: + EventQualitySelection: quality_criteria: - ["multiplicity 4", "np.count_nonzero(HillasReconstructor_telescopes,axis=1) >= 4"] - ["valid classifier", "RandomForestClassifier_is_valid"] diff --git a/src/ctapipe/resources/optimize_cuts.yaml b/src/ctapipe/resources/optimize_cuts.yaml index 2639b1f8373..056b17ad238 100644 --- a/src/ctapipe/resources/optimize_cuts.yaml +++ b/src/ctapipe/resources/optimize_cuts.yaml @@ -17,7 +17,7 @@ EventPreprocessor: geometry_reconstructor: "HillasReconstructor" gammaness_classifier: "RandomForestClassifier" - EventQualityQuery: + EventQualitySelection: quality_criteria: - ["multiplicity 4", "np.count_nonzero(HillasReconstructor_telescopes,axis=1) >= 4"] - ["valid classifier", "RandomForestClassifier_is_valid"] diff --git a/src/ctapipe/tools/tests/test_compute_irf.py b/src/ctapipe/tools/tests/test_compute_irf.py index acd35a50aaf..d6d838dbd63 100644 --- a/src/ctapipe/tools/tests/test_compute_irf.py +++ b/src/ctapipe/tools/tests/test_compute_irf.py @@ -244,7 +244,7 @@ def test_irf_tool_wrong_cuts( "energy_reconstructor": "ExtraTreesRegressor", "geometry_reconstructor": "HillasReconstructor", "gammaness_classifier": "ExtraTreesClassifier", - "EventQualityQuery": { + "EventQualitySelection": { "quality_criteria": [ # No criteria for minimum event multiplicity ("valid classifier", "ExtraTreesClassifier_is_valid"), From 46b5a9b1009a72895ef02b4b04b6986a8ce67403 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Wed, 26 Mar 2025 13:02:01 +0100 Subject: [PATCH 56/58] Add X max calculation --- src/ctapipe/irf/dl3.py | 6 ++- src/ctapipe/irf/preprocessing.py | 64 ++++++++++++++++++++++++++++++-- 2 files changed, 65 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 4d4e0219198..c8b1347b2c8 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -628,6 +628,8 @@ def transform_events_columns_for_gadf_format(self, events: QTable) -> QTable: "reco_core_uncert", "reco_h_max", "reco_h_max_uncert", + "reco_x_max", + "reco_x_max_uncert", ] rename_to_optional = [ "MULTIP", @@ -647,9 +649,11 @@ def transform_events_columns_for_gadf_format(self, events: QTable) -> QTable: "CORE_ERR", "HMAX", "HMAX_ERR", + "XMAX", + "XMAX_ERR", ] - if not self.raise_error_for_optional: + if self.raise_error_for_optional: for i, c in enumerate(rename_from_optional): if c not in events.colnames: self.log.warning( diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index 47b653911a6..e57bc68968c 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -22,12 +22,13 @@ ) from tables import NoSuchNodeError +from ..atmosphere import AtmosphereDensityProfile from ..compat import COPY_IF_NEEDED from ..containers import CoordinateFrameType, PointingMode from ..coordinates import NominalFrame from ..core import Component from ..core.traits import Bool, Unicode -from ..io import TableLoader +from ..io import EventSource, TableLoader from ..version import version as ctapipe_version from .cuts import EventQualitySelection, EventSelection from .spectra import SPECTRA, Spectra @@ -135,9 +136,11 @@ def get_columns_keep_rename_scheme( "reco_core_uncert", "reco_h_max", "reco_h_max_uncert", + "reco_x_max", + "reco_x_max_uncert", ] - if not self.raise_error_for_optional: + if self.raise_error_for_optional: if events is None: raise ValueError( "Require events table to assess existence of optional columns" @@ -250,7 +253,10 @@ def keep_necessary_columns_only(self, events: Table) -> QTable: return QTable(events[keep_columns], copy=COPY_IF_NEEDED) def make_derived_columns( - self, events: QTable, location_subarray: EarthLocation + self, + events: QTable, + location_subarray: EarthLocation, + atmosphere_density_profile: AtmosphereDensityProfile = None, ) -> QTable: """ This function compute all the derived columns necessary either to irf production or DL3 file @@ -316,6 +322,35 @@ def make_derived_columns( + events["reco_core_uncert_y"] ** 2 ) + try: + events[ + "reco_x_max" + ] = atmosphere_density_profile.slant_depth_from_height( + events["reco_h_max"], np.pi * u.rad / 2 - events["reco_alt"] + ).to(u.g / u.cm / u.cm) + uncert_up = ( + atmosphere_density_profile.slant_depth_from_height( + events["reco_h_max"] + events["reco_h_max_uncert"], + np.pi * u.rad / 2 - events["reco_alt"], + ) + - events["reco_x_max"] + ) + uncert_down = events[ + "reco_x_max" + ] - atmosphere_density_profile.slant_depth_from_height( + events["reco_h_max"] + events["reco_h_max_uncert"], + np.pi * u.rad / 2 - events["reco_alt"], + ) + events["reco_x_max_uncert"] = ((uncert_up + uncert_down) / 2.0).to( + u.g / u.cm / u.cm + ) + except Exception as e: + if self.raise_error_for_optional: + raise e + else: + events["reco_x_max"] = np.nan * u.g / u.cm / u.cm + events["reco_x_max_uncert"] = np.nan * u.g / u.cm / u.cm + return events def make_empty_table(self, columns_to_use: list[str]) -> QTable: @@ -469,6 +504,16 @@ def make_empty_table(self, columns_to_use: list[str]) -> QTable: unit=u.m, description="Uncertainty on the reconstructed altitude of the maximum of the shower", ), + Column( + name="reco_x_max", + unit=u.g / u.cm / u.cm, + description="Reconstructed maximum of the shower", + ), + Column( + name="reco_x_max_uncert", + unit=u.g / u.cm / u.cm, + description="Uncertainty on the maximum of the shower", + ), ] # Rearrange in a dict, easier for searching after @@ -525,6 +570,16 @@ def __init__( self.opts_loader = dict(dl2=True, simulated=True, observation_info=True) def load_preselected_events(self, chunk_size: int) -> tuple[QTable, int, dict]: + # Load atmosphere density profile + atmosphere_density_profile = None + if self.epp.optional_dl3_columns: + try: + with EventSource(self.file) as source: + atmosphere_density_profile = source.atmosphere_density_profile + except Exception: + self.log.warning("Unable to read atmosphere density profile") + + # Load event with TableLoader(self.file, parent=self, **self.opts_loader) as load: keep_columns, _, _ = self.epp.get_columns_keep_rename_scheme(None, True) header = self.epp.make_empty_table(keep_columns) @@ -533,8 +588,9 @@ def load_preselected_events(self, chunk_size: int) -> tuple[QTable, int, dict]: chunk_size, **self.opts_loader ): events = self.epp.normalise_column_names(events) + events = self.epp.make_derived_columns( - events, load.subarray.reference_location + events, load.subarray.reference_location, atmosphere_density_profile ) events = self.epp.event_selection.calculate_selection(events) selected = events[events["selected"]] From e099c121007625132071335cdd4dd382340e4aca Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Wed, 26 Mar 2025 14:29:53 +0100 Subject: [PATCH 57/58] Switch time scale to TAI --- src/ctapipe/irf/dl3.py | 34 +++++++++++++++++--------------- src/ctapipe/irf/preprocessing.py | 2 +- 2 files changed, 19 insertions(+), 17 deletions(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index c8b1347b2c8..724e17f2807 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -14,7 +14,7 @@ from ctapipe.compat import COPY_IF_NEEDED from ctapipe.coordinates.tests.test_coordinates import location from ctapipe.core import Component -from ctapipe.core.traits import Bool +from ctapipe.core.traits import AstroTime, Bool from ctapipe.version import version as ctapipe_version @@ -37,6 +37,8 @@ class DL3_Format(Component): help="If true will raise error in the case optional column are missing", ).tag(config=True) + reference_time = AstroTime(default_value=Time("2018-01-01T00:00:00", scale="tai")) + def __init__(self, **kwargs): super().__init__(**kwargs) self._obs_id = None @@ -315,7 +317,7 @@ class DL3_GADF(DL3_Format): def __init__(self, **kwargs): super().__init__(**kwargs) self.file_creation_time = datetime.now(tz=UTC) - self.reference_time = Time(datetime.fromisoformat("1970-01-01T00:00:00+00:00")) + self._reference_time = self.reference_time.tai def write_file(self, path): """ @@ -409,21 +411,21 @@ def get_hdu_header_base_time(self) -> Dict[str, Any]: ) return { - "MJDREFI": int(self.reference_time.mjd), - "MJDREFF": self.reference_time.mjd % 1, + "MJDREFI": int(self._reference_time.mjd), + "MJDREFF": self._reference_time.mjd % 1, "TIMEUNIT": "s", "TIMEREF": "GEOCENTER", - "TIMESYS": "UTC", - "TSTART": (start_time - self.reference_time).to_value(u.s), - "TSTOP": (stop_time - self.reference_time).to_value(u.s), + "TIMESYS": "TAI", + "TSTART": (start_time.tai - self._reference_time).to_value(u.s), + "TSTOP": (stop_time.tai - self._reference_time).to_value(u.s), "ONTIME": ontime.to_value(u.s), "LIVETIME": ontime.to_value(u.s) * self.dead_time_fraction, "DEADC": self.dead_time_fraction, - "TELAPSE": (stop_time - start_time).to_value(u.s), - "DATE-OBS": start_time.fits, - "DATE-BEG": start_time.fits, - "DATE-AVG": (start_time + (stop_time - start_time) / 2.0).fits, - "DATE-END": stop_time.fits, + "TELAPSE": (stop_time.tai - start_time.tai).to_value(u.s), + "DATE-OBS": start_time.tai.fits, + "DATE-BEG": start_time.tai.fits, + "DATE-AVG": (start_time.tai + (stop_time.tai - start_time.tai) / 2.0).fits, + "DATE-END": stop_time.tai.fits, } def get_hdu_header_base_observation_information( @@ -502,7 +504,7 @@ def get_hdu_header_base_pointing(self) -> Dict[str, Any]: np.linspace(gti_table["START"][i], gti_table["STOP"][i], 100) ) delta_time_evaluation = u.Quantity(delta_time_evaluation) - time_evaluation = self.reference_time + TimeDelta(delta_time_evaluation) + time_evaluation = self._reference_time + TimeDelta(delta_time_evaluation) pointing_table = self.create_pointing_table() if self.pointing_mode == "TRACK": @@ -681,10 +683,10 @@ def create_gti_table(self) -> QTable: table_structure = {"START": [], "STOP": []} for gti_interval in self.gti: table_structure["START"].append( - (gti_interval[0] - self.reference_time).to(u.s) + (gti_interval[0].tai - self._reference_time).to(u.s) ) table_structure["STOP"].append( - (gti_interval[1] - self.reference_time).to(u.s) + (gti_interval[1].tai - self._reference_time).to(u.s) ) QTable(table_structure).sort("START") @@ -718,7 +720,7 @@ def create_pointing_table(self) -> QTable: pointing_altaz = pointing[1].transform_to( AltAz(location=location, obstime=time) ) - table_structure["TIME"].append((time - self.reference_time).to(u.s)) + table_structure["TIME"].append((time.tai - self._reference_time).to(u.s)) table_structure["RA_PNT"].append(pointing_icrs.ra.to(u.deg)) table_structure["DEC_PNT"].append(pointing_icrs.dec.to(u.deg)) table_structure["ALT_PNT"].append(pointing_altaz.alt.to(u.deg)) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py index e57bc68968c..e722cd5bddf 100644 --- a/src/ctapipe/irf/preprocessing.py +++ b/src/ctapipe/irf/preprocessing.py @@ -641,7 +641,7 @@ def get_observation_information(self) -> Dict[str, Any]: self.log.warning("Duration of the run is nan, replaced with zero") obs_all_info_table["actual_duration"][mask_nan] = 0.0 * u.s for i in range(len(obs_all_info_table)): - start_time = Time(obs_all_info_table["actual_start_time"][i]) + start_time = Time(obs_all_info_table["actual_start_time"][i]).tai stop_time = start_time + TimeDelta( obs_all_info_table["actual_duration"][i] * obs_all_info_table["actual_duration"].unit From bd20df67a3a4cbfbfa7aa7cb5a677339854c9c36 Mon Sep 17 00:00:00 2001 From: Mathieu de Bony Date: Wed, 26 Mar 2025 14:31:52 +0100 Subject: [PATCH 58/58] Add configurable reference time --- src/ctapipe/irf/dl3.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/irf/dl3.py b/src/ctapipe/irf/dl3.py index 724e17f2807..efe527f73e4 100644 --- a/src/ctapipe/irf/dl3.py +++ b/src/ctapipe/irf/dl3.py @@ -37,7 +37,10 @@ class DL3_Format(Component): help="If true will raise error in the case optional column are missing", ).tag(config=True) - reference_time = AstroTime(default_value=Time("2018-01-01T00:00:00", scale="tai")) + reference_time = AstroTime( + default_value=Time("2018-01-01T00:00:00", scale="tai"), + help="The reference time that will be used in the FITS file", + ).tag(config=True) def __init__(self, **kwargs): super().__init__(**kwargs)