From f5ad42807c773af6b284378152af2676ff9d0fc7 Mon Sep 17 00:00:00 2001 From: Lukas Beiske Date: Thu, 18 Jul 2024 17:44:50 +0200 Subject: [PATCH] Only compute event weights, if a sensitivity is computed --- src/ctapipe/irf/select.py | 26 +++++++++---------- src/ctapipe/irf/tests/test_select.py | 6 +++-- src/ctapipe/tools/make_irf.py | 26 ++++++++++--------- src/ctapipe/tools/optimize_event_selection.py | 21 ++++++++------- 4 files changed, 42 insertions(+), 37 deletions(-) diff --git a/src/ctapipe/irf/select.py b/src/ctapipe/irf/select.py index 8547ea57a8d..316e42677d1 100644 --- a/src/ctapipe/irf/select.py +++ b/src/ctapipe/irf/select.py @@ -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 @@ -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, @@ -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: @@ -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) @@ -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 @@ -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 diff --git a/src/ctapipe/irf/tests/test_select.py b/src/ctapipe/irf/tests/test_select.py index 73d16036aa4..a12fdd4a2d8 100644 --- a/src/ctapipe/irf/tests/test_select.py +++ b/src/ctapipe/irf/tests/test_select.py @@ -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 = [ @@ -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 diff --git a/src/ctapipe/tools/make_irf.py b/src/ctapipe/tools/make_irf.py index 870e6c676a3..0eb73572c93 100644 --- a/src/ctapipe/tools/make_irf.py +++ b/src/ctapipe/tools/make_irf.py @@ -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 diff --git a/src/ctapipe/tools/optimize_event_selection.py b/src/ctapipe/tools/optimize_event_selection.py index a5bc6184a91..7da9b320d0c 100644 --- a/src/ctapipe/tools/optimize_event_selection.py +++ b/src/ctapipe/tools/optimize_event_selection.py @@ -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 `` ``", ).tag(config=True) @@ -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":