Skip to content

Commit

Permalink
Merge pull request #4599 from AtreyeeS/safemaker
Browse files Browse the repository at this point in the history 
Add a SafeMaskMaker at DL3 level
  • Loading branch information
@registerrier
registerrier committed Jan 24, 2024
2 parents 23352df + 8b75dc5 commit a711d11
Show file tree
Hide file tree
Showing 2 changed files with 132 additions and 68 deletions.
174 changes: 106 additions & 68 deletions gammapy/makers/safe.py
View file
Expand Up @@ -46,6 +46,9 @@ class SafeMaskMaker(Maker):
it uses the position of the center of the map by default.
offset_max : str or `~astropy.units.Quantity`
Maximum offset cut.
irfs : {"DL4", "DL3"}
Whether to use reprojected ("DL4") or raw ("DL3") irfs.
Default is "DL4".
"""

tag = "SafeMaskMaker"
Expand All @@ -65,6 +68,7 @@ def __init__(
position=None,
fixed_offset=None,
offset_max="3 deg",
irfs="DL4",
):
methods = set(methods)

Expand All @@ -83,6 +87,12 @@ def __init__(
"`position` and `fixed_offset` attributes are mutually exclusive"
)

if irfs not in ["DL3", "DL4"]:
ValueError(
"Invalid option for irfs: expected 'DL3' or 'DL4', got {irfs} instead."
)
self.irfs = irfs

def make_mask_offset_max(self, dataset, observation):
"""Make maximum offset mask.
Expand Down Expand Up @@ -123,7 +133,7 @@ def make_mask_energy_aeff_default(dataset, observation):
Safe data range mask.
"""
if observation is None:
raise ValueError("Method 'offset-max' requires an observation object.")
raise ValueError("Method 'aeff-default' requires an observation object.")

energy_max = observation.aeff.meta.get("HI_THRES", None)

Expand All @@ -145,6 +155,28 @@ def make_mask_energy_aeff_default(dataset, observation):

return dataset._geom.energy_mask(energy_min=energy_min, energy_max=energy_max)

def _get_offset(self, observation):
offset = self.fixed_offset
if offset is None:
if self.position:
offset = observation.get_pointing_icrs(observation.tmid).separation(
self.position
)
else:
offset = 0.0 * u.deg
return offset

def _get_position(self, observation, geom):
if self.fixed_offset is not None and observation is not None:
pointing = observation.get_pointing_icrs(observation.tmid)
return pointing.directional_offset_by(
position_angle=0 * u.deg, separation=self.fixed_offset
)
elif self.position is not None:
return self.position
else:
return geom.center_skydir

def make_mask_energy_aeff_max(self, dataset, observation=None):
"""Make safe energy mask from effective area maximum value.
Expand All @@ -160,47 +192,50 @@ def make_mask_energy_aeff_max(self, dataset, observation=None):
mask_safe : `~numpy.ndarray`
Safe data range mask.
"""
geom, exposure = dataset._geom, dataset.exposure

if self.fixed_offset is not None:
if observation:
position = observation.get_pointing_icrs(
observation.tmid
).directional_offset_by(
position_angle=0.0 * u.deg, separation=self.fixed_offset
)
else:
raise ValueError(
f"observation argument is mandatory with {self.fixed_offset}"
)
if self.fixed_offset is not None and observation is None:
raise ValueError(
f"{observation} argument is mandatory with {self.fixed_offset}"
)

elif self.position:
position = self.position
else:
position = geom.center_skydir
geom, exposure = dataset._geom, dataset.exposure

aeff = exposure.get_spectrum(position) / exposure.meta["livetime"]
if not np.any(aeff.data > 0.0):
log.warning(
f"Effective area is all zero at [{position.to_string('dms')}]. "
f"No safe energy band can be defined for the dataset '{dataset.name}': "
"setting `mask_safe` to all False."
if self.irfs == "DL3":
offset = self._get_offset(observation)

values = observation.aeff.evaluate(
offset=offset, energy_true=observation.aeff.axes["energy_true"].edges
)
return Map.from_geom(geom, data=False, dtype="bool")
valid = observation.aeff.axes["energy_true"].edges[
values > self.aeff_percent * np.max(values) / 100
]
energy_min = np.min(valid)

model = TemplateSpectralModel.from_region_map(aeff)
else:
position = self._get_position(observation, geom)

aeff = exposure.get_spectrum(position) / exposure.meta["livetime"]
if not np.any(aeff.data > 0.0):
log.warning(
f"Effective area is all zero at [{position.to_string('dms')}]. "
f"No safe energy band can be defined for the dataset '{dataset.name}': "
"setting `mask_safe` to all False."
)
return Map.from_geom(geom, data=False, dtype="bool")

energy_true = model.energy
energy_min = energy_true[np.where(model.values > 0)[0][0]]
energy_max = energy_true[-1]
model = TemplateSpectralModel.from_region_map(aeff)

aeff_thres = (self.aeff_percent / 100) * aeff.quantity.max()
inversion = model.inverse(
aeff_thres, energy_min=energy_min, energy_max=energy_max
)
energy_true = model.energy
energy_min = energy_true[np.where(model.values > 0)[0][0]]
energy_max = energy_true[-1]

aeff_thres = (self.aeff_percent / 100) * aeff.quantity.max()
inversion = model.inverse(
aeff_thres, energy_min=energy_min, energy_max=energy_max
)

if not np.isnan(inversion[0]):
energy_min = inversion[0]
if not np.isnan(inversion[0]):
energy_min = inversion[0]

return geom.energy_mask(energy_min=energy_min)

Expand All @@ -219,39 +254,27 @@ def make_mask_energy_edisp_bias(self, dataset, observation=None):
mask_safe : `~numpy.ndarray`
Safe data range mask.
"""
edisp, geom = dataset.edisp, dataset._geom
position = None

