From 159bff3b2d406b72207709f3bb1122c056e62260 Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Tue, 27 Jan 2026 09:23:28 +0100 Subject: [PATCH 01/16] add EventWeigter classes and test --- src/ctapipe/irf/__init__.py | 10 + src/ctapipe/irf/event_weighter.py | 185 +++++++++++++++++++ src/ctapipe/irf/tests/test_event_weighter.py | 135 ++++++++++++++ 3 files changed, 330 insertions(+) create mode 100644 src/ctapipe/irf/event_weighter.py create mode 100644 src/ctapipe/irf/tests/test_event_weighter.py diff --git a/src/ctapipe/irf/__init__.py b/src/ctapipe/irf/__init__.py index d4abdb43fac..db1852a03ae 100644 --- a/src/ctapipe/irf/__init__.py +++ b/src/ctapipe/irf/__init__.py @@ -18,6 +18,12 @@ check_bins_in_range, make_bins_per_decade, ) +from .event_weighter import ( + EventWeighter, + PolarEventWeighter, + RadialEventWeighter, + SimpleEventWeighter, +) from .irfs import ( BackgroundRate2dMaker, EffectiveArea2dMaker, @@ -53,4 +59,8 @@ "FLUX_UNIT", "check_bins_in_range", "make_bins_per_decade", + "EventWeighter", + "PolarEventWeighter", + "RadialEventWeighter", + "SimpleEventWeighter", ] diff --git a/src/ctapipe/irf/event_weighter.py b/src/ctapipe/irf/event_weighter.py new file mode 100644 index 00000000000..468485ea33b --- /dev/null +++ b/src/ctapipe/irf/event_weighter.py @@ -0,0 +1,185 @@ +#!/usr/bin/env python3 + +"""Components for computing event weights.""" + +from abc import abstractmethod +from collections.abc import Callable +from typing import override + +import numpy as np +from astropy import units as u +from astropy.table import QTable +from pyirf.spectral import ( + calculate_event_weights, +) + +from ..core import Component, traits +from ..core.feature_generator import _shallow_copy_table +from .binning import DefaultFoVOffsetBins, DefaultFoVPhiBins +from .spectra import SPECTRA, Spectra + +__all__ = [ + "EventWeighter", + "SimpleEventWeighter", + "RadialEventWeighter", + "PolarEventWeighter", +] + + +class EventWeighter(Component): + """Compute weights to go from source to target spectra.""" + + target_spectrum_name = traits.UseEnum( + Spectra, + default_value=Spectra.CRAB_HEGRA, + help="Pre-defined source spectrum to reweight to.", + ).tag(config=True) + + is_diffuse = traits.Bool( + default_value=True, help="If True, assume the source is diffuse." + ).tag(config=True) + + energy_column = traits.Unicode( + help="name of energy column", default_value="true_energy" + ).tag(config=True) + weight_column = traits.Unicode( + help="name of output weight column", default_value="weight" + ).tag(config=True) + + def __init__( + self, + source_spectrum: Callable[[u.Quantity], u.Quantity], + config=None, + parent=None, + **kwargs, + ): + """ + Parameters + ---------- + source_spectrum: Callable + initial spectrum of the events to be processed. + target_spectrum: Callable | None + target spectrum to weight to. If None, a pre-defined spectrum + function from the `target_spectrum_name` attribute will be used` + """ + super().__init__(config=config, parent=parent, **kwargs) + self.source_spectrum = source_spectrum + self.target_spectrum = SPECTRA[self.target_spectrum_name] + + @abstractmethod + def _compute_weights(self, events_table: QTable): + raise NotImplementedError( + f"{self.__class__.__name__} weighting is not implemented" + ) + + def __call__(self, events_table: QTable) -> QTable: + """Returns shallow copy of input table with a `weight` column added""" + + table = _shallow_copy_table(events_table) + self._compute_weights(table) + return table + + +class SimpleEventWeighter(EventWeighter): + """Weights all events spectrally with no spatial binning. + + Calling this class adds a column to the output table with the event-wise + spectral weights, with column name ``weight_column``. + """ + + fov_offset_max = traits.AstroQuantity( + help="upper bound of spatial integral applied to source function", + default_value=u.Quantity(10, u.deg), + physical_type=u.physical.angle, + ).tag(config=True) + + @override + def _compute_weights(self, events_table: QTable): + energy = events_table[self.energy_column] + source_spectrum = self.source_spectrum + if self.is_diffuse: + source_spectrum = source_spectrum.integrate_cone( + 0 * u.deg, self.fov_offset_max + ) + weights = calculate_event_weights( + energy, + target_spectrum=self.target_spectrum, + simulated_spectrum=source_spectrum, + ) + events_table[self.weight_column] = weights + + +class RadialEventWeighter(EventWeighter, DefaultFoVOffsetBins): + """ + Weights in radial (FOV) offset bins in the `~ctapipe.coordinates.NominalFrame`. + + Calling this class adds a column to the output table with the event-wise + spectral-spatial weights, with column name `weight_column``. This implementation + additionally adds the column ``output_table["fov_offset_bin"]``, and the + list of offset bin edges in ``output_table.meta["OFFSBINS"]`` + """ + + fov_offset_column = traits.Unicode( + help="name of FOV radial offset column", default_value="true_fov_offset" + ).tag(config=True) + + @override + def _compute_weights(self, events_table: QTable): + offset_bins = self.fov_offset_bins + offset = events_table[self.fov_offset_column].to_value(offset_bins.unit) + energy = events_table[self.energy_column] + weights = np.zeros_like(energy.value) + + # note that the bin i from digitize starts at 1 and means: + # offset_bins[i-1] <= offset < offset_bins[ii]) + r_bin = np.digitize(offset, offset_bins.value) + + for ii in range(0, len(offset_bins)): + self.log.debug( + f"bin {ii} offset=[{offset_bins[ii - 1]}, {offset_bins[ii]})" + ) + mask = r_bin == ii + weights[mask] = calculate_event_weights( + true_energy=energy[mask], + target_spectrum=self.target_spectrum, + simulated_spectrum=self.source_spectrum.integrate_cone( + offset_bins[ii - 1], offset_bins[ii] + ), + ) + + events_table[self.weight_column] = weights + events_table["fov_offset_bin"] = r_bin + + events_table.columns["fov_offset_bin"].description = ( + "Bin i defined as offset[i-1] <= fov_offset < offset[i]. " + "Where offset is `OFFSBINS` array found this table's metadata." + ) + + events_table.meta["OFFSBINS"] = list(offset_bins.to_value("deg")) + + +class PolarEventWeighter(EventWeighter, DefaultFoVOffsetBins, DefaultFoVPhiBins): + """ + Weights in field-of-view polar $(r, \\phi)$ bins in the `~ctapipe.coordinates.NominalFrame`. + + Calling this class adds a column to the output table with the event-wise + spectral-spatial weights, with column name ``weight_column``. This + implementation additionally adds the columns + ``output_table["fov_offset_bin"]``, ``output_table["fov_phi_bin"]``, and the + list of offset and phi bin edges in ``output_table.meta["OFFSBINS"]`` and + ``output_table.meta["PHIBINS"]`` respectively. + """ + + fov_offset_column = traits.Unicode( + help="name of FOV radial offset column", default_value="true_fov_offset" + ).tag(config=True) + + fov_phi_column = traits.Unicode( + help="name of the FOV azimuthal coordinate", default_value="true_fov_phi" + ).tag(config=True) + + @override + def _compute_weights(self, events_table: QTable): + raise NotImplementedError( + f"{self.__class__.__name__} weighting is not implemented" + ) diff --git a/src/ctapipe/irf/tests/test_event_weighter.py b/src/ctapipe/irf/tests/test_event_weighter.py new file mode 100644 index 00000000000..c903d59c44c --- /dev/null +++ b/src/ctapipe/irf/tests/test_event_weighter.py @@ -0,0 +1,135 @@ +#!/usr/bin/env python3 + +import numpy as np +import pytest +from astropy import units as u +from astropy.table import QTable + +from ctapipe.irf.event_weighter import RadialEventWeighter, SimpleEventWeighter + + +@pytest.fixture(scope="function") +def example_weight_table(): + table = QTable( + dict( + true_energy=[1.0, 2.0, 0.5, 0.2] * u.TeV, + true_fov_offset=[0.1, 1.2, 2.2, 3.2] * u.deg, + ) + ) + table.meta["VERSION"] = 1.0 + table.columns["true_energy"].description = "True energy of particle" + return table + + +def test_simple_weights(example_weight_table): + from ctapipe.irf.spectra import PowerLaw + + table = example_weight_table + + source_spectrum = PowerLaw( + normalization=u.Quantity(0.00027, "TeV-1 s-1 sr-1 m-2"), + index=-2.0, + e_ref=1.0 * u.TeV, + ) + + weight1 = SimpleEventWeighter( + source_spectrum=source_spectrum, target_spectrum_name="CRAB_HEGRA" + ) + + table_w1 = weight1(table) + + assert np.all(table_w1["weight"] > 0.0) + assert np.all(table_w1["weight"] <= 1.0) + + +def test_flat_weighting(example_weight_table): + """Check that if source and target spectra are the same, we get 1.0.""" + from ctapipe.irf.spectra import SPECTRA, Spectra + + table = example_weight_table + + weight = SimpleEventWeighter( + source_spectrum=SPECTRA[Spectra.CRAB_HEGRA], + target_spectrum_name=Spectra.CRAB_HEGRA.name, + is_diffuse=False, + ) + + w = weight(table)["weight"] + + assert np.allclose(w, 1.0) + + +def test_radial_weights(example_weight_table): + from ctapipe.irf.spectra import PowerLaw + + table = example_weight_table + + source_spectrum = PowerLaw( + normalization=u.Quantity(0.00027, "TeV-1 s-1 sr-1 m-2"), + index=-2.0, + e_ref=1.0 * u.TeV, + ) + + weight = RadialEventWeighter( + source_spectrum=source_spectrum, + target_spectrum_name="CRAB_HEGRA", + fov_offset_min=0.0 * u.deg, + fov_offset_max=5.0 * u.deg, + fov_offset_n_bins=5, + ) + + assert np.allclose(weight.fov_offset_bins, [0.0, 1.0, 2.0, 3.0, 4.0, 5.0] * u.deg) + + table_w = weight(table) + + assert "fov_offset_bin" in table_w.colnames + assert "OFFSBINS" in table_w.meta + assert table_w.meta["OFFSBINS"] == [0.0, 1.0, 2.0, 3.0, 4.0, 5.0] + + for ii in range(0, 7): + mask = table_w["fov_offset_bin"] == ii + + if ii == 0: + # there should be no 0th bin + assert len(table_w[mask]) == 0 + elif ii == 6: + # for outlier bin 6, weights should be all 0 + assert np.all(table_w["weight"][mask] == 0.0) + else: + # for bins 1-5,weights should be positive in this case + assert np.all(table_w["weight"][mask] >= 0.0) + + +def test_radial_weights_fits(example_weight_table, tmp_path): + """Test round-tripping to FITS""" + from astropy.table import Table + + from ctapipe.irf.spectra import PowerLaw + + table = example_weight_table + + source_spectrum = PowerLaw( + normalization=u.Quantity(0.00027, "TeV-1 s-1 sr-1 m-2"), + index=-2.0, + e_ref=1.0 * u.TeV, + ) + + weight = RadialEventWeighter( + source_spectrum=source_spectrum, + target_spectrum_name="CRAB_HEGRA", + fov_offset_min=0.0 * u.deg, + fov_offset_max=5.0 * u.deg, + fov_offset_n_bins=5, + ) + + table_w = weight(table) + table_w.write(tmp_path / "test.fits") + + assert Table.read(tmp_path / "test.fits").meta["OFFSBINS"] == [ + 0.0, + 1.0, + 2.0, + 3.0, + 4.0, + 5.0, + ] From 8754c5a1ef9dd9e31d2fbbca29a4bb6b62b73821 Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Tue, 27 Jan 2026 09:24:54 +0100 Subject: [PATCH 02/16] make fov_offset_bins a property --- src/ctapipe/irf/binning.py | 31 ++++++++++++++++++++++++++++--- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/src/ctapipe/irf/binning.py b/src/ctapipe/irf/binning.py index 172e159b3d2..fda7fb12653 100644 --- a/src/ctapipe/irf/binning.py +++ b/src/ctapipe/irf/binning.py @@ -17,6 +17,7 @@ "DefaultTrueEnergyBins", "DefaultRecoEnergyBins", "DefaultFoVOffsetBins", + "DefaultFoVPhiBins", ] logger = logging.getLogger(__name__) @@ -179,9 +180,9 @@ class DefaultFoVOffsetBins(Component): default_value=1, ).tag(config=True) - def __init__(self, config=None, parent=None, **kwargs): - super().__init__(config=config, parent=parent, **kwargs) - self.fov_offset_bins = u.Quantity( + @property + def fov_offset_bins(self): + return u.Quantity( np.linspace( self.fov_offset_min.to_value(u.deg), self.fov_offset_max.to_value(u.deg), @@ -189,3 +190,27 @@ def __init__(self, config=None, parent=None, **kwargs): ), u.deg, ) + + +class DefaultFoVPhiBins(Component): + """ + Base class for creating IRFs with phi dependence. + + The range is always assumed to be (0, 360) deg. + """ + + fov_phi_n_bins = Integer( + help="Number of FoV offset bins", + default_value=4, + ).tag(config=True) + + @property + def fov_phi_bins(self): + return u.Quantity( + np.linspace( + 0.0 * u.deg, + 360.0 * u.deg, + self.fov_phi_n_bins + 1, + ), + u.deg, + ) From ce9b7ca1cdf9b6bf7baac5672ea50caa0bb229be Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Tue, 27 Jan 2026 09:26:26 +0100 Subject: [PATCH 03/16] add spectrum_from_simulation_config() helper --- src/ctapipe/irf/__init__.py | 9 ++- src/ctapipe/irf/spectra.py | 114 +++++++++++++++++++++++++++++++++--- 2 files changed, 115 insertions(+), 8 deletions(-) diff --git a/src/ctapipe/irf/__init__.py b/src/ctapipe/irf/__init__.py index db1852a03ae..757861495a8 100644 --- a/src/ctapipe/irf/__init__.py +++ b/src/ctapipe/irf/__init__.py @@ -37,7 +37,13 @@ PointSourceSensitivityOptimizer, ThetaPercentileCutCalculator, ) -from .spectra import ENERGY_FLUX_UNIT, FLUX_UNIT, SPECTRA, Spectra +from .spectra import ( + ENERGY_FLUX_UNIT, + FLUX_UNIT, + SPECTRA, + Spectra, + spectrum_from_simulation_config, +) __all__ = [ "AngularResolution2dMaker", @@ -63,4 +69,5 @@ "PolarEventWeighter", "RadialEventWeighter", "SimpleEventWeighter", + "spectrum_from_simulation_config", ] diff --git a/src/ctapipe/irf/spectra.py b/src/ctapipe/irf/spectra.py index 75106112b97..0e1b8d380f3 100644 --- a/src/ctapipe/irf/spectra.py +++ b/src/ctapipe/irf/spectra.py @@ -1,22 +1,38 @@ """Definition of spectra to be used to calculate event weights for irf computation""" -from enum import Enum +import logging +from collections.abc import Callable +from enum import StrEnum import astropy.units as u -from pyirf.spectral import CRAB_HEGRA, IRFDOC_ELECTRON_SPECTRUM, IRFDOC_PROTON_SPECTRUM +import numpy as np +from astropy.table import Table +from pyirf.simulations import SimulatedEventsInfo +from pyirf.spectral import ( + CRAB_HEGRA, + IRFDOC_ELECTRON_SPECTRUM, + IRFDOC_PROTON_SPECTRUM, + PowerLaw, +) -__all__ = ["ENERGY_FLUX_UNIT", "FLUX_UNIT", "SPECTRA", "Spectra"] +__all__ = [ + "ENERGY_FLUX_UNIT", + "FLUX_UNIT", + "SPECTRA", + "Spectra", + "spectrum_from_simulation_config", +] ENERGY_FLUX_UNIT = (1 * u.erg / u.s / u.cm**2).unit FLUX_UNIT = (1 / u.erg / u.s / u.cm**2).unit -class Spectra(Enum): +class Spectra(StrEnum): """Spectra for calculating event weights""" - CRAB_HEGRA = 1 - IRFDOC_ELECTRON_SPECTRUM = 2 - IRFDOC_PROTON_SPECTRUM = 3 + CRAB_HEGRA = "CRAB_HEGRA" + IRFDOC_ELECTRON_SPECTRUM = "IRFDOC_ELECTRON_SPECTRUM" + IRFDOC_PROTON_SPECTRUM = "IRFDOC_PROTON_SPECTRUM" SPECTRA = { @@ -24,3 +40,87 @@ class Spectra(Enum): Spectra.IRFDOC_ELECTRON_SPECTRUM: IRFDOC_ELECTRON_SPECTRUM, Spectra.IRFDOC_PROTON_SPECTRUM: IRFDOC_PROTON_SPECTRUM, } + + +logger = logging.getLogger(__name__) + + +def spectrum_from_simulation_config( + simulation_configuration_table: Table, + shower_distribution_table: Table | None, + obs_time: u.Quantity, + method: str = "powerlaw", +) -> Callable[[u.Quantity], u.Quantity]: + """ + Return simulated spectrum function from configuration information.. + + Note this currently implements only PowerLaw spectra, but in the future will returnq a + + Parameters + ---------- + simulation_configuration_table: Table + table as read by TableLoader.read_simulation_configuration() + shower_distribution_table: Table + table of simulated shower distribution as read from TableLoader.read_shower_distribution() + obs_time: u.Quantity['s'] + Observation time in a unit convertible to seconds. + method: str + "powerlaw" (for assuming powerlaw distributions), or "histogram" to return an + interpolation function from the underlying distribution histogram. + + Returns + ------- + Callable: + simulated spectrum function. + """ + + if method != "powerlaw": + return NotImplementedError(f"Method '{method}' is not implemented") + + n_showers: int = 0 + if shower_distribution_table: + n_showers = shower_distribution_table["n_entries"].sum() + else: + logger.warning( + "Simulation distributions were not found in the input files, " + "falling back to estimating the number of showers from the " + "simulation configuration." + ) + + # Some tight consistency checks. Eventually we will support using the + # arbitrary shower distribution and non-flat spatial distributions. + # Currently we do not support those, so we raise exceptions here to + # avoid that we incorrectly compute the effective area, which would have + # a high scientific impact. + for itm in [ + "spectral_index", + "energy_range_min", + "energy_range_max", + "max_scatter_range", + "max_viewcone_radius", + "min_viewcone_radius", + ]: + if len(np.unique(simulation_configuration_table[itm])) > 1: + raise NotImplementedError( + f"Unsupported: '{itm}' differs across simulation runs" + ) + + n_showers_config = ( + simulation_configuration_table["n_showers"] + * simulation_configuration_table["shower_reuse"] + ).sum() + if n_showers == 0: + n_showers = n_showers_config + + assert n_showers_config == n_showers, "Inconsistent number of simulated showers" + sim_info = SimulatedEventsInfo( + n_showers=n_showers, + energy_min=simulation_configuration_table["energy_range_min"].quantity[0], + energy_max=simulation_configuration_table["energy_range_max"].quantity[0], + max_impact=simulation_configuration_table["max_scatter_range"].quantity[0], + spectral_index=simulation_configuration_table["spectral_index"][0], + viewcone_max=simulation_configuration_table["max_viewcone_radius"].quantity[0], + viewcone_min=simulation_configuration_table["min_viewcone_radius"].quantity[0], + ) + + return PowerLaw.from_simulation(sim_info, obstime=obs_time) From cbf38f9cebb4e67092c9c50e5dc6464c287dafb4 Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Tue, 27 Jan 2026 09:58:12 +0100 Subject: [PATCH 04/16] add spectrum_from_name() helper --- src/ctapipe/irf/__init__.py | 2 ++ src/ctapipe/irf/event_weighter.py | 4 ++-- src/ctapipe/irf/spectra.py | 9 +++++++++ src/ctapipe/irf/tests/test_event_weighter.py | 4 ++-- 4 files changed, 15 insertions(+), 4 deletions(-) diff --git a/src/ctapipe/irf/__init__.py b/src/ctapipe/irf/__init__.py index 757861495a8..6c851611238 100644 --- a/src/ctapipe/irf/__init__.py +++ b/src/ctapipe/irf/__init__.py @@ -42,6 +42,7 @@ FLUX_UNIT, SPECTRA, Spectra, + spectrum_from_name, spectrum_from_simulation_config, ) @@ -70,4 +71,5 @@ "RadialEventWeighter", "SimpleEventWeighter", "spectrum_from_simulation_config", + "spectrum_from_name", ] diff --git a/src/ctapipe/irf/event_weighter.py b/src/ctapipe/irf/event_weighter.py index 468485ea33b..06e94f2c2f4 100644 --- a/src/ctapipe/irf/event_weighter.py +++ b/src/ctapipe/irf/event_weighter.py @@ -16,7 +16,7 @@ from ..core import Component, traits from ..core.feature_generator import _shallow_copy_table from .binning import DefaultFoVOffsetBins, DefaultFoVPhiBins -from .spectra import SPECTRA, Spectra +from .spectra import Spectra, spectrum_from_name __all__ = [ "EventWeighter", @@ -64,7 +64,7 @@ def __init__( """ super().__init__(config=config, parent=parent, **kwargs) self.source_spectrum = source_spectrum - self.target_spectrum = SPECTRA[self.target_spectrum_name] + self.target_spectrum = spectrum_from_name(self.target_spectrum_name) @abstractmethod def _compute_weights(self, events_table: QTable): diff --git a/src/ctapipe/irf/spectra.py b/src/ctapipe/irf/spectra.py index 0e1b8d380f3..496caa129f0 100644 --- a/src/ctapipe/irf/spectra.py +++ b/src/ctapipe/irf/spectra.py @@ -21,6 +21,7 @@ "SPECTRA", "Spectra", "spectrum_from_simulation_config", + "spectrum_from_name", ] ENERGY_FLUX_UNIT = (1 * u.erg / u.s / u.cm**2).unit @@ -42,6 +43,14 @@ class Spectra(StrEnum): } +def spectrum_from_name(name: str | Spectra) -> Callable: + """ + Return spectrum function given a name or Spectra enum value. + """ + spectrum_value = Spectra(name) + return SPECTRA[spectrum_value] + + logger = logging.getLogger(__name__) diff --git a/src/ctapipe/irf/tests/test_event_weighter.py b/src/ctapipe/irf/tests/test_event_weighter.py index c903d59c44c..3e4bcb3746e 100644 --- a/src/ctapipe/irf/tests/test_event_weighter.py +++ b/src/ctapipe/irf/tests/test_event_weighter.py @@ -44,12 +44,12 @@ def test_simple_weights(example_weight_table): def test_flat_weighting(example_weight_table): """Check that if source and target spectra are the same, we get 1.0.""" - from ctapipe.irf.spectra import SPECTRA, Spectra + from ctapipe.irf.spectra import Spectra, spectrum_from_name table = example_weight_table weight = SimpleEventWeighter( - source_spectrum=SPECTRA[Spectra.CRAB_HEGRA], + source_spectrum=spectrum_from_name(Spectra.CRAB_HEGRA), target_spectrum_name=Spectra.CRAB_HEGRA.name, is_diffuse=False, ) From 9d1fc22c3980f911c4136f9d74544944637340d2 Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Tue, 27 Jan 2026 10:04:23 +0100 Subject: [PATCH 05/16] add changelog --- docs/changes/2927.feature.rst | 40 +++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 docs/changes/2927.feature.rst diff --git a/docs/changes/2927.feature.rst b/docs/changes/2927.feature.rst new file mode 100644 index 00000000000..6e65960c4bd --- /dev/null +++ b/docs/changes/2927.feature.rst @@ -0,0 +1,40 @@ +Introduces a new Component `~ctapipe.irf.EventWeighter` that has different implementations of spectral event weighting: + +* `SimpleEventWeigher`: spectral weighting for the full FOV +* `RadialEventWeighter`: spectral weighing in radial bins in the FOV +* `PolarEventWeighter`: spectral weighting in radial+azimuthal (r,phi) bins in the FOV + +They operate on a pre-processed table of DL2 information, where you specify the energy and FOV coordinate columns to use. The interface is as follows: + +.. code-block:: python + + from astropy.table import QTable + + from ctapipe.irf import RadialEventWeighter spectrum_from_name + + table = QTable( + dict( + true_energy=[1.0, 2.0, 0.5, 0.2] * u.TeV, + true_fov_offset=[0.1, 1.2, 2.2, 3.2] * u.deg, + ) + ) + weighter = RadialEventWeighter( + source_spectrum=spectrum_from_name("IRFDOC_ELECTRON_SPECTRUM"), + target_spectrum_name="CRAB_HEGRA", + fov_offset_max=5.0*u.deg, + fov_offset_n_bins=5, + ) + + table_with_weights = weighter(table) + print(table_with_weights) + + +:: + + true_energy true_fov_offset weight fov_offset_bin + TeV deg + ----------- --------------- ------------------ -------------- + 1.0 0.1 12.399508446433442 1 + 2.0 1.2 7.247055871572714 2 + 0.5 2.2 1.4149219056909543 3 + 0.2 3.2 0.4812883464910382 4 From 08670a5e78ac7e589e7a6f5c9a8cdeb326602cfc Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Tue, 27 Jan 2026 12:02:06 +0100 Subject: [PATCH 06/16] docstring fixes --- src/ctapipe/irf/event_weighter.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/ctapipe/irf/event_weighter.py b/src/ctapipe/irf/event_weighter.py index 06e94f2c2f4..6d309b62038 100644 --- a/src/ctapipe/irf/event_weighter.py +++ b/src/ctapipe/irf/event_weighter.py @@ -58,9 +58,6 @@ def __init__( ---------- source_spectrum: Callable initial spectrum of the events to be processed. - target_spectrum: Callable | None - target spectrum to weight to. If None, a pre-defined spectrum - function from the `target_spectrum_name` attribute will be used` """ super().__init__(config=config, parent=parent, **kwargs) self.source_spectrum = source_spectrum @@ -73,7 +70,7 @@ def _compute_weights(self, events_table: QTable): ) def __call__(self, events_table: QTable) -> QTable: - """Returns shallow copy of input table with a `weight` column added""" + """Returns shallow copy of input table with a ``weight`` column added""" table = _shallow_copy_table(events_table) self._compute_weights(table) From 5a2f6b5e00ede5ea3235e2bfc43cb70ac99a8baa Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Wed, 28 Jan 2026 10:46:01 +0100 Subject: [PATCH 07/16] added implementation of PolarEventWeighter --- src/ctapipe/irf/event_weighter.py | 47 +++++++++++++++-- src/ctapipe/irf/tests/test_event_weighter.py | 54 +++++++++++++++++++- 2 files changed, 96 insertions(+), 5 deletions(-) diff --git a/src/ctapipe/irf/event_weighter.py b/src/ctapipe/irf/event_weighter.py index 6d309b62038..d40b58fe4ff 100644 --- a/src/ctapipe/irf/event_weighter.py +++ b/src/ctapipe/irf/event_weighter.py @@ -4,6 +4,7 @@ from abc import abstractmethod from collections.abc import Callable +from itertools import product from typing import override import numpy as np @@ -111,7 +112,7 @@ class RadialEventWeighter(EventWeighter, DefaultFoVOffsetBins): Weights in radial (FOV) offset bins in the `~ctapipe.coordinates.NominalFrame`. Calling this class adds a column to the output table with the event-wise - spectral-spatial weights, with column name `weight_column``. This implementation + spectral-spatial weights, with column name ``weight_column``. This implementation additionally adds the column ``output_table["fov_offset_bin"]``, and the list of offset bin edges in ``output_table.meta["OFFSBINS"]`` """ @@ -131,7 +132,7 @@ def _compute_weights(self, events_table: QTable): # offset_bins[i-1] <= offset < offset_bins[ii]) r_bin = np.digitize(offset, offset_bins.value) - for ii in range(0, len(offset_bins)): + for ii in range(1, len(offset_bins)): self.log.debug( f"bin {ii} offset=[{offset_bins[ii - 1]}, {offset_bins[ii]})" ) @@ -177,6 +178,44 @@ class PolarEventWeighter(EventWeighter, DefaultFoVOffsetBins, DefaultFoVPhiBins) @override def _compute_weights(self, events_table: QTable): - raise NotImplementedError( - f"{self.__class__.__name__} weighting is not implemented" + offset_bins = self.fov_offset_bins + phi_bins = self.fov_phi_bins + offset = events_table[self.fov_offset_column].to_value(offset_bins.unit) + phi = events_table[self.fov_phi_column].to_value(phi_bins.unit) + energy = events_table[self.energy_column] + weights = np.zeros_like(energy.value) + + r_bin = np.digitize(offset, offset_bins.value) + phi_bin = np.digitize(phi, phi_bins.value) + + for i_r, i_phi in product(range(1, len(offset_bins)), range(1, len(phi_bins))): + mask = (r_bin == i_r) & (phi_bin == i_phi) + phi_solid_angle_fraction = (phi_bins[i_phi] - phi_bins[i_phi - 1]) / ( + 360.0 * u.deg + ) + print(i_phi, phi_solid_angle_fraction) + weights[mask] = calculate_event_weights( + true_energy=energy[mask], + target_spectrum=self.target_spectrum, + simulated_spectrum=lambda E: self.source_spectrum.integrate_cone( + offset_bins[i_r - 1], offset_bins[i_r] + )(E) + * phi_solid_angle_fraction, + ) + + events_table[self.weight_column] = weights + events_table["fov_offset_bin"] = r_bin + events_table["fov_phi_bin"] = phi_bin + + events_table.columns["fov_offset_bin"].description = ( + "Bin i defined as offset[i-1] <= fov_offset < offset[i]. " + "Where offset is `OFFSBINS` array found this table's metadata." + ) + + events_table.columns["fov_phi_bin"].description = ( + "Bin i defined as phi[i-1] <= fov_phi < phi[i]. " + "Where phi is `PHIBINS` array found this table's metadata." ) + + events_table.meta["OFFSBINS"] = list(offset_bins.to_value("deg")) + events_table.meta["PHIBINS"] = list(phi_bins.to_value("deg")) diff --git a/src/ctapipe/irf/tests/test_event_weighter.py b/src/ctapipe/irf/tests/test_event_weighter.py index 3e4bcb3746e..305cb1c8232 100644 --- a/src/ctapipe/irf/tests/test_event_weighter.py +++ b/src/ctapipe/irf/tests/test_event_weighter.py @@ -5,7 +5,11 @@ from astropy import units as u from astropy.table import QTable -from ctapipe.irf.event_weighter import RadialEventWeighter, SimpleEventWeighter +from ctapipe.irf.event_weighter import ( + PolarEventWeighter, + RadialEventWeighter, + SimpleEventWeighter, +) @pytest.fixture(scope="function") @@ -14,6 +18,7 @@ def example_weight_table(): dict( true_energy=[1.0, 2.0, 0.5, 0.2] * u.TeV, true_fov_offset=[0.1, 1.2, 2.2, 3.2] * u.deg, + true_fov_phi=[45.0, 45.0, 268.0, 120.0] * u.deg, ) ) table.meta["VERSION"] = 1.0 @@ -133,3 +138,50 @@ def test_radial_weights_fits(example_weight_table, tmp_path): 4.0, 5.0, ] + + +def test_polar_weights(example_weight_table): + from ctapipe.irf.spectra import PowerLaw + + table = example_weight_table + + source_spectrum = PowerLaw( + normalization=u.Quantity(0.00027, "TeV-1 s-1 sr-1 m-2"), + index=-2.0, + e_ref=1.0 * u.TeV, + ) + + weight = PolarEventWeighter( + source_spectrum=source_spectrum, + target_spectrum_name="CRAB_HEGRA", + fov_offset_min=0.0 * u.deg, + fov_offset_max=5.0 * u.deg, + fov_offset_n_bins=5, + fov_phi_n_bins=4, + ) + + assert np.allclose(weight.fov_offset_bins, [0.0, 1.0, 2.0, 3.0, 4.0, 5.0] * u.deg) + assert np.allclose(weight.fov_phi_bins, [0.0, 90.0, 180.0, 270.0, 360.0] * u.deg) + + table_w = weight(table) + + assert "fov_offset_bin" in table_w.colnames + assert "fov_phi_bin" in table_w.colnames + assert "OFFSBINS" in table_w.meta + assert "PHIBINS" in table_w.meta + assert table_w.meta["OFFSBINS"] == [0.0, 1.0, 2.0, 3.0, 4.0, 5.0] + assert table_w.meta["PHIBINS"] == [0.0, 90.0, 180.0, 270.0, 360.0] + + for ii in range(0, 7): + for jj in range(0, 4): + mask = (table_w["fov_offset_bin"] == ii) & (table_w["fov_phi_bin"] == jj) + + if ii == 0 or jj == 0: + # there should be no 0th bin + assert len(table_w[mask]) == 0 + elif ii == 6 or jj == 5: + # for outlier bin 6,5, weights should be all 0 + assert np.all(table_w["weight"][mask] == 0.0) + else: + # for other bins, should have a weight + assert np.all(table_w["weight"][mask] >= 0.0) From 6fc63facbe5e382af018d105ad75c1325d4342e4 Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Wed, 28 Jan 2026 13:32:24 +0100 Subject: [PATCH 08/16] drop polar weighter, since not needed currently --- docs/changes/2927.feature.rst | 7 +- src/ctapipe/irf/__init__.py | 2 - src/ctapipe/irf/event_weighter.py | 69 +------------------- src/ctapipe/irf/tests/test_event_weighter.py | 48 -------------- 4 files changed, 5 insertions(+), 121 deletions(-) diff --git a/docs/changes/2927.feature.rst b/docs/changes/2927.feature.rst index 6e65960c4bd..a37be7c1566 100644 --- a/docs/changes/2927.feature.rst +++ b/docs/changes/2927.feature.rst @@ -1,10 +1,11 @@ -Introduces a new Component `~ctapipe.irf.EventWeighter` that has different implementations of spectral event weighting: +Introduces a new Component `~ctapipe.irf.EventWeighter` that has different +implementations of spectral event weighting: * `SimpleEventWeigher`: spectral weighting for the full FOV * `RadialEventWeighter`: spectral weighing in radial bins in the FOV -* `PolarEventWeighter`: spectral weighting in radial+azimuthal (r,phi) bins in the FOV -They operate on a pre-processed table of DL2 information, where you specify the energy and FOV coordinate columns to use. The interface is as follows: +They operate on a pre-processed table of DL2 information, where you specify the +energy and FOV coordinate columns to use. The interface is as follows: .. code-block:: python diff --git a/src/ctapipe/irf/__init__.py b/src/ctapipe/irf/__init__.py index 6c851611238..8acf1ce2747 100644 --- a/src/ctapipe/irf/__init__.py +++ b/src/ctapipe/irf/__init__.py @@ -20,7 +20,6 @@ ) from .event_weighter import ( EventWeighter, - PolarEventWeighter, RadialEventWeighter, SimpleEventWeighter, ) @@ -67,7 +66,6 @@ "check_bins_in_range", "make_bins_per_decade", "EventWeighter", - "PolarEventWeighter", "RadialEventWeighter", "SimpleEventWeighter", "spectrum_from_simulation_config", diff --git a/src/ctapipe/irf/event_weighter.py b/src/ctapipe/irf/event_weighter.py index d40b58fe4ff..98c756a4ccc 100644 --- a/src/ctapipe/irf/event_weighter.py +++ b/src/ctapipe/irf/event_weighter.py @@ -4,7 +4,6 @@ from abc import abstractmethod from collections.abc import Callable -from itertools import product from typing import override import numpy as np @@ -16,14 +15,13 @@ from ..core import Component, traits from ..core.feature_generator import _shallow_copy_table -from .binning import DefaultFoVOffsetBins, DefaultFoVPhiBins +from .binning import DefaultFoVOffsetBins from .spectra import Spectra, spectrum_from_name __all__ = [ "EventWeighter", "SimpleEventWeighter", "RadialEventWeighter", - "PolarEventWeighter", ] @@ -154,68 +152,3 @@ def _compute_weights(self, events_table: QTable): ) events_table.meta["OFFSBINS"] = list(offset_bins.to_value("deg")) - - -class PolarEventWeighter(EventWeighter, DefaultFoVOffsetBins, DefaultFoVPhiBins): - """ - Weights in field-of-view polar $(r, \\phi)$ bins in the `~ctapipe.coordinates.NominalFrame`. - - Calling this class adds a column to the output table with the event-wise - spectral-spatial weights, with column name ``weight_column``. This - implementation additionally adds the columns - ``output_table["fov_offset_bin"]``, ``output_table["fov_phi_bin"]``, and the - list of offset and phi bin edges in ``output_table.meta["OFFSBINS"]`` and - ``output_table.meta["PHIBINS"]`` respectively. - """ - - fov_offset_column = traits.Unicode( - help="name of FOV radial offset column", default_value="true_fov_offset" - ).tag(config=True) - - fov_phi_column = traits.Unicode( - help="name of the FOV azimuthal coordinate", default_value="true_fov_phi" - ).tag(config=True) - - @override - def _compute_weights(self, events_table: QTable): - offset_bins = self.fov_offset_bins - phi_bins = self.fov_phi_bins - offset = events_table[self.fov_offset_column].to_value(offset_bins.unit) - phi = events_table[self.fov_phi_column].to_value(phi_bins.unit) - energy = events_table[self.energy_column] - weights = np.zeros_like(energy.value) - - r_bin = np.digitize(offset, offset_bins.value) - phi_bin = np.digitize(phi, phi_bins.value) - - for i_r, i_phi in product(range(1, len(offset_bins)), range(1, len(phi_bins))): - mask = (r_bin == i_r) & (phi_bin == i_phi) - phi_solid_angle_fraction = (phi_bins[i_phi] - phi_bins[i_phi - 1]) / ( - 360.0 * u.deg - ) - print(i_phi, phi_solid_angle_fraction) - weights[mask] = calculate_event_weights( - true_energy=energy[mask], - target_spectrum=self.target_spectrum, - simulated_spectrum=lambda E: self.source_spectrum.integrate_cone( - offset_bins[i_r - 1], offset_bins[i_r] - )(E) - * phi_solid_angle_fraction, - ) - - events_table[self.weight_column] = weights - events_table["fov_offset_bin"] = r_bin - events_table["fov_phi_bin"] = phi_bin - - events_table.columns["fov_offset_bin"].description = ( - "Bin i defined as offset[i-1] <= fov_offset < offset[i]. " - "Where offset is `OFFSBINS` array found this table's metadata." - ) - - events_table.columns["fov_phi_bin"].description = ( - "Bin i defined as phi[i-1] <= fov_phi < phi[i]. " - "Where phi is `PHIBINS` array found this table's metadata." - ) - - events_table.meta["OFFSBINS"] = list(offset_bins.to_value("deg")) - events_table.meta["PHIBINS"] = list(phi_bins.to_value("deg")) diff --git a/src/ctapipe/irf/tests/test_event_weighter.py b/src/ctapipe/irf/tests/test_event_weighter.py index 305cb1c8232..fe817aca17a 100644 --- a/src/ctapipe/irf/tests/test_event_weighter.py +++ b/src/ctapipe/irf/tests/test_event_weighter.py @@ -6,7 +6,6 @@ from astropy.table import QTable from ctapipe.irf.event_weighter import ( - PolarEventWeighter, RadialEventWeighter, SimpleEventWeighter, ) @@ -138,50 +137,3 @@ def test_radial_weights_fits(example_weight_table, tmp_path): 4.0, 5.0, ] - - -def test_polar_weights(example_weight_table): - from ctapipe.irf.spectra import PowerLaw - - table = example_weight_table - - source_spectrum = PowerLaw( - normalization=u.Quantity(0.00027, "TeV-1 s-1 sr-1 m-2"), - index=-2.0, - e_ref=1.0 * u.TeV, - ) - - weight = PolarEventWeighter( - source_spectrum=source_spectrum, - target_spectrum_name="CRAB_HEGRA", - fov_offset_min=0.0 * u.deg, - fov_offset_max=5.0 * u.deg, - fov_offset_n_bins=5, - fov_phi_n_bins=4, - ) - - assert np.allclose(weight.fov_offset_bins, [0.0, 1.0, 2.0, 3.0, 4.0, 5.0] * u.deg) - assert np.allclose(weight.fov_phi_bins, [0.0, 90.0, 180.0, 270.0, 360.0] * u.deg) - - table_w = weight(table) - - assert "fov_offset_bin" in table_w.colnames - assert "fov_phi_bin" in table_w.colnames - assert "OFFSBINS" in table_w.meta - assert "PHIBINS" in table_w.meta - assert table_w.meta["OFFSBINS"] == [0.0, 1.0, 2.0, 3.0, 4.0, 5.0] - assert table_w.meta["PHIBINS"] == [0.0, 90.0, 180.0, 270.0, 360.0] - - for ii in range(0, 7): - for jj in range(0, 4): - mask = (table_w["fov_offset_bin"] == ii) & (table_w["fov_phi_bin"] == jj) - - if ii == 0 or jj == 0: - # there should be no 0th bin - assert len(table_w[mask]) == 0 - elif ii == 6 or jj == 5: - # for outlier bin 6,5, weights should be all 0 - assert np.all(table_w["weight"][mask] == 0.0) - else: - # for other bins, should have a weight - assert np.all(table_w["weight"][mask] >= 0.0) From 74116b7edad647d66dfd9e84fb227595ab95802f Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Wed, 28 Jan 2026 14:31:21 +0100 Subject: [PATCH 09/16] use PyIRF binning conventions and function --- src/ctapipe/irf/event_weighter.py | 25 +++++++++++------ src/ctapipe/irf/tests/test_event_weighter.py | 29 ++++++++++---------- 2 files changed, 32 insertions(+), 22 deletions(-) diff --git a/src/ctapipe/irf/event_weighter.py b/src/ctapipe/irf/event_weighter.py index 98c756a4ccc..d7e5e325546 100644 --- a/src/ctapipe/irf/event_weighter.py +++ b/src/ctapipe/irf/event_weighter.py @@ -9,6 +9,7 @@ import numpy as np from astropy import units as u from astropy.table import QTable +from pyirf.binning import calculate_bin_indices from pyirf.spectral import ( calculate_event_weights, ) @@ -110,9 +111,12 @@ class RadialEventWeighter(EventWeighter, DefaultFoVOffsetBins): Weights in radial (FOV) offset bins in the `~ctapipe.coordinates.NominalFrame`. Calling this class adds a column to the output table with the event-wise - spectral-spatial weights, with column name ``weight_column``. This implementation - additionally adds the column ``output_table["fov_offset_bin"]``, and the - list of offset bin edges in ``output_table.meta["OFFSBINS"]`` + spectral-spatial weights, with column name ``weight_column``. This + implementation additionally adds the column + ``output_table["fov_offset_bin"]`` and + ``output_table["fov_offset_is_valid"]`` following the conventions of + ``pyirf.binning.calculate_bin_indices``, and the list of offset bin edges in + ``output_table.meta["OFFSBINS"]`` """ fov_offset_column = traits.Unicode( @@ -122,33 +126,38 @@ class RadialEventWeighter(EventWeighter, DefaultFoVOffsetBins): @override def _compute_weights(self, events_table: QTable): offset_bins = self.fov_offset_bins - offset = events_table[self.fov_offset_column].to_value(offset_bins.unit) + offset = events_table[self.fov_offset_column].to(offset_bins.unit) energy = events_table[self.energy_column] weights = np.zeros_like(energy.value) # note that the bin i from digitize starts at 1 and means: # offset_bins[i-1] <= offset < offset_bins[ii]) - r_bin = np.digitize(offset, offset_bins.value) + r_bin, r_valid = calculate_bin_indices(offset, offset_bins) - for ii in range(1, len(offset_bins)): + for ii in range(len(offset_bins) - 1): self.log.debug( - f"bin {ii} offset=[{offset_bins[ii - 1]}, {offset_bins[ii]})" + f"bin {ii} offset=[{offset_bins[ii]}, {offset_bins[ii + 1]})" ) mask = r_bin == ii weights[mask] = calculate_event_weights( true_energy=energy[mask], target_spectrum=self.target_spectrum, simulated_spectrum=self.source_spectrum.integrate_cone( - offset_bins[ii - 1], offset_bins[ii] + offset_bins[ii], offset_bins[ii + 1] ), ) events_table[self.weight_column] = weights events_table["fov_offset_bin"] = r_bin + events_table["fov_offset_is_valid"] = r_valid events_table.columns["fov_offset_bin"].description = ( "Bin i defined as offset[i-1] <= fov_offset < offset[i]. " "Where offset is `OFFSBINS` array found this table's metadata." ) + events_table.columns[ + "fov_offset_is_valid" + ].description = "True if event's offset was inside the binning range." + events_table.meta["OFFSBINS"] = list(offset_bins.to_value("deg")) diff --git a/src/ctapipe/irf/tests/test_event_weighter.py b/src/ctapipe/irf/tests/test_event_weighter.py index fe817aca17a..a847b837061 100644 --- a/src/ctapipe/irf/tests/test_event_weighter.py +++ b/src/ctapipe/irf/tests/test_event_weighter.py @@ -15,9 +15,8 @@ def example_weight_table(): table = QTable( dict( - true_energy=[1.0, 2.0, 0.5, 0.2] * u.TeV, - true_fov_offset=[0.1, 1.2, 2.2, 3.2] * u.deg, - true_fov_phi=[45.0, 45.0, 268.0, 120.0] * u.deg, + true_energy=[1.0, 2.0, 0.5, 0.2, 4.0, 3.2] * u.TeV, + true_fov_offset=[0.1, 1.2, 2.2, 3.2, 10.0, 4.1] * u.deg, ) ) table.meta["VERSION"] = 1.0 @@ -64,6 +63,8 @@ def test_flat_weighting(example_weight_table): def test_radial_weights(example_weight_table): + from pyirf.binning import OVERFLOW_INDEX + from ctapipe.irf.spectra import PowerLaw table = example_weight_table @@ -87,21 +88,21 @@ def test_radial_weights(example_weight_table): table_w = weight(table) assert "fov_offset_bin" in table_w.colnames + assert "fov_offset_is_valid" in table_w.colnames assert "OFFSBINS" in table_w.meta assert table_w.meta["OFFSBINS"] == [0.0, 1.0, 2.0, 3.0, 4.0, 5.0] - for ii in range(0, 7): - mask = table_w["fov_offset_bin"] == ii + assert np.all( + np.unique(table_w["fov_offset_bin"].data) + == np.array([0, 1, 2, 3, 4, OVERFLOW_INDEX]) + ) - if ii == 0: - # there should be no 0th bin - assert len(table_w[mask]) == 0 - elif ii == 6: - # for outlier bin 6, weights should be all 0 - assert np.all(table_w["weight"][mask] == 0.0) - else: - # for bins 1-5,weights should be positive in this case - assert np.all(table_w["weight"][mask] >= 0.0) + for ii in range(0, 5): + mask = table_w["fov_offset_bin"] == ii + assert np.all(table_w["weight"][mask] >= 0.0) + assert np.all( + weight.fov_offset_bins[ii] < table_w["true_fov_offset"][mask] + ) & np.all(table_w["true_fov_offset"][mask] <= weight.fov_offset_bins[ii + 1]) def test_radial_weights_fits(example_weight_table, tmp_path): From a871d9a43d0d6235d9e31e2bfa5baba3d39d3648 Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Wed, 28 Jan 2026 14:39:24 +0100 Subject: [PATCH 10/16] add over/underflow metadata --- src/ctapipe/irf/event_weighter.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/ctapipe/irf/event_weighter.py b/src/ctapipe/irf/event_weighter.py index d7e5e325546..3fb4bfc0612 100644 --- a/src/ctapipe/irf/event_weighter.py +++ b/src/ctapipe/irf/event_weighter.py @@ -9,7 +9,7 @@ import numpy as np from astropy import units as u from astropy.table import QTable -from pyirf.binning import calculate_bin_indices +from pyirf.binning import OVERFLOW_INDEX, UNDERFLOW_INDEX, calculate_bin_indices from pyirf.spectral import ( calculate_event_weights, ) @@ -161,3 +161,5 @@ def _compute_weights(self, events_table: QTable): ].description = "True if event's offset was inside the binning range." events_table.meta["OFFSBINS"] = list(offset_bins.to_value("deg")) + events_table.meta["BINOVER"] = OVERFLOW_INDEX + events_table.meta["BINUNDR"] = UNDERFLOW_INDEX From 7d4e19e1009f18131d11d6a9b7294c1a5247af11 Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Wed, 28 Jan 2026 16:58:28 +0100 Subject: [PATCH 11/16] fix changelog links --- docs/changes/2927.feature.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/changes/2927.feature.rst b/docs/changes/2927.feature.rst index a37be7c1566..9454f54e507 100644 --- a/docs/changes/2927.feature.rst +++ b/docs/changes/2927.feature.rst @@ -1,8 +1,8 @@ Introduces a new Component `~ctapipe.irf.EventWeighter` that has different implementations of spectral event weighting: -* `SimpleEventWeigher`: spectral weighting for the full FOV -* `RadialEventWeighter`: spectral weighing in radial bins in the FOV +* `~ctapipe.irf.SimpleEventWeighter`: spectral weighting for the full FOV +* `~ctapipe.irf.RadialEventWeighter`: spectral weighing in radial bins in the FOV They operate on a pre-processed table of DL2 information, where you specify the energy and FOV coordinate columns to use. The interface is as follows: From 8b67c7805da489c8aaef40a268ce541e8edf8ed0 Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Wed, 28 Jan 2026 17:26:26 +0100 Subject: [PATCH 12/16] use new shallow_copy_table --- src/ctapipe/irf/event_weighter.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/ctapipe/irf/event_weighter.py b/src/ctapipe/irf/event_weighter.py index 3fb4bfc0612..fdd41af7b7d 100644 --- a/src/ctapipe/irf/event_weighter.py +++ b/src/ctapipe/irf/event_weighter.py @@ -8,14 +8,14 @@ import numpy as np from astropy import units as u -from astropy.table import QTable +from astropy.table import QTable, Table from pyirf.binning import OVERFLOW_INDEX, UNDERFLOW_INDEX, calculate_bin_indices from pyirf.spectral import ( calculate_event_weights, ) from ..core import Component, traits -from ..core.feature_generator import _shallow_copy_table +from ..core.feature_generator import shallow_copy_table from .binning import DefaultFoVOffsetBins from .spectra import Spectra, spectrum_from_name @@ -69,10 +69,10 @@ def _compute_weights(self, events_table: QTable): f"{self.__class__.__name__} weighting is not implemented" ) - def __call__(self, events_table: QTable) -> QTable: + def __call__(self, events_table: Table | QTable) -> QTable: """Returns shallow copy of input table with a ``weight`` column added""" - table = _shallow_copy_table(events_table) + table = shallow_copy_table(events_table, output_cls=QTable) self._compute_weights(table) return table From e0e693104e9507a2a6a2cbeafb30e8d52a5a26d2 Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Fri, 30 Jan 2026 10:20:08 +0100 Subject: [PATCH 13/16] added example to module docs --- docs/api-reference/irf/event_weighter.rst | 12 +++++ docs/api-reference/irf/index.rst | 1 + docs/api-reference/irf/irfs.rst | 2 +- src/ctapipe/irf/event_weighter.py | 54 ++++++++++++++++++++++- 4 files changed, 67 insertions(+), 2 deletions(-) create mode 100644 docs/api-reference/irf/event_weighter.rst diff --git a/docs/api-reference/irf/event_weighter.rst b/docs/api-reference/irf/event_weighter.rst new file mode 100644 index 00000000000..3c2df2c0295 --- /dev/null +++ b/docs/api-reference/irf/event_weighter.rst @@ -0,0 +1,12 @@ +.. _event_weighter: + +*************** +Event Weighting +*************** + + +Reference/ API +============== + +.. automodapi:: ctapipe.irf.event_weighter + :no-inheritance-diagram: diff --git a/docs/api-reference/irf/index.rst b/docs/api-reference/irf/index.rst index fab2810178f..95d2c4d6b56 100644 --- a/docs/api-reference/irf/index.rst +++ b/docs/api-reference/irf/index.rst @@ -34,6 +34,7 @@ Submodules benchmarks binning spectra + event_weighter Reference/API diff --git a/docs/api-reference/irf/irfs.rst b/docs/api-reference/irf/irfs.rst index 9755f91ec70..cad2cc44fd2 100644 --- a/docs/api-reference/irf/irfs.rst +++ b/docs/api-reference/irf/irfs.rst @@ -1,7 +1,7 @@ .. _irfs: ************** -IRF components +IRF generation ************** diff --git a/src/ctapipe/irf/event_weighter.py b/src/ctapipe/irf/event_weighter.py index fdd41af7b7d..000e6fc121b 100644 --- a/src/ctapipe/irf/event_weighter.py +++ b/src/ctapipe/irf/event_weighter.py @@ -1,6 +1,58 @@ #!/usr/bin/env python3 -"""Components for computing event weights.""" +""" +Components for computing event weights. + +`~ctapipe.irf.EventWeighter` is a base Component with implementations of spectral event weighting: + +* `~ctapipe.irf.SimpleEventWeighter`: spectral weighting for the full FOV +* `~ctapipe.irf.RadialEventWeighter`: spectral weighing in radial bins in the FOV + +They operate on a pre-processed table of DL2 information, where you specify the +energy and FOV coordinate columns to use. The interface is as follows: + +.. code-block:: python + + from astropy.table import QTable + + from ctapipe.irf import RadialEventWeighter spectrum_from_name + + table = QTable( + dict( + true_energy=[1.0, 2.0, 0.5, 0.2] * u.TeV, + true_fov_offset=[0.1, 1.2, 2.2, 3.2] * u.deg, + ) + ) + + # the source spectrum can be anything, but is usually loaded + # with ctapipe.irf.spectrum_from_simulation_config() to + # weight simulations. Here we will just use the electron + # spectrum as a demo: + source_spectrum = spectrum_from_name("IRFDOC_ELECTRON_SPECTRUM") + + weighter = RadialEventWeighter( + source_spectrum=source_spectrum + target_spectrum_name="CRAB_HEGRA", + fov_offset_max=5.0*u.deg, + fov_offset_n_bins=5, + ) + + table_with_weights = weighter(table) + print(table_with_weights) + + +:: + + true_energy true_fov_offset weight fov_offset_bin + TeV deg + ----------- --------------- ------------------ -------------- + 1.0 0.1 12.399508446433442 1 + 2.0 1.2 7.247055871572714 2 + 0.5 2.2 1.4149219056909543 3 + 0.2 3.2 0.4812883464910382 4 + + +""" from abc import abstractmethod from collections.abc import Callable From 3f55e5211a9dc52f11b9d03b289f453202fa74e3 Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Tue, 3 Feb 2026 14:11:24 +0100 Subject: [PATCH 14/16] add better docstring --- src/ctapipe/irf/event_weighter.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/ctapipe/irf/event_weighter.py b/src/ctapipe/irf/event_weighter.py index 000e6fc121b..dffb4d35233 100644 --- a/src/ctapipe/irf/event_weighter.py +++ b/src/ctapipe/irf/event_weighter.py @@ -132,6 +132,9 @@ def __call__(self, events_table: Table | QTable) -> QTable: class SimpleEventWeighter(EventWeighter): """Weights all events spectrally with no spatial binning. + This should be used for point-like signal gammas and diffuse background + events (protons, electrons). + Calling this class adds a column to the output table with the event-wise spectral weights, with column name ``weight_column``. """ @@ -162,6 +165,9 @@ class RadialEventWeighter(EventWeighter, DefaultFoVOffsetBins): """ Weights in radial (FOV) offset bins in the `~ctapipe.coordinates.NominalFrame`. + This should only be used for diffuse signal (gamma) events where you want to + weight them to a point-source spectrum in radial offset bins. + Calling this class adds a column to the output table with the event-wise spectral-spatial weights, with column name ``weight_column``. This implementation additionally adds the column From 408bb299aae3edd70abd8dfbbffb1f67826e9d76 Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Tue, 3 Feb 2026 14:18:05 +0100 Subject: [PATCH 15/16] simplify uniqueness check --- src/ctapipe/irf/spectra.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/ctapipe/irf/spectra.py b/src/ctapipe/irf/spectra.py index 496caa129f0..ae670172cc3 100644 --- a/src/ctapipe/irf/spectra.py +++ b/src/ctapipe/irf/spectra.py @@ -101,18 +101,19 @@ def spectrum_from_simulation_config( # Currently we do not support those, so we raise exceptions here to # avoid that we incorrectly compute the effective area, which would have # a high scientific impact. - for itm in [ + cols = [ "spectral_index", "energy_range_min", "energy_range_max", "max_scatter_range", "max_viewcone_radius", "min_viewcone_radius", - ]: - if len(np.unique(simulation_configuration_table[itm])) > 1: - raise NotImplementedError( - f"Unsupported: '{itm}' differs across simulation runs" - ) + ] + + if len(np.unique(simulation_configuration_table[cols], axis=0)) > 1: + raise NotImplementedError( + f"Unsupported: {', '.join(cols)} must not differ across simulation runs" + ) n_showers_config = ( simulation_configuration_table["n_showers"] From f5ce1b15e3ab77695df970b6f857e230e6fccedc Mon Sep 17 00:00:00 2001 From: Karl Kosack Date: Wed, 4 Mar 2026 09:48:42 +0100 Subject: [PATCH 16/16] bump