Skip to content

Commit 6af3bc3

Browse files
kosackkosack
committed
Refactor event weighting into irf.EventWeighter
Weighting should not be in ``ctapipe.io``, as it's really IRF specific. So I moved it to ctapipe.irf and created a EventWeighter class hierarchy for the different methods.
1 parent 949cd82 commit 6af3bc3

5 files changed

Lines changed: 442 additions & 11 deletions

File tree

‎src/ctapipe/irf/__init__.py‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,13 @@
3131
PointSourceSensitivityOptimizer,
3232
ThetaPercentileCutCalculator,
3333
)
34-
from .spectra import ENERGY_FLUX_UNIT, FLUX_UNIT, SPECTRA, Spectra
34+
from .spectra import (
35+
ENERGY_FLUX_UNIT,
36+
FLUX_UNIT,
37+
SPECTRA,
38+
Spectra,
39+
spectrum_from_simulation_config,
40+
)
3541

3642
__all__ = [
3743
"AngularResolution2dMaker",
@@ -53,4 +59,5 @@
5359
"FLUX_UNIT",
5460
"check_bins_in_range",
5561
"make_bins_per_decade",
62+
"spectrum_from_simulation_config",
5663
]

‎src/ctapipe/irf/binning.py‎

Lines changed: 7 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -179,9 +179,13 @@ class DefaultFoVOffsetBins(Component):
179179
default_value=1,
180180
).tag(config=True)
181181

