Irf maker and cut optimiser

docs/changes/2473.feature.rst

@@ -0,0 +1,3 @@
+Add a ``ctapipe-make-irf`` tool and a able to produce irfs given a cut-selection file and gamma, proton, and electron DL2 input files.
+
+Add a ``ctapipe-optimize-event-selection`` tool to produce cut-selection files.

docs/conf.py

@@ -118,6 +118,7 @@ def setup(app):
  ("py:class", "t.Type"),
  ("py:class", "t.List"),
  ("py:class", "t.Tuple"),
+ ("py:class", "t.Sequence"),
  ("py:class", "Config"),
  ("py:class", "traitlets.config.configurable.Configurable"),
  ("py:class", "traitlets.traitlets.HasTraits"),
@@ -132,35 +133,34 @@ def setup(app):
  ("py:class", "traitlets.config.application.Application"),
  ("py:class", "traitlets.utils.sentinel.Sentinel"),
  ("py:class", "traitlets.traitlets.ObserveHandler"),
+ ("py:class", "traitlets.traitlets.T"),
+ ("py:class", "traitlets.traitlets.G"),
+ ("py:class", "Sentinel"),
+ ("py:class", "ObserveHandler"),
  ("py:class", "dict[K, V]"),
  ("py:class", "G"),
  ("py:class", "K"),
  ("py:class", "V"),
- ("py:class", "t.Sequence"),
  ("py:class", "StrDict"),
  ("py:class", "ClassesType"),
- ("py:class", "traitlets.traitlets.G"),
+ ("py:class", "re.Pattern"),
+ ("py:class", "re.Pattern[t.Any]"),
+ ("py:class", "astropy.coordinates.baseframe.BaseCoordinateFrame"),
+ ("py:class", "astropy.table.table.Table"),
+ ("py:class", "eventio.simtel.simtelfile.SimTelFile"),
+ ("py:class", "ctapipe.compat.StrEnum"),
+ ("py:class", "ctapipe.compat.StrEnum"),
+ ("py:obj", "traitlets.traitlets.T"),
  ("py:obj", "traitlets.traitlets.G"),
  ("py:obj", "traitlets.traitlets.S"),
- ("py:obj", "traitlets.traitlets.T"),
- ("py:class", "traitlets.traitlets.T"),
- ("py:class", "re.Pattern[t.Any]"),
- ("py:class", "re.Pattern"),
- ("py:class", "Sentinel"),
- ("py:class", "ObserveHandler"),
  ("py:obj", "traitlets.config.boolean_flag"),
  ("py:obj", "traitlets.TraitError"),
  ("py:obj", "-v"), # fix for wrong syntax in a traitlets docstring
+ ("py:obj", "cls"),
+ ("py:obj", "name"),
  ("py:meth", "MetaHasDescriptors.__init__"),
  ("py:meth", "HasTraits.__new__"),
  ("py:meth", "BaseDescriptor.instance_init"),
- ("py:obj", "cls"),
- ("py:obj", "name"),
- ("py:class", "astropy.coordinates.baseframe.BaseCoordinateFrame"),
- ("py:class", "astropy.table.table.Table"),
- ("py:class", "eventio.simtel.simtelfile.SimTelFile"),
- ("py:class", "ctapipe.compat.StrEnum"),
- ("py:class", "ctapipe.compat.StrEnum"),
 ]
 
 # Sphinx gallery config

environment.yml

@@ -25,6 +25,7 @@ dependencies:
  - pypandoc
  - pre-commit
  - psutil
+ - pyirf
  - pytables
  - pytest
  - pytest-cov

pyproject.toml

