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/docs/changes/2927.feature.rst b/docs/changes/2927.feature.rst new file mode 100644 index 00000000000..9454f54e507 --- /dev/null +++ b/docs/changes/2927.feature.rst @@ -0,0 +1,41 @@ +Introduces a new Component `~ctapipe.irf.EventWeighter` that has different +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, + ) + ) + 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 diff --git a/src/ctapipe/irf/__init__.py b/src/ctapipe/irf/__init__.py index d4abdb43fac..8acf1ce2747 100644 --- a/src/ctapipe/irf/__init__.py +++ b/src/ctapipe/irf/__init__.py @@ -18,6 +18,11 @@ check_bins_in_range, make_bins_per_decade, ) +from .event_weighter import ( + EventWeighter, + RadialEventWeighter, + SimpleEventWeighter, +) from .irfs import ( BackgroundRate2dMaker, EffectiveArea2dMaker, @@ -31,7 +36,14 @@ PointSourceSensitivityOptimizer, ThetaPercentileCutCalculator, ) -from .spectra import ENERGY_FLUX_UNIT, FLUX_UNIT, SPECTRA, Spectra +from .spectra import ( + ENERGY_FLUX_UNIT, + FLUX_UNIT, + SPECTRA, + Spectra, + spectrum_from_name, + spectrum_from_simulation_config, +) __all__ = [ "AngularResolution2dMaker", @@ -53,4 +65,9 @@ "FLUX_UNIT", "check_bins_in_range", "make_bins_per_decade", + "EventWeighter", + "RadialEventWeighter", + "SimpleEventWeighter", + "spectrum_from_simulation_config", + "spectrum_from_name", ] 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, + ) diff --git a/src/ctapipe/irf/event_weighter.py b/src/ctapipe/irf/event_weighter.py new file mode 100644 index 00000000000..dffb4d35233 --- /dev/null +++ b/src/ctapipe/irf/event_weighter.py @@ -0,0 +1,223 @@ +#!/usr/bin/env python3 + +""" +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 typing import override + +import numpy as np +from astropy import units as u +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 .binning import DefaultFoVOffsetBins +from .spectra import Spectra, spectrum_from_name + +__all__ = [ + "EventWeighter", + "SimpleEventWeighter", + "RadialEventWeighter", +] + + +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. + """ + super().__init__(config=config, parent=parent, **kwargs) + self.source_spectrum = source_spectrum + self.target_spectrum = spectrum_from_name(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: Table | QTable) -> QTable: + """Returns shallow copy of input table with a ``weight`` column added""" + + table = shallow_copy_table(events_table, output_cls=QTable) + self._compute_weights(table) + return table + + +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``. + """ + + 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`. + + 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 + ``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( + 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(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, r_valid = calculate_bin_indices(offset, offset_bins) + + for ii in range(len(offset_bins) - 1): + self.log.debug( + 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], 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")) + events_table.meta["BINOVER"] = OVERFLOW_INDEX + events_table.meta["BINUNDR"] = UNDERFLOW_INDEX diff --git a/src/ctapipe/irf/spectra.py b/src/ctapipe/irf/spectra.py index 75106112b97..ae670172cc3 100644 --- a/src/ctapipe/irf/spectra.py +++ b/src/ctapipe/irf/spectra.py @@ -1,22 +1,39 @@ """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", + "spectrum_from_name", +] 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 +41,96 @@ class Spectra(Enum): Spectra.IRFDOC_ELECTRON_SPECTRUM: IRFDOC_ELECTRON_SPECTRUM, Spectra.IRFDOC_PROTON_SPECTRUM: IRFDOC_PROTON_SPECTRUM, } + + +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__) + + +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. + 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[cols], axis=0)) > 1: + raise NotImplementedError( + f"Unsupported: {', '.join(cols)} must not differ 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) 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..a847b837061 --- /dev/null +++ b/src/ctapipe/irf/tests/test_event_weighter.py @@ -0,0 +1,140 @@ +#!/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, 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 + 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, spectrum_from_name + + table = example_weight_table + + weight = SimpleEventWeighter( + source_spectrum=spectrum_from_name(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 pyirf.binning import OVERFLOW_INDEX + + 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 "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] + + assert np.all( + np.unique(table_w["fov_offset_bin"].data) + == np.array([0, 1, 2, 3, 4, OVERFLOW_INDEX]) + ) + + 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): + """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, + ]