Skip to content

Commit

Permalink
Rework event weight calculation
Browse files Browse the repository at this point in the history
  • Loading branch information
LukasBeiske committed Jul 24, 2024
1 parent f39847c commit 62dbcb8
Show file tree
Hide file tree
Showing 2 changed files with 34 additions and 21 deletions.
33 changes: 24 additions & 9 deletions src/ctapipe/irf/select.py
Original file line number Diff line number Diff line change
Expand Up @@ -122,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 @@ -258,26 +257,42 @@ def make_derived_columns(self, events: QTable, obs_conf: Table) -> QTable:
reco_nominal = reco.transform_to(nominal)
events["reco_fov_lon"] = -reco_nominal.fov_lon # minus for GADF
events["reco_fov_lat"] = reco_nominal.fov_lat
events["weight"] = 1.0
return events

def make_event_weights(
self,
events: QTable,
spectrum: PowerLaw,
fov_range: tuple[u.Quantity, u.Quantity],
fov_offset_bins: u.Quantity | None = None,
) -> QTable:
if (
self.kind == "gammas"
and self.target_spectrum.normalization.unit.is_equivalent(
spectrum.normalization.unit * u.sr
)
):
spectrum = spectrum.integrate_cone(fov_range[0], fov_range[1])
if fov_offset_bins is None:
raise ValueError(
"gamma_target_spectrum is point-like, but no fov offset bins "
"for the integration of the simulated diffuse spectrum was given."
)

events["weight"] = 1.0

for low, high in zip(fov_offset_bins[:-1], fov_offset_bins[1:]):
fov_mask = events["true_source_fov_offset"] >= low
fov_mask &= events["true_source_fov_offset"] < high

events[fov_mask]["weight"] = calculate_event_weights(
events[fov_mask]["true_energy"],
target_spectrum=self.target_spectrum,
simulated_spectrum=spectrum.integrate_cone(low, high),
)
else:
events["weight"] = calculate_event_weights(
events["true_energy"],
target_spectrum=self.target_spectrum,
simulated_spectrum=spectrum,
)

events["weight"] = calculate_event_weights(
events["true_energy"],
target_spectrum=self.target_spectrum,
simulated_spectrum=spectrum,
)
return events
22 changes: 10 additions & 12 deletions src/ctapipe/tools/make_irf.py
Original file line number Diff line number Diff line change
Expand Up @@ -484,19 +484,17 @@ def start(self):

self.log.debug("%s Precuts: %s" % (sel.kind, sel.epp.quality_criteria))
evs, cnt, meta = sel.load_preselected_events(self.chunk_size, self.obs_time)
# Only calculate real event weights if sensitivity or background should be computed
if self.do_benchmarks or self.do_background:
fov_conf = self.bkg if self.do_background else self.sens
for i in range(fov_conf.fov_offset_n_bins):
low = fov_conf.fov_offset_bins[i]
high = fov_conf.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),
# Only calculate event weights if background or sensitivity should be calculated.
if self.do_background:
# Sensitivity is only calculated, if do_background and do_benchmarks is true.
if self.do_benchmarks:
evs = sel.make_event_weights(
evs, meta["spectrum"], self.sens.fov_offset_bins
)
# If only background should be calculated,
# only calculate weights for protons and electrons.
elif sel.kind in ("protons", "electrons"):
evs = sel.make_event_weights(evs, meta["spectrum"])

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

0 comments on commit 62dbcb8

Please sign in to comment.