if self.fixed_offset is not None:
if observation:
pointing = observation.get_pointing_icrs(observation.tmid)
position = pointing.directional_offset_by(
position_angle=0 * u.deg, separation=self.fixed_offset
)
else:
raise ValueError(
f"{observation} argument is mandatory with {self.fixed_offset}"
)
if self.fixed_offset is not None and observation is None:
raise ValueError(
f"{observation} argument is mandatory with {self.fixed_offset}"
)

if isinstance(edisp, EDispKernelMap):
if position:
edisp = edisp.get_edisp_kernel(position=position)
else:
edisp = edisp.get_edisp_kernel(position=self.position)
else:
e_reco = dataset._geom.axes["energy"]
if position:
edisp = edisp.get_edisp_kernel(position=position, energy_axis=e_reco)
else:
edisp = edisp.get_edisp_kernel(
position=self.position, energy_axis=e_reco
)
edisp, geom = dataset.edisp, dataset._geom

energy_min = edisp.get_bias_energy(self.bias_percent / 100)
return geom.energy_mask(energy_min=energy_min[0])
if self.irfs == "DL3":
offset = self._get_offset(observation)
edisp = observation.edisp.to_edisp_kernel(offset)
else:
kwargs = dict()
kwargs["position"] = self._get_position(observation, geom)
if not isinstance(edisp, EDispKernelMap):
kwargs["energy_axis"] = dataset._geom.axes["energy"]
edisp = edisp.get_edisp_kernel(**kwargs)
energy_min = edisp.get_bias_energy(self.bias_percent / 100)[0]
return geom.energy_mask(energy_min=energy_min)

@staticmethod
def make_mask_energy_bkg_peak(dataset):
def make_mask_energy_bkg_peak(self, dataset, observation=None):
"""Make safe energy mask based on the binned background.
The energy threshold is defined as the lower edge of the energy
Expand All @@ -264,18 +287,28 @@ def make_mask_energy_bkg_peak(dataset):
----------
dataset : `~gammapy.datasets.MapDataset` or `~gammapy.datasets.SpectrumDataset`
Dataset to compute mask for.
observation: `~gammapy.data.Observation`
Observation to compute mask for. It is a mandatory argument when DL3 irfs are used.
Returns
-------
mask_safe : `~numpy.ndarray`
Safe data range mask.
"""
geom = dataset._geom
background_spectrum = dataset.npred_background().get_spectrum()
if self.irfs == "DL3":
bkg = observation.bkg.to_2d()
background_spectrum = np.ravel(
bkg.integral(axis_name="offset", offset=bkg.axes["offset"].bounds[1])
)
energy_axis = bkg.axes["energy"]
else:
background_spectrum = dataset.npred_background().get_spectrum()
energy_axis = geom.axes["energy"]

