diff --git a/docs/api-reference/io/index.rst b/docs/api-reference/io/index.rst index 667c77a87dc..aaa5c22cf23 100644 --- a/docs/api-reference/io/index.rst +++ b/docs/api-reference/io/index.rst @@ -281,6 +281,9 @@ Reference/API .. automodapi:: ctapipe.io.metadata :no-inheritance-diagram: +.. automodapi:: ctapipe.io.dl2_tables_preprocessing + :no-inheritance-diagram: + .. automodapi:: ctapipe.io.eventsource :no-inheritance-diagram: diff --git a/docs/api-reference/irf/index.rst b/docs/api-reference/irf/index.rst index a8f2cc41f3b..fab2810178f 100644 --- a/docs/api-reference/irf/index.rst +++ b/docs/api-reference/irf/index.rst @@ -33,7 +33,6 @@ Submodules irfs benchmarks binning - preprocessing spectra diff --git a/docs/api-reference/irf/preprocessing.rst b/docs/api-reference/irf/preprocessing.rst deleted file mode 100644 index 9d57445cbe3..00000000000 --- a/docs/api-reference/irf/preprocessing.rst +++ /dev/null @@ -1,12 +0,0 @@ -.. _preprocessing: - -******************************* -Event Loading and Preprocessing -******************************* - - -Reference/ API -============== - -.. automodapi:: ctapipe.irf.preprocessing - :no-inheritance-diagram: diff --git a/docs/changes/2791.feature.rst b/docs/changes/2791.feature.rst new file mode 100644 index 00000000000..7f8e6fecbcb --- /dev/null +++ b/docs/changes/2791.feature.rst @@ -0,0 +1,2 @@ +Generalise the DL2 table processing implemented in IRF tools in order to be used by other usecases that ingest, filter and merge DL2 tables. +The module was moved from IRF to IO. diff --git a/src/ctapipe/conftest.py b/src/ctapipe/conftest.py index ab0c3b5ce32..2f3de4da273 100644 --- a/src/ctapipe/conftest.py +++ b/src/ctapipe/conftest.py @@ -905,11 +905,11 @@ def irf_event_loader_test_config(): return Config( { - "EventPreprocessor": { + "DL2EventPreprocessor": { "energy_reconstructor": "ExtraTreesRegressor", "geometry_reconstructor": "HillasReconstructor", "gammaness_classifier": "ExtraTreesClassifier", - "EventQualityQuery": { + "DL2EventQualityQuery": { "quality_criteria": [ ( "multiplicity 4", @@ -936,15 +936,19 @@ def event_loader_config_path(irf_event_loader_test_config, irf_tmp_path): @pytest.fixture(scope="session") def irf_events_table(): - from ctapipe.irf import EventPreprocessor + from ctapipe.io import DL2EventPreprocessor N1 = 1000 N2 = 100 N = N1 + N2 - epp = EventPreprocessor() + epp = DL2EventPreprocessor() tab = epp.make_empty_table() - ids, bulk, unitless = tab.colnames[:2], tab.colnames[2:-2], tab.colnames[-2:] + ids = ["obs_id", "event_id"] + unitless = set( + [colname for colname in tab.colnames if tab[colname].unit is None] + ) - set(ids) + bulk = set(tab.colnames) - set(ids) - set(unitless) id_tab = QTable( data=np.zeros((N, len(ids)), dtype=np.uint64), @@ -956,16 +960,21 @@ def irf_events_table(): names=bulk, units={c: tab[c].unit for c in bulk}, ) - # Setting values following pyirf test in pyirf/irf/tests/test_background.py - bulk_tab["reco_energy"] = np.append(np.full(N1, 1), np.full(N2, 2)) * u.TeV - bulk_tab["true_energy"] = np.append(np.full(N1, 0.9), np.full(N2, 2.1)) * u.TeV - bulk_tab["reco_source_fov_offset"] = ( - np.append(np.full(N1, 0.1), np.full(N2, 0.05)) * u.deg + bulk_tab.replace_column( + "reco_energy", np.append(np.full(N1, 1), np.full(N2, 2)) * u.TeV + ) + bulk_tab.replace_column( + "true_energy", np.append(np.full(N1, 0.9), np.full(N2, 2.1)) * u.TeV ) - bulk_tab["true_source_fov_offset"] = ( - np.append(np.full(N1, 0.11), np.full(N2, 0.04)) * u.deg + bulk_tab.replace_column( + "reco_source_fov_offset", np.append(np.full(N1, 0.1), np.full(N2, 0.05)) * u.deg ) + bulk_tab.replace_column( + "true_source_fov_offset", + np.append(np.full(N1, 0.11), np.full(N2, 0.04)) * u.deg, + ) + for name in unitless: bulk_tab.add_column( Column(name=name, unit=tab[name].unit, data=np.zeros(N) * np.nan) @@ -975,3 +984,40 @@ def irf_events_table(): ev = vstack([e_tab, tab], join_type="exact", metadata_conflicts="silent") return ev + + +@pytest.fixture(scope="function") +def test_config(): + return { + "DL2EventLoader": {"event_reader_function": "read_telescope_events_chunked"}, + "DL2EventPreprocessor": { + "energy_reconstructor": "ExtraTreesRegressor", + "gammaness_classifier": "ExtraTreesClassifier", + "columns_to_rename": {}, + "output_table_schema": [ + Column( + name="obs_id", dtype=np.uint64, description="Observation Block ID" + ), + Column(name="event_id", dtype=np.uint64, description="Array event ID"), + Column(name="tel_id", dtype=np.uint64, description="Telescope ID"), + Column( + name="ExtraTreesRegressor_tel_energy", + unit=u.TeV, + description="Reconstructed energy", + ), + Column( + name="ExtraTreesRegressor_tel_energy_uncert", + unit=u.TeV, + description="Reconstructed energy uncertainty", + ), + ], + "apply_derived_columns": False, + # "disable_column_renaming": True, + "allow_unsupported_pointing_frames": True, + }, + "DL2EventQualityQuery": { + "quality_criteria": [ + ("valid reco", "ExtraTreesRegressor_tel_is_valid"), + ] + }, + } diff --git a/src/ctapipe/io/__init__.py b/src/ctapipe/io/__init__.py index e3a43ec3f34..e0e5db26ad3 100644 --- a/src/ctapipe/io/__init__.py +++ b/src/ctapipe/io/__init__.py @@ -6,6 +6,7 @@ """ from .astropy_helpers import read_table, write_table # noqa: I001 from .datalevels import DataLevel +from .dl2_tables_preprocessing import DL2EventPreprocessor, DL2EventLoader from .eventsource import EventSource from .eventseeker import EventSeeker from .tableio import TableReader, TableWriter @@ -41,5 +42,7 @@ "DataWriter", "DATA_MODEL_VERSION", "get_hdf5_datalevels", + "DL2EventPreprocessor", + "DL2EventLoader", "get_hdf5_monitoring_types", ] diff --git a/src/ctapipe/io/dl2_tables_preprocessing.py b/src/ctapipe/io/dl2_tables_preprocessing.py new file mode 100644 index 00000000000..6838ad700d4 --- /dev/null +++ b/src/ctapipe/io/dl2_tables_preprocessing.py @@ -0,0 +1,476 @@ +"""Module containing classes related to event loading and preprocessing""" + +from pathlib import Path + +import astropy.units as u +import numpy as np +from astropy.coordinates import AltAz, SkyCoord +from astropy.table import Column, QTable, Table, vstack + +try: + from pyirf.simulations import SimulatedEventsInfo + from pyirf.spectral import ( + DIFFUSE_FLUX_UNIT, + POINT_SOURCE_FLUX_UNIT, + PowerLaw, + calculate_event_weights, + ) + from pyirf.utils import calculate_source_fov_offset, calculate_theta + + has_pyirf = True +except ModuleNotFoundError: + has_pyirf = False + +from tables import NoSuchNodeError +from traitlets import default + +from ..compat import COPY_IF_NEEDED +from ..containers import CoordinateFrameType +from ..coordinates import NominalFrame +from ..core import Component, QualityQuery +from ..core.traits import Bool, Dict, List, Tuple, Unicode +from .tableloader import TableLoader + +__all__ = ["DL2EventLoader", "DL2EventPreprocessor", "DL2EventQualityQuery"] + + +class DL2EventQualityQuery(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) + + +class DL2EventPreprocessor(Component): + """Defines pre-selection cuts and the necessary renaming of columns.""" + + classes = [DL2EventQualityQuery] + + energy_reconstructor = Unicode( + default_value="RandomForestRegressor", + help="Prefix of the reco `_energy` column", + ).tag(config=True) + + geometry_reconstructor = Unicode( + default_value="HillasReconstructor", + help="Prefix of the `_alt` and `_az` reco geometry columns", + ).tag(config=True) + + gammaness_classifier = Unicode( + default_value="RandomForestClassifier", + help="Prefix of the classifier `_prediction` column", + ).tag(config=True) + + apply_derived_columns = Bool( + default_value=True, help="Whether to compute derived columns" + ).tag(config=True) + + allow_unsupported_pointing_frames = Bool( + default_value=False, + help=( + "Check whether the pointing is supported." + "For the moment, only pointing in altaz is supported." + "Divergent pointing is also not supported." + ), + ).tag(config=True) + + columns_to_rename = Dict( + key_trait=Unicode(), + value_trait=Unicode(), + help=( + "Dictionary of columns to rename. " + "Leave unset to apply default renaming. " + "Set to an empty dictionary to disable renaming entirely. " + "Set to a partial dictionary to override only some names." + ), + ).tag(config=True) + + output_table_schema = List( + default_value=[ + Column(name="obs_id", dtype=np.uint64, description="Observation Block ID"), + Column(name="event_id", dtype=np.uint64, description="Array event ID"), + Column(name="true_energy", unit=u.TeV, description="Simulated energy"), + Column(name="true_az", unit=u.deg, description="Simulated azimuth"), + Column(name="true_alt", unit=u.deg, description="Simulated altitude"), + Column(name="reco_energy", unit=u.TeV, description="Reconstructed energy"), + Column(name="reco_az", unit=u.deg, description="Reconstructed azimuth"), + Column(name="reco_alt", unit=u.deg, description="Reconstructed altitude"), + Column(name="pointing_az", unit=u.deg, description="Pointing azimuth"), + Column(name="pointing_alt", unit=u.deg, description="Pointing altitude"), + Column( + name="gh_score", + dtype=np.float64, + description="prediction of the classifier, defined between [0,1]," + " where values close to 1 mean that the positive class" + " (e.g. gamma in gamma-ray analysis) is more likely", + ), + ], + help="Schema definition for output event QTable", + ).tag(config=True) + + def __init__(self, config=None, parent=None, **kwargs): + super().__init__(config=config, parent=parent, **kwargs) + self.quality_query = DL2EventQualityQuery(parent=self) + + @default("columns_to_rename") + def _default_columns_to_rename(self): + return { + f"{self.energy_reconstructor}_energy": "reco_energy", + f"{self.geometry_reconstructor}_az": "reco_az", + f"{self.geometry_reconstructor}_alt": "reco_alt", + f"{self.gammaness_classifier}_prediction": "gh_score", + "subarray_pointing_lat": "pointing_alt", + "subarray_pointing_lon": "pointing_az", + } + + def normalise_column_names(self, events: QTable) -> QTable: + """ + Rename column names according to configuration. + + Parameters + ---------- + events : QTable + Input event table. + + Returns + ------- + QTable + Table with selected and renamed columns. + + Raises + ------ + NotImplementedError + If pointing is not AltAz or varies too much. + ValueError + If required columns are missing. + """ + if not self.allow_unsupported_pointing_frames: + if 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." + ) + + columns_to_keep = [col.name for col in self.output_table_schema] + + rename_dict = self.columns_to_rename + rename_from = list(rename_dict.keys()) + rename_to = list(rename_dict.values()) + + fixed_columns = list(set(columns_to_keep) - set(rename_to)) + columns_to_read = fixed_columns + rename_from + for col in columns_to_read: + if col not in events.colnames: + raise ValueError( + f"Input files must conform to the ctapipe DL2 data model. " + f"Required column {col} is missing." + ) + + events = QTable(events[columns_to_read], copy=COPY_IF_NEEDED) + if rename_from and rename_to: + events.rename_columns(rename_from, rename_to) + return events + + def make_empty_table(self) -> QTable: + """ + Create an empty event table based on the configured output schema. + """ + schema = list( + self.output_table_schema + ) # make a shallow copy to extend the schema with derived columns + + if self.apply_derived_columns: + schema.extend( + [ + Column( + name="reco_fov_lat", + unit=u.deg, + description="Reconstructed FOV lat", + ), + Column( + name="reco_fov_lon", + unit=u.deg, + description="Reconstructed FOV lon", + ), + Column( + name="theta", + unit=u.deg, + description="Angular offset from source", + ), + Column( + name="true_source_fov_offset", + unit=u.deg, + description="Simulated angular offset from pointing direction", + ), + Column( + name="reco_source_fov_offset", + unit=u.deg, + description="Reconstructed angular offset from pointing direction", + ), + Column( + name="weight", + dtype=np.float64, + description="Event weight", + ), + ] + ) + + return QTable( + names=[col.name for col in schema], + dtype=[col.dtype for col in schema], + units=[col.unit for col in schema] + if any(col.unit for col in schema) + else None, + meta={}, + ) + + +class DL2EventLoader(Component): + """ + Component for loading events and simulation metadata, applying preselection and optional derived column logic. + """ + + classes = [DL2EventPreprocessor] + + # User-selectable event reading function and kwargs + event_reader_function = Unicode( + default_value="read_subarray_events_chunked", + help=( + "Function of TableLoader used to read event chunks. " + "E.g., 'read_subarray_events_chunked' or 'read_telescope_events_chunked'." + ), + ).tag(config=True) + + event_reader_kwargs = Dict( + default_value={}, + help="Extra keyword arguments passed to the event reading function, e.g., {'path': '/dl2/event/telescope/Reconstructor'}", + ).tag(config=True) + + def __init__(self, file: Path, target_spectrum: "Spectra", **kwargs): # noqa: F821 + from ..irf.spectra import SPECTRA + + super().__init__(**kwargs) + self.epp = DL2EventPreprocessor(parent=self) + 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]: + """ + Load and filter events from the file. + + Parameters + ---------- + chunk_size : int + Size of chunks to read from the file. + obs_time : Quantity + Observation time to scale weights. + + Returns + ------- + table : QTable + Filtered and processed event table. + n_raw_events : int + Number of events before selection. + meta : dict + Metadata dictionary with simulation info and input spectrum. + """ + + opts = dict(dl2=True, simulated=True, observation_info=True) + + with TableLoader(self.file, parent=self, **opts) as loader: + table_template = self.epp.make_empty_table() + sim_info, spectrum = self.get_simulation_information(loader, obs_time) + meta = {"sim_info": sim_info, "spectrum": spectrum} + event_chunks = [table_template] + n_raw_events = 0 + reader_func = getattr(loader, self.event_reader_function) + table_reader = reader_func(chunk_size, **opts, **self.event_reader_kwargs) + for _, _, events in table_reader: + selected = events[self.epp.quality_query.get_table_mask(events)] + selected = self.epp.normalise_column_names(selected) + if self.epp.apply_derived_columns: + selected = self.make_derived_columns(selected) + event_chunks.append(selected) + n_raw_events += len(events) + + event_chunks.append( + table_template + ) # Putting it last ensures the correct metadata is used + table = vstack(event_chunks, join_type="exact", metadata_conflicts="silent") + return table, n_raw_events, meta + + def get_simulation_information( + self, loader: TableLoader, obs_time: u.Quantity + ) -> tuple["SimulatedEventsInfo", "PowerLaw"]: + """ + Extract simulation information from the input file. + + Parameters + ---------- + loader : TableLoader + Loader object for reading from the input file. + obs_time : Quantity + Total observation time. + + Returns + ------- + sim_info : SimulatedEventsInfo + Metadata about the simulated events. + spectrum : PowerLaw + Power-law model derived from simulation configuration. + + Raises + ------ + NotImplementedError + If simulation parameters vary across runs. + """ + from ..exceptions import OptionalDependencyMissing + + if not has_pyirf: + raise OptionalDependencyMissing("pyirf") + + sim = loader.read_simulation_configuration() + try: + show = loader.read_shower_distribution() + except NoSuchNodeError: + self.log.warning( + "Simulation distributions were not found in the input files, falling back to estimating the number of showers from the simulation configuration." + ) + 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) + + def make_derived_columns(self, events: QTable) -> QTable: + """ + Add derived quantities (e.g., theta, FOV offsets) to the table. + + Parameters + ---------- + events : QTable + Table containing normalized events. + + Returns + ------- + QTable + Table with added derived columns. + """ + 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" + ) + + 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) + events["weight"] = 1.0 # defer calculation of proper weights to later + return events + + def make_event_weights( + self, + events: QTable, + spectrum: "PowerLaw", + kind: str, + fov_offset_bins: u.Quantity | None = None, + ) -> QTable: + """ + Compute event weights to match the target spectrum. + + Parameters + ---------- + events : QTable + Input events. + spectrum : PowerLaw + Spectrum from simulation. + kind : str + Type of events ("gammas", etc.). + fov_offset_bins : Quantity, optional + Offset bins for integrating the diffuse flux into point source bins. + + Returns + ------- + QTable + Table with updated weights. + + Raises + ------ + ValueError + If ``fov_offset_bins`` is required but not provided. + """ + if ( + kind == "gammas" + and self.target_spectrum.normalization.unit.is_equivalent( + POINT_SOURCE_FLUX_UNIT + ) + and spectrum.normalization.unit.is_equivalent(DIFFUSE_FLUX_UNIT) + ): + if fov_offset_bins is None: + raise ValueError( + "gamma_target_spectrum is point-like, but no fov offset bins " + "for the integration of the simulated diffuse spectrum were given." + ) + + for low, high in zip(fov_offset_bins[:-1], fov_offset_bins[1:]): + fov_mask = events["true_source_fov_offset"] >= low + fov_mask &= events["true_source_fov_offset"] < high + + events["weight"][fov_mask] = calculate_event_weights( + events[fov_mask]["true_energy"], + target_spectrum=self.target_spectrum, + simulated_spectrum=spectrum.integrate_cone(low, high), + ) + else: + events["weight"] = calculate_event_weights( + events["true_energy"], + target_spectrum=self.target_spectrum, + simulated_spectrum=spectrum, + ) + + return events diff --git a/src/ctapipe/io/tests/test_preprocessing.py b/src/ctapipe/io/tests/test_preprocessing.py new file mode 100644 index 00000000000..46879ac6350 --- /dev/null +++ b/src/ctapipe/io/tests/test_preprocessing.py @@ -0,0 +1,236 @@ +import astropy.units as u +import numpy as np +import pytest +from astropy.table import Column, QTable, Table +from traitlets.config import Config + + +@pytest.fixture(scope="function") +def dummy_table(): + """Dummy table to test column renaming.""" + return Table( + { + "obs_id": [1, 1, 1, 2, 3, 3], + "event_id": [1, 2, 3, 1, 1, 2], + "true_energy": [0.99, 10, 0.37, 2.1, 73.4, 1] * u.TeV, + "dummy_energy": [1, 10, 0.4, 2.5, 73, 1] * u.TeV, + "classifier_prediction": [1, 0.3, 0.87, 0.93, 0, 0.1], + "true_alt": [60, 60, 60, 60, 60, 60] * u.deg, + "geom_alt": [58.5, 61.2, 59, 71.6, 60, 62] * u.deg, + "true_az": [13, 13, 13, 13, 13, 13] * u.deg, + "geom_az": [12.5, 13, 11.8, 15.1, 14.7, 12.8] * u.deg, + "subarray_pointing_frame": np.zeros(6), + "subarray_pointing_lat": np.full(6, 20) * u.deg, + "subarray_pointing_lon": np.full(6, 0) * u.deg, + } + ) + + +def test_normalise_column_names(dummy_table): + from ctapipe.io.dl2_tables_preprocessing import DL2EventPreprocessor + + output_table_schema = [ + Column(name="obs_id", dtype=np.uint64, description="Observation Block ID"), + Column(name="event_id", dtype=np.uint64, description="Array event ID"), + Column(name="true_energy", unit=u.TeV, description="Simulated energy"), + Column(name="true_az", unit=u.deg, description="Simulated azimuth"), + Column(name="true_alt", unit=u.deg, description="Simulated altitude"), + Column(name="reco_energy", unit=u.TeV, description="Reconstructed energy"), + Column(name="reco_az", unit=u.deg, description="Reconstructed azimuth"), + Column(name="reco_alt", unit=u.deg, description="Reconstructed altitude"), + Column(name="pointing_alt", unit=u.deg, description="Pointing latitude"), + Column(name="pointing_az", unit=u.deg, description="Pointing longitude"), + Column( + name="gh_score", + dtype=np.float64, + description="prediction of the classifier, defined between [0,1]," + " where values close to 1 mean that the positive class" + " (e.g. gamma in gamma-ray analysis) is more likely", + ), + ] + epp = DL2EventPreprocessor( + energy_reconstructor="dummy", + geometry_reconstructor="geom", + gammaness_classifier="classifier", + output_table_schema=output_table_schema, + ) + norm_table = epp.normalise_column_names(dummy_table) + + needed_cols = [ + "obs_id", + "event_id", + "true_energy", + "true_alt", + "true_az", + "reco_energy", + "reco_alt", + "reco_az", + "gh_score", + "pointing_alt", + "pointing_az", + ] + for c in needed_cols: + assert c in norm_table.colnames + + with pytest.raises(ValueError, match="Required column geom_alt is missing."): + dummy_table.rename_column("geom_alt", "alt_geom") + epp = DL2EventPreprocessor( + energy_reconstructor="dummy", + geometry_reconstructor="geom", + gammaness_classifier="classifier", + output_table_schema=output_table_schema, + ) + _ = epp.normalise_column_names(dummy_table) + + +def test_event_loader(gamma_diffuse_full_reco_file, irf_event_loader_test_config): + pytest.importorskip("pyirf", reason="pyirf is an optional dependency") + from pyirf.simulations import SimulatedEventsInfo + from pyirf.spectral import PowerLaw + + from ctapipe.io.dl2_tables_preprocessing import DL2EventLoader + from ctapipe.irf import Spectra + + loader = DL2EventLoader( + config=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), + ) + + columns = [ + "obs_id", + "event_id", + "true_energy", + "true_az", + "true_alt", + "reco_energy", + "reco_az", + "reco_alt", + "reco_fov_lat", + "reco_fov_lon", + "gh_score", + "pointing_az", + "pointing_alt", + "theta", + "true_source_fov_offset", + "reco_source_fov_offset", + "weight", + ] + + assert sorted(columns) == sorted(events.colnames) + assert isinstance(count, int) + assert isinstance(meta["sim_info"], SimulatedEventsInfo) + assert isinstance(meta["spectrum"], PowerLaw) + + events = loader.make_event_weights( + events, meta["spectrum"], "gammas", (0 * u.deg, 1 * u.deg) + ) + + assert "weight" in events.colnames + + +def test_preprocessor_tel_table_with_custom_reconstructor(tmp_path, test_config): + from ctapipe.io.dl2_tables_preprocessing import DL2EventPreprocessor + + # Create a test table with required columns + table = QTable( + { + "obs_id": [1, 1, 2], + "event_id": [100, 101, 102], + "tel_id": [1, 1, 1], + "ExtraTreesRegressor_tel_energy": [1.0, 2.0, 3.0] * u.TeV, + "ExtraTreesRegressor_tel_is_valid": [True, False, True], + "ExtraTreesRegressor_tel_energy_uncert": [0.1, 0.2, 0.1], + "ExtraTreesRegressor_tel_goodness_of_fit": [0.9, 0.8, 0.95], + "subarray_pointing_lat": [80.0, 80.0, 80.0] * u.deg, + "subarray_pointing_lon": [0.0, 0.0, 0.0] * u.deg, + "true_energy": [1.1, 2.1, 3.1] * u.TeV, + "true_az": [42.0, 43.0, 44.0] * u.deg, + "true_alt": [70.0, 71.0, 72.0] * u.deg, + } + ) + + # Set up config + config = test_config + + # Create preprocessor with config + preprocessor = DL2EventPreprocessor(config=Config(config)) + + # Apply quality query and preprocessing + mask = preprocessor.quality_query.get_table_mask(table) + filtered = table[mask] + + # Apply renaming and derived column generation + processed = preprocessor.normalise_column_names(filtered) + + # Check expected column names after renaming + assert "ExtraTreesRegressor_tel_energy" in processed.colnames + assert "obs_id" in processed.colnames # might exist depending on classifier config + assert "tel_id" in processed.colnames + + # Check the number of surviving rows (only valid events) + assert len(processed) == 2 + assert np.all(processed["ExtraTreesRegressor_tel_energy"] > 0 * u.TeV) + + +def test_name_overriding(dummy_table): + from ctapipe.io.dl2_tables_preprocessing import DL2EventPreprocessor + + epp = DL2EventPreprocessor( + energy_reconstructor="dummy", + geometry_reconstructor="geom", + gammaness_classifier="classifier", + columns_to_rename={"true_energy": "false_energy"}, + output_table_schema=[ + Column(name="obs_id", dtype=np.uint64, description="Observation Block ID"), + Column(name="event_id", dtype=np.uint64, description="Array event ID"), + Column(name="false_energy", unit=u.TeV, description="Simulated energy"), + Column(name="true_az", unit=u.deg, description="Simulated azimuth"), + Column(name="true_alt", unit=u.deg, description="Simulated altitude"), + Column(name="dummy_energy", unit=u.TeV, description="Reconstructed energy"), + Column(name="geom_az", unit=u.deg, description="Reconstructed azimuth"), + Column(name="geom_alt", unit=u.deg, description="Reconstructed altitude"), + Column( + name="subarray_pointing_frame", + unit=u.dimensionless_unscaled, + description="Pointing frame", + ), + Column( + name="subarray_pointing_lat", + unit=u.deg, + description="Pointing latitude", + ), + Column( + name="subarray_pointing_lon", + unit=u.deg, + description="Pointing longitude", + ), + Column( + name="classifier_prediction", + unit=u.dimensionless_unscaled, + description="prediction of the classifier, defined between [0,1]," + " where values close to 1 mean that the positive class" + " (e.g. gamma in gamma-ray analysis) is more likely", + ), + ], + ) + norm_table = epp.normalise_column_names(dummy_table) + columns = [ + "obs_id", + "event_id", + "false_energy", + "true_az", + "true_alt", + "dummy_energy", + "classifier_prediction", + "geom_alt", + "geom_az", + "subarray_pointing_frame", + "subarray_pointing_lat", + "subarray_pointing_lon", + ] + assert sorted(columns) == sorted(norm_table.colnames) diff --git a/src/ctapipe/irf/__init__.py b/src/ctapipe/irf/__init__.py index dbbd812e1c1..d4abdb43fac 100644 --- a/src/ctapipe/irf/__init__.py +++ b/src/ctapipe/irf/__init__.py @@ -31,7 +31,6 @@ PointSourceSensitivityOptimizer, ThetaPercentileCutCalculator, ) -from .preprocessing import EventLoader, EventPreprocessor from .spectra import ENERGY_FLUX_UNIT, FLUX_UNIT, SPECTRA, Spectra __all__ = [ @@ -46,8 +45,6 @@ "OptimizationResult", "PointSourceSensitivityOptimizer", "PercentileCuts", - "EventLoader", - "EventPreprocessor", "Spectra", "GhPercentileCutCalculator", "ThetaPercentileCutCalculator", diff --git a/src/ctapipe/irf/optimize.py b/src/ctapipe/irf/optimize.py index fd8f9e9ad52..3aa99bf21ac 100644 --- a/src/ctapipe/irf/optimize.py +++ b/src/ctapipe/irf/optimize.py @@ -13,8 +13,8 @@ from ..core import Component, QualityQuery from ..core.traits import AstroQuantity, Float, Integer, Path +from ..io.dl2_tables_preprocessing import DL2EventQualityQuery from .binning import DefaultRecoEnergyBins, ResultValidRange -from .preprocessing import EventQualityQuery __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: DL2EventQualityQuery, clf_prefix: str, ) -> OptimizationResult: """ @@ -180,7 +180,7 @@ def __call__( events: dict[str, astropy.table.QTable] 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 + quality_query: ctapipe.io.DL2EventPreprocessor ``ctapipe.core.QualityQuery`` subclass containing preselection criteria for events. clf_prefix: str @@ -305,7 +305,7 @@ def __init__(self, config=None, parent=None, **kwargs): def __call__( self, events: dict[str, QTable], - quality_query: EventQualityQuery, + quality_query: DL2EventQualityQuery, 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: DL2EventQualityQuery, clf_prefix: str, ) -> OptimizationResult: self._check_events(events) diff --git a/src/ctapipe/irf/preprocessing.py b/src/ctapipe/irf/preprocessing.py deleted file mode 100644 index 21653693793..00000000000 --- a/src/ctapipe/irf/preprocessing.py +++ /dev/null @@ -1,333 +0,0 @@ -"""Module containing classes related to event loading and preprocessing""" - -from pathlib import Path - -import astropy.units as u -import numpy as np -from astropy.coordinates import AltAz, SkyCoord -from astropy.table import Column, QTable, Table, vstack -from pyirf.simulations import SimulatedEventsInfo -from pyirf.spectral import ( - DIFFUSE_FLUX_UNIT, - POINT_SOURCE_FLUX_UNIT, - PowerLaw, - calculate_event_weights, -) -from pyirf.utils import calculate_source_fov_offset, calculate_theta -from tables import NoSuchNodeError - -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 ..io import TableLoader -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) - - -class EventPreprocessor(Component): - """Defines pre-selection cuts and the necessary renaming of columns.""" - - classes = [EventQualityQuery] - - energy_reconstructor = Unicode( - default_value="RandomForestRegressor", - help="Prefix of the reco `_energy` column", - ).tag(config=True) - - geometry_reconstructor = Unicode( - default_value="HillasReconstructor", - help="Prefix of the `_alt` and `_az` reco geometry columns", - ).tag(config=True) - - gammaness_classifier = Unicode( - default_value="RandomForestClassifier", - help="Prefix of the classifier `_prediction` column", - ).tag(config=True) - - def __init__(self, config=None, parent=None, **kwargs): - super().__init__(config=config, parent=parent, **kwargs) - self.quality_query = EventQualityQuery(parent=self) - - def normalise_column_names(self, events: Table) -> QTable: - if 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 = [ - "obs_id", - "event_id", - "true_energy", - "true_az", - "true_alt", - ] - rename_from = [ - f"{self.energy_reconstructor}_energy", - f"{self.geometry_reconstructor}_az", - f"{self.geometry_reconstructor}_alt", - f"{self.gammaness_classifier}_prediction", - "subarray_pointing_lat", - "subarray_pointing_lon", - ] - rename_to = [ - "reco_energy", - "reco_az", - "reco_alt", - "gh_score", - "pointing_alt", - "pointing_az", - ] - keep_columns.extend(rename_from) - for c in keep_columns: - if c not in events.colnames: - raise ValueError( - "Input files must conform to the ctapipe DL2 data model. " - f"Required column {c} is missing." - ) - - events = QTable(events[keep_columns], copy=COPY_IF_NEEDED) - events.rename_columns(rename_from, rename_to) - return events - - def make_empty_table(self) -> QTable: - """ - This function defines the columns later functions expect to be present - in the event table. - """ - columns = [ - Column(name="obs_id", dtype=np.uint64, description="Observation block ID"), - Column(name="event_id", dtype=np.uint64, description="Array event ID"), - Column( - name="true_energy", - unit=u.TeV, - description="Simulated energy", - ), - Column( - name="true_az", - unit=u.deg, - description="Simulated azimuth", - ), - Column( - name="true_alt", - unit=u.deg, - description="Simulated altitude", - ), - Column( - name="reco_energy", - unit=u.TeV, - description="Reconstructed energy", - ), - Column( - name="reco_az", - unit=u.deg, - description="Reconstructed azimuth", - ), - Column( - name="reco_alt", - unit=u.deg, - description="Reconstructed altitude", - ), - Column( - name="reco_fov_lat", - unit=u.deg, - description="Reconstructed field of view lat", - ), - Column( - name="reco_fov_lon", - unit=u.deg, - description="Reconstructed field of view lon", - ), - Column(name="pointing_az", unit=u.deg, description="Pointing azimuth"), - Column(name="pointing_alt", unit=u.deg, description="Pointing altitude"), - Column( - name="theta", - unit=u.deg, - description="Reconstructed angular offset from source position", - ), - Column( - name="true_source_fov_offset", - unit=u.deg, - description="Simulated angular offset from pointing direction", - ), - Column( - name="reco_source_fov_offset", - unit=u.deg, - description="Reconstructed angular offset from pointing direction", - ), - Column( - name="gh_score", - unit=u.dimensionless_unscaled, - description="prediction of the classifier, defined between [0,1]," - " where values close to 1 mean that the positive class" - " (e.g. gamma in gamma-ray analysis) is more likely", - ), - Column( - name="weight", - unit=u.dimensionless_unscaled, - description="Event weight", - ), - ] - - return QTable(columns) - - -class EventLoader(Component): - """ - Contains functions to load events and simulation information from a file - and derive some additional columns needed for irf calculation. - """ - - classes = [EventPreprocessor] - - def __init__(self, file: Path, target_spectrum: Spectra, **kwargs): - super().__init__(**kwargs) - - self.epp = EventPreprocessor(parent=self) - 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]: - 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 - - def get_simulation_information( - self, loader: TableLoader, 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], - ) - - return sim_info, PowerLaw.from_simulation(sim_info, obstime=obs_time) - - 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, - spectrum: PowerLaw, - kind: str, - fov_offset_bins: u.Quantity | None = None, - ) -> QTable: - if ( - kind == "gammas" - and self.target_spectrum.normalization.unit.is_equivalent( - POINT_SOURCE_FLUX_UNIT - ) - and spectrum.normalization.unit.is_equivalent(DIFFUSE_FLUX_UNIT) - ): - if fov_offset_bins is None: - raise ValueError( - "gamma_target_spectrum is point-like, but no fov offset bins " - "for the integration of the simulated diffuse spectrum were given." - ) - - for low, high in zip(fov_offset_bins[:-1], fov_offset_bins[1:]): - fov_mask = events["true_source_fov_offset"] >= low - fov_mask &= events["true_source_fov_offset"] < high - - events["weight"][fov_mask] = calculate_event_weights( - events[fov_mask]["true_energy"], - target_spectrum=self.target_spectrum, - simulated_spectrum=spectrum.integrate_cone(low, high), - ) - else: - events["weight"] = calculate_event_weights( - events["true_energy"], - target_spectrum=self.target_spectrum, - simulated_spectrum=spectrum, - ) - - return events diff --git a/src/ctapipe/irf/tests/test_benchmarks.py b/src/ctapipe/irf/tests/test_benchmarks.py index 8264003c9ee..b8957bdb978 100644 --- a/src/ctapipe/irf/tests/test_benchmarks.py +++ b/src/ctapipe/irf/tests/test_benchmarks.py @@ -83,10 +83,11 @@ def test_make_2d_ang_res(irf_events_table): def test_make_2d_sensitivity( gamma_diffuse_full_reco_file, proton_full_reco_file, irf_event_loader_test_config ): - from ctapipe.irf import EventLoader, Sensitivity2dMaker, Spectra + from ctapipe.io import DL2EventLoader + from ctapipe.irf import Sensitivity2dMaker, Spectra from ctapipe.irf.tests.test_irfs import _check_boundaries_in_hdu - gamma_loader = EventLoader( + gamma_loader = DL2EventLoader( config=irf_event_loader_test_config, file=gamma_diffuse_full_reco_file, target_spectrum=Spectra.CRAB_HEGRA, @@ -95,7 +96,7 @@ def test_make_2d_sensitivity( chunk_size=10000, obs_time=u.Quantity(50, u.h), ) - proton_loader = EventLoader( + proton_loader = DL2EventLoader( config=irf_event_loader_test_config, file=proton_full_reco_file, target_spectrum=Spectra.IRFDOC_PROTON_SPECTRUM, diff --git a/src/ctapipe/irf/tests/test_optimize.py b/src/ctapipe/irf/tests/test_optimize.py index 366c7ca334a..fb5b99694d0 100644 --- a/src/ctapipe/irf/tests/test_optimize.py +++ b/src/ctapipe/irf/tests/test_optimize.py @@ -8,14 +8,14 @@ def test_optimization_result(tmp_path, irf_event_loader_test_config): + from ctapipe.io import DL2EventPreprocessor from ctapipe.irf import ( - EventPreprocessor, OptimizationResult, ResultValidRange, ) result_path = tmp_path / "result.h5" - epp = EventPreprocessor(irf_event_loader_test_config) + epp = DL2EventPreprocessor(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"], @@ -87,9 +87,10 @@ def test_cut_optimizer( proton_full_reco_file, irf_event_loader_test_config, ): - from ctapipe.irf import EventLoader, OptimizationResult, Spectra + from ctapipe.io import DL2EventLoader + from ctapipe.irf import OptimizationResult, Spectra - gamma_loader = EventLoader( + gamma_loader = DL2EventLoader( config=irf_event_loader_test_config, file=gamma_diffuse_full_reco_file, target_spectrum=Spectra.CRAB_HEGRA, @@ -98,7 +99,7 @@ def test_cut_optimizer( chunk_size=10000, obs_time=u.Quantity(50, u.h), ) - proton_loader = EventLoader( + proton_loader = DL2EventLoader( config=irf_event_loader_test_config, file=proton_full_reco_file, target_spectrum=Spectra.IRFDOC_PROTON_SPECTRUM, diff --git a/src/ctapipe/irf/tests/test_preprocessing.py b/src/ctapipe/irf/tests/test_preprocessing.py deleted file mode 100644 index 71cedc22ade..00000000000 --- a/src/ctapipe/irf/tests/test_preprocessing.py +++ /dev/null @@ -1,106 +0,0 @@ -import astropy.units as u -import numpy as np -import pytest -from astropy.table import Table -from pyirf.simulations import SimulatedEventsInfo -from pyirf.spectral import PowerLaw - - -@pytest.fixture(scope="module") -def dummy_table(): - """Dummy table to test column renaming.""" - return Table( - { - "obs_id": [1, 1, 1, 2, 3, 3], - "event_id": [1, 2, 3, 1, 1, 2], - "true_energy": [0.99, 10, 0.37, 2.1, 73.4, 1] * u.TeV, - "dummy_energy": [1, 10, 0.4, 2.5, 73, 1] * u.TeV, - "classifier_prediction": [1, 0.3, 0.87, 0.93, 0, 0.1], - "true_alt": [60, 60, 60, 60, 60, 60] * u.deg, - "geom_alt": [58.5, 61.2, 59, 71.6, 60, 62] * u.deg, - "true_az": [13, 13, 13, 13, 13, 13] * u.deg, - "geom_az": [12.5, 13, 11.8, 15.1, 14.7, 12.8] * u.deg, - "subarray_pointing_frame": np.zeros(6), - "subarray_pointing_lat": np.full(6, 20) * u.deg, - "subarray_pointing_lon": np.full(6, 0) * u.deg, - } - ) - - -def test_normalise_column_names(dummy_table): - from ctapipe.irf import EventPreprocessor - - epp = EventPreprocessor( - energy_reconstructor="dummy", - geometry_reconstructor="geom", - gammaness_classifier="classifier", - ) - norm_table = epp.normalise_column_names(dummy_table) - - needed_cols = [ - "obs_id", - "event_id", - "true_energy", - "true_alt", - "true_az", - "reco_energy", - "reco_alt", - "reco_az", - "gh_score", - "pointing_alt", - "pointing_az", - ] - for c in needed_cols: - assert c in norm_table.colnames - - with pytest.raises(ValueError, match="Required column geom_alt is missing."): - dummy_table.rename_column("geom_alt", "alt_geom") - epp = EventPreprocessor( - energy_reconstructor="dummy", - geometry_reconstructor="geom", - gammaness_classifier="classifier", - ) - _ = epp.normalise_column_names(dummy_table) - - -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, - ) - events, count, meta = loader.load_preselected_events( - chunk_size=10000, - obs_time=u.Quantity(50, u.h), - ) - - columns = [ - "obs_id", - "event_id", - "true_energy", - "true_az", - "true_alt", - "reco_energy", - "reco_az", - "reco_alt", - "reco_fov_lat", - "reco_fov_lon", - "gh_score", - "pointing_az", - "pointing_alt", - "theta", - "true_source_fov_offset", - "reco_source_fov_offset", - ] - assert columns.sort() == events.colnames.sort() - - assert isinstance(count, int) - assert isinstance(meta["sim_info"], SimulatedEventsInfo) - assert isinstance(meta["spectrum"], PowerLaw) - - events = loader.make_event_weights( - events, meta["spectrum"], "gammas", (0 * u.deg, 1 * u.deg) - ) - assert "weight" in events.colnames diff --git a/src/ctapipe/resources/compute_irf.yaml b/src/ctapipe/resources/compute_irf.yaml index 806f73d160b..411e1782440 100644 --- a/src/ctapipe/resources/compute_irf.yaml +++ b/src/ctapipe/resources/compute_irf.yaml @@ -18,12 +18,12 @@ IrfTool: energy_bias_resolution_maker_name: "EnergyBiasResolution2dMaker" sensitivity_maker_name: "Sensitivity2dMaker" -EventPreprocessor: +DL2EventPreprocessor: energy_reconstructor: "RandomForestRegressor" geometry_reconstructor: "HillasReconstructor" gammaness_classifier: "RandomForestClassifier" - EventQualityQuery: + DL2EventQualityQuery: 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..e99ed0fe1df 100644 --- a/src/ctapipe/resources/optimize_cuts.yaml +++ b/src/ctapipe/resources/optimize_cuts.yaml @@ -12,12 +12,12 @@ EventSelectionOptimizer: obs_time: 50 hour optimization_algorithm: "PointSourceSensitivityOptimizer" # Alternative: "PercentileCuts" -EventPreprocessor: +DL2EventPreprocessor: energy_reconstructor: "RandomForestRegressor" geometry_reconstructor: "HillasReconstructor" gammaness_classifier: "RandomForestClassifier" - EventQualityQuery: + DL2EventQualityQuery: quality_criteria: - ["multiplicity 4", "np.count_nonzero(HillasReconstructor_telescopes,axis=1) >= 4"] - ["valid classifier", "RandomForestClassifier_is_valid"] diff --git a/src/ctapipe/tools/compute_irf.py b/src/ctapipe/tools/compute_irf.py index 3c5a06641d3..dd47f8755c5 100644 --- a/src/ctapipe/tools/compute_irf.py +++ b/src/ctapipe/tools/compute_irf.py @@ -19,8 +19,11 @@ from ..core import Provenance, Tool, ToolConfigurationError, traits from ..core.traits import AstroQuantity, Bool, Integer, classes_with_traits, flag +from ..io.dl2_tables_preprocessing import ( + DL2EventLoader, + DL2EventQualityQuery, +) from ..irf import ( - EventLoader, OptimizationResult, Spectra, check_bins_in_range, @@ -36,7 +39,6 @@ EnergyDispersionMakerBase, PSFMakerBase, ) -from ..irf.preprocessing import EventQualityQuery __all__ = ["IrfTool"] @@ -214,7 +216,7 @@ class IrfTool(Tool): classes = ( [ - EventLoader, + DL2EventLoader, ] + classes_with_traits(BackgroundRateMakerBase) + classes_with_traits(EffectiveAreaMakerBase) @@ -267,7 +269,7 @@ def setup(self): raise_error=self.range_check_error, ) self.event_loaders = { - "gammas": EventLoader( + "gammas": DL2EventLoader( parent=self, file=self.gamma_file, target_spectrum=self.gamma_target_spectrum, @@ -281,13 +283,13 @@ def setup(self): "At least a proton file required when specifying `do_background`." ) - self.event_loaders["protons"] = EventLoader( + self.event_loaders["protons"] = DL2EventLoader( parent=self, file=self.proton_file, target_spectrum=self.proton_target_spectrum, ) if self.electron_file and self.electron_file.exists(): - self.event_loaders["electrons"] = EventLoader( + self.event_loaders["electrons"] = DL2EventLoader( parent=self, file=self.electron_file, target_spectrum=self.electron_target_spectrum, @@ -507,7 +509,7 @@ def start(self): ], ) ) - loader.epp.quality_query = EventQualityQuery( + loader.epp.quality_query = DL2EventQualityQuery( parent=loader, quality_criteria=self.opt_result.quality_query.quality_criteria, ) diff --git a/src/ctapipe/tools/optimize_event_selection.py b/src/ctapipe/tools/optimize_event_selection.py index 186b8a34fb4..0cfc57435f8 100644 --- a/src/ctapipe/tools/optimize_event_selection.py +++ b/src/ctapipe/tools/optimize_event_selection.py @@ -5,7 +5,8 @@ from ..core import Provenance, Tool, traits from ..core.traits import AstroQuantity, Integer, classes_with_traits -from ..irf import EventLoader, Spectra +from ..io import DL2EventLoader +from ..irf import Spectra from ..irf.optimize import CutOptimizerBase __all__ = ["EventSelectionOptimizer"] @@ -102,7 +103,7 @@ class EventSelectionOptimizer(Tool): "chunk_size": "EventSelectionOptimizer.chunk_size", } - classes = [EventLoader] + classes_with_traits(CutOptimizerBase) + classes = [DL2EventLoader] + classes_with_traits(CutOptimizerBase) def setup(self): """ @@ -112,7 +113,7 @@ def setup(self): self.optimization_algorithm, parent=self ) self.event_loaders = { - "gammas": EventLoader( + "gammas": DL2EventLoader( parent=self, file=self.gamma_file, target_spectrum=self.gamma_target_spectrum, @@ -127,13 +128,13 @@ def setup(self): f"using {self.optimization_algorithm}." ) - self.event_loaders["protons"] = EventLoader( + self.event_loaders["protons"] = DL2EventLoader( parent=self, file=self.proton_file, target_spectrum=self.proton_target_spectrum, ) if self.electron_file and self.electron_file.exists(): - self.event_loaders["electrons"] = EventLoader( + self.event_loaders["electrons"] = DL2EventLoader( parent=self, file=self.electron_file, target_spectrum=self.electron_target_spectrum, diff --git a/src/ctapipe/tools/tests/test_compute_irf.py b/src/ctapipe/tools/tests/test_compute_irf.py index acd35a50aaf..8a272bcaad3 100644 --- a/src/ctapipe/tools/tests/test_compute_irf.py +++ b/src/ctapipe/tools/tests/test_compute_irf.py @@ -240,11 +240,11 @@ def test_irf_tool_wrong_cuts( with config_path.open("w") as f: json.dump( { - "EventPreprocessor": { + "DL2EventPreprocessor": { "energy_reconstructor": "ExtraTreesRegressor", "geometry_reconstructor": "HillasReconstructor", "gammaness_classifier": "ExtraTreesClassifier", - "EventQualityQuery": { + "DL2EventQualityQuery": { "quality_criteria": [ # No criteria for minimum event multiplicity ("valid classifier", "ExtraTreesClassifier_is_valid"),