Skip to content

Commit

Permalink
Only compute event weights, if a sensitivity is computed
Browse files Browse the repository at this point in the history
  • Loading branch information
LukasBeiske committed Jul 18, 2024
1 parent 0e879d3 commit f5ad428
Show file tree
Hide file tree
Showing 4 changed files with 42 additions and 37 deletions.
26 changes: 12 additions & 14 deletions src/ctapipe/irf/select.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,6 @@
from ..core import Component, QualityQuery
from ..core.traits import List, Tuple, Unicode
from ..io import TableLoader
from .binning import ResultValidRange
from .spectra import SPECTRA, Spectra


Expand Down Expand Up @@ -123,7 +122,6 @@ def make_empty_table(self) -> QTable:
"theta",
"true_source_fov_offset",
"reco_source_fov_offset",
"weight",
]
units = {
"true_energy": u.TeV,
Expand Down Expand Up @@ -174,7 +172,7 @@ def __init__(self, kind: str, file: Path, target_spectrum: Spectra, **kwargs):
self.file = file

def load_preselected_events(
self, chunk_size: int, obs_time: u.Quantity, valid_fov
self, chunk_size: int, obs_time: u.Quantity
) -> tuple[QTable, int, dict]:
opts = dict(dl2=True, simulated=True)
with TableLoader(self.file, parent=self, **opts) as load:
Expand All @@ -186,9 +184,7 @@ def load_preselected_events(
for _, _, events in load.read_subarray_events_chunked(chunk_size, **opts):
selected = events[self.epp.get_table_mask(events)]
selected = self.epp.normalise_column_names(selected)
selected = self.make_derived_columns(
selected, spectrum, obs_conf, valid_fov
)
selected = self.make_derived_columns(selected, obs_conf)
bits.append(selected)
n_raw_events += len(events)

Expand Down Expand Up @@ -225,9 +221,7 @@ def get_metadata(
obs,
)

def make_derived_columns(
self, events: QTable, spectrum: PowerLaw, obs_conf: Table, valid_fov
) -> QTable:
def make_derived_columns(self, events: QTable, obs_conf: Table) -> QTable:
if obs_conf["subarray_pointing_lat"].std() < 1e-3:
assert all(obs_conf["subarray_pointing_frame"] == 0)
# Lets suppose 0 means ALTAZ
Expand Down Expand Up @@ -264,21 +258,25 @@ def make_derived_columns(
events["reco_fov_lon"] = -reco_nominal.fov_lon # minus for GADF
events["reco_fov_lat"] = reco_nominal.fov_lat

return events

def make_event_weights(
self,
events: QTable,
spectrum: PowerLaw,
fov_range: tuple[u.Quantity, u.Quantity],
) -> QTable:
if (
self.kind == "gammas"
and self.target_spectrum.normalization.unit.is_equivalent(
spectrum.normalization.unit * u.sr
)
):
if isinstance(valid_fov, ResultValidRange):
spectrum = spectrum.integrate_cone(valid_fov.min, valid_fov.max)
else:
spectrum = spectrum.integrate_cone(valid_fov[0], valid_fov[-1])
spectrum = spectrum.integrate_cone(fov_range[0], fov_range[1])

events["weight"] = calculate_event_weights(
events["true_energy"],
target_spectrum=self.target_spectrum,
simulated_spectrum=spectrum,
)

return events
6 changes: 4 additions & 2 deletions src/ctapipe/irf/tests/test_select.py
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,6 @@ def test_events_loader(gamma_diffuse_full_reco_file):
events, count, meta = loader.load_preselected_events(
chunk_size=10000,
obs_time=u.Quantity(50, u.h),
valid_fov=u.Quantity([0, 1], u.deg),
)

columns = [
Expand All @@ -163,9 +162,12 @@ def test_events_loader(gamma_diffuse_full_reco_file):
"theta",
"true_source_fov_offset",
"reco_source_fov_offset",
"weight",
]
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"], (0 * u.deg, 1 * u.deg))
assert "weight" in events.colnames
26 changes: 14 additions & 12 deletions src/ctapipe/tools/make_irf.py
Original file line number Diff line number Diff line change
Expand Up @@ -483,19 +483,21 @@ def start(self):
sel.epp.gammaness_classifier = self.opt_result.gh_cuts.meta["CLFNAME"]

self.log.debug("%s Precuts: %s" % (sel.kind, sel.epp.quality_criteria))
# TODO: This fov range is only used for the event weights for the sensitivity calculation.
# This should only be done if `do_benchmarks == True` and for each fov bin
# for which the sensitivity is calculated.
if self.do_benchmarks:
valid_fov = [self.sens.fov_offset_min, self.sens.fov_offset_max]
else:
valid_fov = [0, 5] * u.deg
evs, cnt, meta = sel.load_preselected_events(self.chunk_size, self.obs_time)
# Only calculate event weights if sensitivity should be computed
if self.do_benchmarks and self.do_background:
evs["weight"] = 1.0
for i in range(len(self.sens.fov_offset_bins) - 1):
low = self.sens.fov_offset_bins[i]
high = self.sens.fov_offset_bins[i + 1]
fov_mask = evs["true_source_fov_offset"] >= low
fov_mask &= evs["true_source_fov_offset"] < high
evs[fov_mask] = sel.make_event_weights(
evs[fov_mask],
meta["spectrum"],
(low, high),
)

evs, cnt, meta = sel.load_preselected_events(
chunk_size=self.chunk_size,
obs_time=self.obs_time,
valid_fov=valid_fov,
)
reduced_events[sel.kind] = evs
reduced_events[f"{sel.kind}_count"] = cnt
reduced_events[f"{sel.kind}_meta"] = meta
Expand Down
21 changes: 12 additions & 9 deletions src/ctapipe/tools/optimize_event_selection.py
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ class IrfEventSelector(Tool):
).tag(config=True)

obs_time = AstroQuantity(
default_value=50.0 * u.hour,
default_value=u.Quantity(50, u.hour),
physical_type=u.physical.time,
help="Observation time in the form ``<value> <unit>``",
).tag(config=True)
Expand Down Expand Up @@ -122,16 +122,19 @@ def setup(self):
]

def start(self):
# TODO: this event loading code seems to be largely repeated between all the tools,
# try to refactor to a common solution

reduced_events = dict()
for sel in self.particles:
evs, cnt, meta = sel.load_preselected_events(
self.chunk_size,
self.obs_time,
[self.optimizer.min_bkg_fov_offset, self.optimizer.max_bkg_fov_offset],
)
evs, cnt, meta = sel.load_preselected_events(self.chunk_size, self.obs_time)
if self.optimization_algorithm == "PointSourceSensitivityOptimizer":
evs = sel.make_event_weights(
evs,
meta["spectrum"],
(
self.optimizer.min_bkg_fov_offset,
self.optimizer.max_bkg_fov_offset,
),
)

reduced_events[sel.kind] = evs
reduced_events[f"{sel.kind}_count"] = cnt
if sel.kind == "gammas":
Expand Down

0 comments on commit f5ad428

Please sign in to comment.