idx = np.argmax(background_spectrum.data, axis=0)
energy_axis = geom.axes["energy"]
energy_min = energy_axis.edges[idx]
return geom.energy_mask(energy_min=energy_min)
return geom.energy_mask(energy_min=energy_axis.edges[idx])

@staticmethod
def make_mask_bkg_invalid(dataset):
Expand Down Expand Up @@ -314,6 +347,11 @@ def run(self, dataset, observation=None):
dataset : `Dataset`
Dataset with defined safe range mask.
"""

if self.irfs == "DL3":
if observation is None:
raise ValueError("observation argument is mandatory with DL3 irfs")

if dataset.mask_safe:
mask_safe = dataset.mask_safe.data
else:
Expand All @@ -336,7 +374,7 @@ def run(self, dataset, observation=None):
mask_safe &= self.make_mask_energy_edisp_bias(dataset, observation)

if "bkg-peak" in self.methods:
mask_safe &= self.make_mask_energy_bkg_peak(dataset)
mask_safe &= self.make_mask_energy_bkg_peak(dataset, observation)

dataset.mask_safe = Map.from_geom(dataset._geom, data=mask_safe, dtype=bool)
return dataset
26 changes: 26 additions & 0 deletions gammapy/makers/tests/test_safe.py
View file
Expand Up @@ -157,6 +157,10 @@ def test_safe_mask_maker_aeff_max_fixed_observation(
)
assert_allclose(mask_aeff_max_bis.data.sum(), 0)

maker = SafeMaskMaker(irfs="DL3")
with pytest.raises(ValueError):
maker.run(dataset)


@requires_data()
def test_safe_mask_maker_aeff_max_fixed_offset(dataset, observation_cta_1dc):
Expand All @@ -172,6 +176,13 @@ def test_safe_mask_maker_aeff_max_fixed_offset(dataset, observation_cta_1dc):
with pytest.raises(ValueError):
mask_aeff_max = safe_mask_maker.make_mask_energy_aeff_max(dataset)

safe_mask_maker1 = SafeMaskMaker(
methods=["aeff-max"], aeff_percent=20, fixed_offset=None, irfs="DL3"
)

mask1 = safe_mask_maker1.make_mask_energy_aeff_max(dataset, observation_cta_1dc)
assert_allclose(mask1.data.sum(), 847)


@requires_data()
def test_safe_mask_maker_offset_max_fixed_offset(dataset, observation_cta_1dc):
Expand Down Expand Up @@ -215,6 +226,15 @@ def test_safe_mask_maker_edisp_bias_fixed_offset(dataset, observation_cta_1dc):
)
assert_allclose(mask_edisp_bias_offset.data.sum(), 1694)

safe_mask_maker1 = SafeMaskMaker(
irfs="DL3", bias_percent=0.02, fixed_offset=1.5 * u.deg
)

mask_edisp_bias_offset = safe_mask_maker1.make_mask_energy_edisp_bias(
dataset, observation_cta_1dc
)
assert_allclose(mask_edisp_bias_offset.data.sum(), 1694)


@requires_data()
def test_safe_mask_maker_bkg_peak(dataset, observation_cta_1dc):
Expand All @@ -224,6 +244,12 @@ def test_safe_mask_maker_bkg_peak(dataset, observation_cta_1dc):
mask_bkg_peak = safe_mask_maker.make_mask_energy_bkg_peak(dataset)
assert_allclose(mask_bkg_peak.data.sum(), 1936)

safe_mask_maker1 = SafeMaskMaker(fixed_offset=1.0 * u.deg, irfs="DL3")
mask_bkg_peak = safe_mask_maker1.make_mask_energy_bkg_peak(
dataset, observation_cta_1dc
)
assert_allclose(mask_bkg_peak.data.sum(), 1936)


@requires_data()
def test_safe_mask_maker_bkg_peak_first_bin(dataset, observation_cta_1dc):
Expand Down

0 comments on commit a711d11

Please sign in to comment.