182-
def __init__(self, config=None, parent=None, **kwargs):
183-
super().__init__(config=config, parent=parent, **kwargs)
184-
self.fov_offset_bins = u.Quantity(
182+
# def __init__(self, config=None, parent=None, **kwargs):
183+
# super().__init__(config=config, parent=parent, **kwargs)
184+
#
185+
#
186+
@property
187+
def fov_offset_bins(self):
188+
return u.Quantity(
185189
np.linspace(
186190
self.fov_offset_min.to_value(u.deg),
187191
self.fov_offset_max.to_value(u.deg),

‎src/ctapipe/irf/event_weighter.py‎

Lines changed: 185 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,185 @@
1+
#!/usr/bin/env python3
2+
3+
"""Components for computing event weights."""
4+
5+
from abc import abstractmethod
6+
from collections.abc import Callable
7+
from typing import override
8+
9+
import numpy as np
10+
from astropy import units as u
11+
from astropy.table import QTable
12+
from pyirf.spectral import (
13+
calculate_event_weights,
14+
)
15+
16+
from ..core import Component, traits
17+
from ..core.feature_generator import _shallow_copy_table
18+
from .binning import DefaultFoVOffsetBins
19+
from .spectra import SPECTRA, Spectra
20+
21+
__all__ = [
22+
"EventWeighter",
23+
"SimpleEventWeighter",
24+
"RadialEventWeighter",
25+
"PolarEventWeighter",
26+
]
27+
28+
29+
class EventWeighter(Component):
30+
"""Compute weights to go from source to target spectra."""
31+
32+
target_spectrum_name = traits.UseEnum(
33+
Spectra,
34+
default_value=Spectra.CRAB_HEGRA,
35+
help="Pre-defined source spectrum to reweight to.",
36+
).tag(config=True)
37+
38+
is_diffuse = traits.Bool(
39+
default_value=True, help="If True, assume the source is diffuse."
40+
).tag(config=True)
41+
42+
energy_column = traits.Unicode(
43+
help="name of energy column", default_value="true_energy"
44+
).tag(config=True)
45+
weight_column = traits.Unicode(
46+
help="name of output weight column", default_value="weight"
47+
).tag(config=True)
48+
49+
def __init__(
50+
self,
51+
source_spectrum: Callable[[u.Quantity], u.Quantity],
52+
config=None,
53+
parent=None,
54+
**kwargs,
55+
):
56+
"""
57+
Parameters
58+
----------
59+
source_spectrum: Callable
60+
initial spectrum of the events to be processed.
61+
target_spectrum: Callable | None
62+
target spectrum to weight to. If None, a pre-defined spectrum
63+
function from the `target_spectrum_name` attribute will be used`
64+
"""
65+
super().__init__(config=config, parent=parent, **kwargs)
66+
self.source_spectrum = source_spectrum
67+
self.target_spectrum = SPECTRA[self.target_spectrum_name]
68+
69+
@abstractmethod
70+
def _compute_weights(self, events_table: QTable):
71+
raise NotImplementedError(
72+
f"{self.__class__.__name__} weighting is not implemented"
73+
)
74+
75+
def __call__(self, events_table: QTable) -> QTable:
76+
"""Returns shallow copy of input table with a `weight` column added"""
77+
78+
table = _shallow_copy_table(events_table)
79+
self._compute_weights(table)
80+
return table
81+
82+
83+
class SimpleEventWeighter(EventWeighter):
84+
"""Weights all events spectrally with no spatial binning.
85+
86+
Calling this class adds a column to the output table with the event-wise
87+
spectral weights, with column name ``weight_column``.
88+
"""
89+
90+
fov_offset_max = traits.AstroQuantity(
91+
help="upper bound of spatial integral applied to source function",
92+
default_value=u.Quantity(10, u.deg),
93+
physical_type=u.physical.angle,
94+
).tag(config=True)
95+
96+
@override
97+
def _compute_weights(self, events_table: QTable):
98+
energy = events_table[self.energy_column]
99+
source_spectrum = self.source_spectrum
100+
if self.is_diffuse:
101+
source_spectrum = source_spectrum.integrate_cone(
102+
0 * u.deg, self.fov_offset_max
103+
)
104+
weights = calculate_event_weights(
105+
energy,
106+
target_spectrum=self.target_spectrum,
107+
simulated_spectrum=source_spectrum,
108+
)
109+
events_table[self.weight_column] = weights
110+
111+
112+
class RadialEventWeighter(EventWeighter, DefaultFoVOffsetBins):
113+
"""
114+
Weights in radial (FOV) offset bins in the `~ctapipe.coordinates.NominalFrame`.
115+
116+
Calling this class adds a column to the output table with the event-wise
117+
spectral-spatial weights, with column name `weight_column``. This implementation
118+
additionally adds the column ``output_table["fov_offset_bin"]``, and the
119+
list of offset bin edges in ``output_table.meta["OFFSBINS"]``
120+
"""
121+
122+
fov_offset_column = traits.Unicode(
123+
help="name of FOV radial offset column", default_value="true_fov_offset"
124+
).tag(config=True)
125+
126+
@override
127+
def _compute_weights(self, events_table: QTable):
128+
offset_bins = self.fov_offset_bins
129+
offset = events_table[self.fov_offset_column].to_value(offset_bins.unit)
130+
energy = events_table[self.energy_column]
131+
weights = np.zeros_like(energy.value)
132+
133+
# note that the bin i from digitize starts at 1 and means:
134+
# offset_bins[i-1] <= offset < offset_bins[ii])
135+
r_bin = np.digitize(offset, offset_bins.value)
136+
137+
for ii in range(0, len(offset_bins)):
138+
self.log.debug(
139+
f"bin {ii} offset=[{offset_bins[ii - 1]}, {offset_bins[ii]})"
140+
)
141+
mask = r_bin == ii
142+
weights[mask] = calculate_event_weights(
143+
true_energy=energy[mask],
144+
target_spectrum=self.target_spectrum,
145+
simulated_spectrum=self.source_spectrum.integrate_cone(
146+
offset_bins[ii - 1], offset_bins[ii]
147+
),
148+
)
149+
150+
events_table[self.weight_column] = weights
151+
events_table["fov_offset_bin"] = r_bin
152+
153+
events_table.columns["fov_offset_bin"].description = (
154+
"Bin i defined as offset[i-1] <= fov_offset < offset[i]. "
155+
"Where offset is `OFFSBINS` array found this table's metadata."
156+
)
157+
158+
events_table.meta["OFFSBINS"] = list(offset_bins.to_value("deg"))
159+
160+
161+
class PolarEventWeighter(EventWeighter):
162+
"""
163+
Weights in field-of-view polar $(r, \\phi)$ bins in the `~ctapipe.coordinates.NominalFrame`.
164+
165+
Calling this class adds a column to the output table with the event-wise
166+
spectral-spatial weights, with column name `weight_column``. This
167+
implementation additionally adds the columns
168+
``output_table["fov_offset_bin"]``, ``output_table["fov_phi_bin"]``, and the
169+
list of offset and phi bin edges in ``output_table.meta["OFFSBINS"]`` and
170+
``output_table.meta["PHIBINS"]`` respectively.
171+
"""
172+
173+
fov_offset_column = traits.Unicode(
174+
help="name of FOV radial offset column", default_value="true_fov_offset"
175+
).tag(config=True)
176+
177+
fov_phi_column = traits.Unicode(
178+
help="name of the FOV azimuthal coordinate", default_value="true_fov_phi"
179+
).tag(config=True)
180+
181+
@override
182+
def _compute_weights(self, events_table: QTable):
183+
raise NotImplementedError(
184+
f"{self.__class__.__name__} weighting is not implemented"
185+
)

‎src/ctapipe/irf/spectra.py‎

Lines changed: 107 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1,26 +1,126 @@
11
"""Definition of spectra to be used to calculate event weights for irf computation"""
22

3-
from enum import Enum
3+
import logging
4+
from collections.abc import Callable
5+
from enum import StrEnum
46

57
import astropy.units as u
6-
from pyirf.spectral import CRAB_HEGRA, IRFDOC_ELECTRON_SPECTRUM, IRFDOC_PROTON_SPECTRUM
8+
import numpy as np
9+
from astropy.table import Table
10+
from pyirf.simulations import SimulatedEventsInfo
11+
from pyirf.spectral import (
12+
CRAB_HEGRA,
13+
IRFDOC_ELECTRON_SPECTRUM,
14+
IRFDOC_PROTON_SPECTRUM,
15+
PowerLaw,
16+
)
717

8-
__all__ = ["ENERGY_FLUX_UNIT", "FLUX_UNIT", "SPECTRA", "Spectra"]
18+
__all__ = [
19+
"ENERGY_FLUX_UNIT",
20+
"FLUX_UNIT",
21+
"SPECTRA",
22+
"Spectra",
23+
"spectrum_from_simulation_config",
24+
]
925

1026
ENERGY_FLUX_UNIT = (1 * u.erg / u.s / u.cm**2).unit
1127
FLUX_UNIT = (1 / u.erg / u.s / u.cm**2).unit
1228

1329

14-
class Spectra(Enum):
30+
class Spectra(StrEnum):
1531
"""Spectra for calculating event weights"""
1632

17-
CRAB_HEGRA = 1
18-
IRFDOC_ELECTRON_SPECTRUM = 2
19-
IRFDOC_PROTON_SPECTRUM = 3
33+
CRAB_HEGRA = "CRAB_HEGRA"
34+
IRFDOC_ELECTRON_SPECTRUM = "IRFDOC_ELECTRON_SPECTRUM"
35+
IRFDOC_PROTON_SPECTRUM = "IRFDOC_PROTON_SPECTRUM"
2036

2137

2238
SPECTRA = {
2339
Spectra.CRAB_HEGRA: CRAB_HEGRA,
2440
Spectra.IRFDOC_ELECTRON_SPECTRUM: IRFDOC_ELECTRON_SPECTRUM,
2541
Spectra.IRFDOC_PROTON_SPECTRUM: IRFDOC_PROTON_SPECTRUM,
2642
}
43+
44+
45+
logger = logging.getLogger(__name__)
46+
47+
48+
def spectrum_from_simulation_config(
49+
simulation_configuration_table: Table,
50+
shower_distribution_table: Table | None,
51+
obs_time: u.Quantity,
52+
method: str = "powerlaw",
53+
) -> Callable[[u.Quantity], u.Quantity]:
54+
"""
55+
Return simulated spectrum function from configuration information..
56+
57+
Note this currently implements only PowerLaw spectra, but in the future will returnq a
58+
59+
Parameters
60+
----------
61+
simulation_configuration_table: Table
62+
table as read by TableLoader.read_simulation_configuration()
63+
shower_distribution_table: Table
64+
table of simulated shower distribution as read from TableLoader.read_shower_distribution()
65+
obs_time: u.Quantity['s']
66+
Observation time in a unit convertible to seconds.
67+
method: str
68+
"powerlaw" (for assuming powerlaw distributions), or "histogram" to return an
69+
interpolation function from the underlying distribution histogram.
70+
71+
Returns
72+
-------
73+
Callable:
74+
simulated spectrum function.
75+
"""
76+
77+
if method != "powerlaw":
78+
return NotImplementedError(f"Method '{method}' is not implemented")
79+
80+
n_showers: int = 0
81+
if shower_distribution_table:
82+
n_showers = shower_distribution_table["n_entries"].sum()
83+
else:
84+
logger.warning(
85+
"Simulation distributions were not found in the input files, "
86+
"falling back to estimating the number of showers from the "
87+
"simulation configuration."
88+
)
89+
90+
# Some tight consistency checks. Eventually we will support using the
91+
# arbitrary shower distribution and non-flat spatial distributions.
92+
# Currently we do not support those, so we raise exceptions here to
93+
# avoid that we incorrectly compute the effective area, which would have
94+
# a high scientific impact.
95+
for itm in [
96+
"spectral_index",
97+
"energy_range_min",
98+
"energy_range_max",
99+
"max_scatter_range",
100+
"max_viewcone_radius",
101+
"min_viewcone_radius",
102+
]:
103+
if len(np.unique(simulation_configuration_table[itm])) > 1:
104+
raise NotImplementedError(
105+
f"Unsupported: '{itm}' differs across simulation runs"
106+
)
107+
108+
n_showers_config = (
109+
simulation_configuration_table["n_showers"]
110+
* simulation_configuration_table["shower_reuse"]
111+
).sum()
112+
if n_showers == 0:
113+
n_showers = n_showers_config
114+
115+
assert n_showers_config == n_showers, "Inconsistent number of simulated showers"
116+
sim_info = SimulatedEventsInfo(
117+
n_showers=n_showers,
118+
energy_min=simulation_configuration_table["energy_range_min"].quantity[0],
119+
energy_max=simulation_configuration_table["energy_range_max"].quantity[0],
120+
max_impact=simulation_configuration_table["max_scatter_range"].quantity[0],
121+
spectral_index=simulation_configuration_table["spectral_index"][0],
122+
viewcone_max=simulation_configuration_table["max_viewcone_radius"].quantity[0],
123+
viewcone_min=simulation_configuration_table["min_viewcone_radius"].quantity[0],
124+
)
125+
126+
return PowerLaw.from_simulation(sim_info, obstime=obs_time)

0 commit comments

Comments
 (0)