@@ -40,6 +40,7 @@ dependencies = [
  "numba >=0.56",
  "numpy ~=1.16",
  "psutil",
+ "pyirf >0.10.1",
  "pyyaml >=5.1",
  "requests",
  "scikit-learn !=1.4.0", # 1.4.0 breaks with astropy tables, before and after works
@@ -96,6 +97,8 @@ ctapipe-dump-instrument = "ctapipe.tools.dump_instrument:main"
 ctapipe-display-dl1 = "ctapipe.tools.display_dl1:main"
 ctapipe-process = "ctapipe.tools.process:main"
 ctapipe-merge = "ctapipe.tools.merge:main"
+ctapipe-optimize-event-selection = "ctapipe.tools.optimize_event_selection:main"
+ctapipe-make-irf = "ctapipe.tools.make_irf:main"
 ctapipe-fileinfo = "ctapipe.tools.fileinfo:main"
 ctapipe-quickstart = "ctapipe.tools.quickstart:main"
 ctapipe-train-energy-regressor = "ctapipe.tools.train_energy_regressor:main"

src/ctapipe/conftest.py

@@ -1,6 +1,7 @@
 """
 common pytest fixtures for tests in ctapipe
 """
+
 import shutil
 from copy import deepcopy
 
@@ -9,7 +10,7 @@
 import pytest
 import tables
 from astropy.coordinates import EarthLocation
-from astropy.table import Table
+from astropy.table import QTable, Table, vstack
 from pytest_astropy_header.display import PYTEST_HEADER_MODULES
 
 from ctapipe.core import run_tool
@@ -684,3 +685,108 @@ def dl1_mon_pointing_file(dl1_file, dl1_tmp_path):
  f.remove_node("/configuration/telescope/pointing", recursive=True)
 
  return path
+
+
+@pytest.fixture(scope="session")
+def gamma_diffuse_full_reco_file(
+ gamma_train_clf,
+ particle_classifier_path,
+ model_tmp_path,
+):
+ """
+ Energy reconstruction and geometric origin reconstruction have already been done.
+ """
+ from ctapipe.tools.apply_models import ApplyModels
+
+ output_path = model_tmp_path / "gamma_diffuse_full_reco.dl2.h5"
+ run_tool(
+ ApplyModels(),
+ argv=[
+ f"--input={gamma_train_clf}",
+ f"--output={output_path}",
+ f"--reconstructor={particle_classifier_path}",
+ "--no-dl1-parameters",
+ "--StereoMeanCombiner.weights=konrad",
+ ],
+ raises=True,
+ )
+ return output_path
+
+
+@pytest.fixture(scope="session")
+def proton_full_reco_file(
+ proton_train_clf,
+ particle_classifier_path,
+ model_tmp_path,
+):
+ """
+ Energy reconstruction and geometric origin reconstruction have already been done.
+ """
+ from ctapipe.tools.apply_models import ApplyModels
+
+ output_path = model_tmp_path / "proton_full_reco.dl2.h5"
+ run_tool(
+ ApplyModels(),
+ argv=[
+ f"--input={proton_train_clf}",
+ f"--output={output_path}",
+ f"--reconstructor={particle_classifier_path}",
+ "--no-dl1-parameters",
+ "--StereoMeanCombiner.weights=konrad",
+ ],
+ raises=True,
+ )
+ return output_path
+
+
+@pytest.fixture(scope="session")
+def irf_events_loader_test_config():
+ from traitlets.config import Config
+
+ return Config(
+ {
+ "EventPreProcessor": {
+ "energy_reconstructor": "ExtraTreesRegressor",
+ "geometry_reconstructor": "HillasReconstructor",
+ "gammaness_classifier": "ExtraTreesClassifier",
+ "quality_criteria": [
+ (
+ "multiplicity 4",
+ "np.count_nonzero(tels_with_trigger,axis=1) >= 4",
+ ),
+ ("valid classifier", "ExtraTreesClassifier_is_valid"),
+ ("valid geom reco", "HillasReconstructor_is_valid"),
+ ("valid energy reco", "ExtraTreesRegressor_is_valid"),
+ ],
+ }
+ }
+ )
+
+
+@pytest.fixture(scope="session")
+def irf_events_table():
+ from ctapipe.irf import EventPreProcessor
+
+ N1 = 1000
+ N2 = 100
+ N = N1 + N2
+ epp = EventPreProcessor()
+ tab = epp.make_empty_table()
+ # Add fake weight column
+ tab.add_column((), name="weight")
+ units = {c: tab[c].unit for c in tab.columns}
+
+ empty = np.zeros((len(tab.columns), N)) * np.nan
+ e_tab = QTable(data=empty.T, names=tab.colnames, units=units)
+ # Setting values following pyirf test in pyirf/irf/tests/test_background.py
+ e_tab["reco_energy"] = np.append(np.full(N1, 1), np.full(N2, 2)) * u.TeV
+ e_tab["true_energy"] = np.append(np.full(N1, 0.9), np.full(N2, 2.1)) * u.TeV
+ e_tab["reco_source_fov_offset"] = (
+ np.append(np.full(N1, 0.1), np.full(N2, 0.05)) * u.deg
+ )
+ e_tab["true_source_fov_offset"] = (
+ np.append(np.full(N1, 0.11), np.full(N2, 0.04)) * u.deg
+ )
+
+ ev = vstack([e_tab, tab], join_type="exact", metadata_conflicts="silent")
+ return ev

src/ctapipe/irf/__init__.py

@@ -0,0 +1,51 @@
+"""Top level module for the irf functionality"""
+
+from .benchmarks import (
+ AngularResolution2dMaker,
+ EnergyBiasResolution2dMaker,
+ Sensitivity2dMaker,
+)
+from .binning import (
+ ResultValidRange,
+ check_bins_in_range,
+ make_bins_per_decade,
+)
+from .irfs import (
+ BackgroundRate2dMaker,
+ EffectiveArea2dMaker,
+ EnergyMigration2dMaker,
+ Psf3dMaker,
+)
+from .optimize import (
+ GhPercentileCutCalculator,
+ OptimizationResult,
+ OptimizationResultStore,
+ PercentileCuts,
+ PointSourceSensitivityOptimizer,
+ ThetaPercentileCutCalculator,
+)
+from .select import EventPreProcessor, EventsLoader
+from .spectra import SPECTRA, Spectra
+
+__all__ = [
+ "AngularResolution2dMaker",
+ "EnergyBiasResolution2dMaker",
+ "Sensitivity2dMaker",
+ "Psf3dMaker",
+ "BackgroundRate2dMaker",
+ "EnergyMigration2dMaker",
+ "EffectiveArea2dMaker",
+ "ResultValidRange",
+ "OptimizationResult",
+ "OptimizationResultStore",
+ "PointSourceSensitivityOptimizer",
+ "PercentileCuts",
+ "EventsLoader",
+ "EventPreProcessor",
+ "Spectra",
+ "GhPercentileCutCalculator",
+ "ThetaPercentileCutCalculator",
+ "SPECTRA",
+ "check_bins_in_range",
+ "make_bins_per_decade",
+]