From 7c04020cd35400594ccdeb230440e796e6cd140c Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Tue, 20 Jun 2023 19:44:51 +0200 Subject: [PATCH 01/11] add SafeMaker Signed-off-by: Atreyee Sinha --- gammapy/makers/__init__.py | 3 +- gammapy/makers/safe.py | 243 ++++++++++++++++++++++++++++++ gammapy/makers/tests/test_safe.py | 28 +++- 3 files changed, 272 insertions(+), 2 deletions(-) diff --git a/gammapy/makers/__init__.py b/gammapy/makers/__init__.py index d57f3d92d2..0a8f6d344d 100644 --- a/gammapy/makers/__init__.py +++ b/gammapy/makers/__init__.py @@ -13,7 +13,7 @@ from .core import Maker from .map import MapDatasetMaker from .reduce import DatasetsMaker -from .safe import SafeMaskMaker +from .safe import SafeMakerDL3, SafeMaskMaker from .spectrum import SpectrumDatasetMaker MAKER_REGISTRY = Registry( @@ -44,6 +44,7 @@ "RegionsFinder", "RingBackgroundMaker", "SafeMaskMaker", + "SafeMakerDL3", "SpectrumDatasetMaker", "WobbleRegionsFinder", ] diff --git a/gammapy/makers/safe.py b/gammapy/makers/safe.py index f3b42a33f2..17fd48d23d 100644 --- a/gammapy/makers/safe.py +++ b/gammapy/makers/safe.py @@ -340,3 +340,246 @@ def run(self, dataset, observation=None): dataset.mask_safe = Map.from_geom(dataset._geom, data=mask_safe, dtype=bool) return dataset + + +class SafeMakerDL3(Maker): + """Make safe data range mask for a given observation, directly from the DL3 + IRFs + + + .. warning:: + + Currently some methods computing a safe energy range ("aeff-default", + "aeff-max" and "edisp-bias") determine a true energy range and apply + it to reconstructed energy, effectively neglecting the energy dispersion. + + Parameters + ---------- + methods : {"aeff-default", "aeff-max", "offset-max", "bkg-peak"} + Method to use for the safe energy range. Can be a + list with a combination of those. Resulting masks + are combined with logical `and`. "aeff-default" + uses the energy ranged specified in the DL3 data + files, if available. + aeff_percent : float + Percentage of the maximal effective area to be used + as lower energy threshold for method "aeff-max". + bias_percent : float + Percentage of the energy bias to be used as lower + energy threshold for method "edisp-bias" + position : `~astropy.coordinates.SkyCoord` + Position at which the `aeff_percent` or `bias_percent` are computed. + fixed_offset : `~astropy.coordinates.Angle` + offset, calculated from the pointing position, at which + the `aeff_percent` or `bias_percent` are computed. + If neither the position nor fixed_offset is specified, + it uses the position of the center of the map by default. + offset_max : str or `~astropy.units.Quantity` + Maximum offset cut. + """ + + tag = "SafeMakerDL3" + available_methods = { + "aeff-default", + "aeff-max", + "edisp-bias", + "offset-max", + "bkg-peak", + } + + def __init__( + self, + methods=("aeff-default",), + aeff_percent=10, + bias_percent=10, + position=None, + fixed_offset=None, + offset_max="3 deg", + ): + methods = set(methods) + + if not methods.issubset(self.available_methods): + difference = methods.difference(self.available_methods) + raise ValueError(f"{difference} is not a valid method.") + + self.methods = methods + self.aeff_percent = aeff_percent / 100.0 + self.bias_percent = bias_percent / 100.0 + self.position = position + self.fixed_offset = fixed_offset + self.offset_max = Angle(offset_max) + if self.position and self.fixed_offset: + raise ValueError( + "`position` and `fixed_offset` attributes are mutually exclusive" + ) + + @staticmethod + def make_energy_aeff_default(observation): + """Select events based on from aeff default energies. + + Parameters + ---------- + observation: `~gammapy.data.Observation` + Observation to apply cuts on. + + Returns + ------- + events: `~gammapy.data.EventList` + Event list with applied cuts. + """ + + energy_max = observation.aeff.meta.get("HI_THRES", None) + + if energy_max: + energy_max = energy_max * u.TeV + else: + log.warning( + f"No default upper safe energy threshold defined for obs {observation.obs_id}" + ) + energy_max = observation.aeff.axes["energy_true"].edges[-1] + + energy_min = observation.aeff.meta.get("LO_THRES", None) + + if energy_min: + energy_min = energy_min * u.TeV + else: + log.warning( + f"No default lower safe energy threshold defined for obs {observation.obs_id}" + ) + energy_min = observation.aeff.axes["energy_true"].edges[0] + # TODO: Apply edisp + return observation.events.select_energy([energy_min, energy_max]) + + def make_energy_aeff_max(self, observation): + """Apply selection based on a minimal effective area + + TODO: Should this be offset dependent? + + Parameters + ---------- + observation: `~gammapy.data.Observation` + Observation to apply cuts on. + + Returns + ------- + events: `~gammapy.data.EventList` + Event list with applied cuts. + """ + + 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 + + values = observation.aeff.evaluate( + offset=offset, energy_true=observation.aeff.axes["energy_true"].edges + ) + valid = observation.aeff.axes["energy_true"].edges[ + values > self.aeff_percent * np.max(values) + ] + energy_min = np.min(valid) + energy_max = observation.aeff.axes["energy_true"].edges[-1] + # TODO: Apply edisp + return observation.events.select_energy([energy_min, energy_max]) + + def make_energy_bkg_peak(self, observation): + """Apply selection based the peak of the background IRF + + Parameters + ---------- + observation: `~gammapy.data.Observation` + Observation to apply cuts on. + + Returns + ------- + events: `~gammapy.data.EventList` + Event list with applied cuts. + """ + bkg = observation.bkg.to_2d() + + values = np.ravel( + bkg.integral(axis_name="offset", offset=bkg.axes["offset"].bounds[1]) + ) + energy_min = bkg.axes["energy"].center[values == np.max(values)] + energy_max = bkg.axes["energy"].center[-1] + return observation.events.select_energy([energy_min, energy_max]) + + def make_offset_max(self, observation): + """Apply selection based on maximum offset cut + + Parameters + ---------- + observation: `~gammapy.data.Observation` + Observation to apply cuts on. + + Returns + ------- + events: `~gammapy.data.EventList` + Event list with applied cuts. + """ + return observation.events.select_offset([0 * u.deg, self.offset_max]) + + def make_energy_edisp_bias(self, observation=None): + """Apply selection based on energy dispersion bias + + Parameters + ---------- + observation: `~gammapy.data.Observation` + Observation to apply cuts on. + + Returns + ------- + events: `~gammapy.data.EventList` + Event list with applied cuts. + """ + 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 + + edisp = observation.edisp.to_edisp_kernel(offset=offset) + energy_min = edisp.get_bias_energy(self.bias_percent)[0] + energy_max = observation.edisp.axes["energy_true"].center[-1] + return observation.events.select_energy([energy_min, energy_max]) + + def run(self, observation=None): + """Make safe data range mask. + + Parameters + ---------- + dataset : `~gammapy.datasets.MapDataset` or `~gammapy.datasets.SpectrumDataset` + Dataset to compute mask for. + observation: `~gammapy.data.Observation` + Observation to compute mask for. + + Returns + ------- + observation : `~gammapy.data.Observation + Observation with selected events + """ + + if "offset-max" in self.methods: + event_list = self.make_offset_max(observation) + observation = observation.copy(in_memory=True, events=event_list) + if "aeff-default" in self.methods: + event_list = self.make_energy_aeff_default(observation) + observation = observation.copy(in_memory=True, events=event_list) + if "aeff-max" in self.methods: + event_list = self.make_energy_aeff_max(observation) + observation = observation.copy(in_memory=True, events=event_list) + if "bkg-peak" in self.methods: + event_list = self.make_energy_bkg_peak(observation) + observation = observation.copy(in_memory=True, events=event_list) + if "edisp-bias" in self.methods: + event_list = self.make_energy_edisp_bias(observation) + observation = observation.copy(in_memory=True, events=event_list) + + return observation diff --git a/gammapy/makers/tests/test_safe.py b/gammapy/makers/tests/test_safe.py index e81fba69f7..0375494912 100644 --- a/gammapy/makers/tests/test_safe.py +++ b/gammapy/makers/tests/test_safe.py @@ -7,7 +7,7 @@ from regions import CircleSkyRegion from gammapy.data import DataStore from gammapy.datasets import MapDataset, SpectrumDatasetOnOff -from gammapy.makers import MapDatasetMaker, SafeMaskMaker +from gammapy.makers import MapDatasetMaker, SafeMakerDL3, SafeMaskMaker from gammapy.maps import MapAxis, RegionGeom, WcsGeom from gammapy.utils.testing import requires_data @@ -282,3 +282,29 @@ def test_safe_mask_maker_bkg_invalid(observations_hess_dl3): dataset = safe_mask_maker_nonan.run(dataset, obs) assert_allclose(dataset.mask_safe, mask_nonan) + + +@requires_data() +def test_safe_maker_dl3(observations_hess_dl3): + obs = observations_hess_dl3[0] + maker = SafeMakerDL3( + methods=["offset-max", "aeff-max", "aeff-default", "bkg-peak", "edisp-bias"] + ) + + evt = maker.make_energy_aeff_default(obs) + assert_allclose(len(evt.table), 4132) + + evt = maker.make_energy_aeff_max(obs) + assert_allclose(len(evt.table), 6903) + + evt = maker.make_energy_bkg_peak(obs) + assert_allclose(len(evt.table), 5633) + + evt = maker.make_offset_max(obs) + assert_allclose(len(evt.table), 6920) + + evt = maker.make_energy_edisp_bias(obs) + assert_allclose(len(evt.table), 3966) + + obs6 = maker.run(obs) + assert_allclose(len(obs6.events.table), 3327) From f98fedaae17f58410cad7ee575d83e3a66046bdf Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Thu, 6 Jul 2023 12:56:56 +0530 Subject: [PATCH 02/11] add docstring Signed-off-by: Atreyee Sinha --- gammapy/makers/safe.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/gammapy/makers/safe.py b/gammapy/makers/safe.py index 17fd48d23d..b7d7fbb697 100644 --- a/gammapy/makers/safe.py +++ b/gammapy/makers/safe.py @@ -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 : str, "DL3" or "DL4" + Whether to use reprojected ("DL4") or raw ("DL3") irfs + Default is "DL3" """ tag = "SafeMaskMaker" From 6f9f8f8ca950d72d2c4b748748e35d40d69babb6 Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Thu, 6 Jul 2023 21:56:36 +0530 Subject: [PATCH 03/11] add dl3 option Signed-off-by: Atreyee Sinha --- gammapy/makers/__init__.py | 3 +- gammapy/makers/safe.py | 434 +++++++++--------------------- gammapy/makers/tests/test_safe.py | 54 ++-- 3 files changed, 152 insertions(+), 339 deletions(-) diff --git a/gammapy/makers/__init__.py b/gammapy/makers/__init__.py index 0a8f6d344d..d57f3d92d2 100644 --- a/gammapy/makers/__init__.py +++ b/gammapy/makers/__init__.py @@ -13,7 +13,7 @@ from .core import Maker from .map import MapDatasetMaker from .reduce import DatasetsMaker -from .safe import SafeMakerDL3, SafeMaskMaker +from .safe import SafeMaskMaker from .spectrum import SpectrumDatasetMaker MAKER_REGISTRY = Registry( @@ -44,7 +44,6 @@ "RegionsFinder", "RingBackgroundMaker", "SafeMaskMaker", - "SafeMakerDL3", "SpectrumDatasetMaker", "WobbleRegionsFinder", ] diff --git a/gammapy/makers/safe.py b/gammapy/makers/safe.py index b7d7fbb697..bfd838d035 100644 --- a/gammapy/makers/safe.py +++ b/gammapy/makers/safe.py @@ -48,7 +48,7 @@ class SafeMaskMaker(Maker): Maximum offset cut. irfs : str, "DL3" or "DL4" Whether to use reprojected ("DL4") or raw ("DL3") irfs - Default is "DL3" + Default is "DL4" """ tag = "SafeMaskMaker" @@ -68,6 +68,7 @@ def __init__( position=None, fixed_offset=None, offset_max="3 deg", + irfs="DL4", ): methods = set(methods) @@ -86,6 +87,10 @@ def __init__( "`position` and `fixed_offset` attributes are mutually exclusive" ) + if irfs not in ["DL3", "DL4"]: + raise IOError("Invalid option fro irfs. Choose 'DL3' or 'DL4'.") + self.irfs = irfs + def make_mask_offset_max(self, dataset, observation): """Make maximum offset mask. @@ -126,7 +131,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) @@ -165,45 +170,64 @@ def make_mask_energy_aeff_max(self, dataset, observation=None): """ 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.irfs == "DL3": + 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 + + values = observation.aeff.evaluate( + offset=offset, energy_true=observation.aeff.axes["energy_true"].edges + ) + valid = observation.aeff.axes["energy_true"].edges[ + values > self.aeff_percent * np.max(values) / 100 + ] + energy_min = np.min(valid) - elif self.position: - position = self.position else: - position = geom.center_skydir - - 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") + 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}" + ) + + elif self.position: + position = self.position + else: + position = geom.center_skydir + + 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") - model = TemplateSpectralModel.from_region_map(aeff) + model = TemplateSpectralModel.from_region_map(aeff) - energy_true = model.energy - energy_min = energy_true[np.where(model.values > 0)[0][0]] - energy_max = energy_true[-1] + 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 - ) + 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) @@ -223,38 +247,54 @@ def make_mask_energy_edisp_bias(self, dataset, observation=None): 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.irfs == "DL3": + 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 + + edisp = observation.edisp.to_edisp_kernel(offset=offset) + energy_min = edisp.get_bias_energy(self.bias_percent / 100)[-1] - 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) + 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 isinstance(edisp, EDispKernelMap): + if position: + edisp = edisp.get_edisp_kernel(position=position) + else: + edisp = edisp.get_edisp_kernel(position=self.position) else: - edisp = edisp.get_edisp_kernel( - position=self.position, energy_axis=e_reco - ) - - energy_min = edisp.get_bias_energy(self.bias_percent / 100) - return geom.energy_mask(energy_min=energy_min[0]) + 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 + ) + + 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 @@ -267,6 +307,9 @@ 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 ------- @@ -274,10 +317,19 @@ def make_mask_energy_bkg_peak(dataset): Safe data range mask. """ geom = dataset._geom - background_spectrum = dataset.npred_background().get_spectrum() - idx = np.argmax(background_spectrum.data, axis=0) - energy_axis = geom.axes["energy"] - energy_min = energy_axis.edges[idx] + 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_min = bkg.axes["energy"].center[np.argmax(background_spectrum)] + + else: + background_spectrum = dataset.npred_background().get_spectrum() + 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) @staticmethod @@ -317,6 +369,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: @@ -339,250 +396,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 - - -class SafeMakerDL3(Maker): - """Make safe data range mask for a given observation, directly from the DL3 - IRFs - - - .. warning:: - - Currently some methods computing a safe energy range ("aeff-default", - "aeff-max" and "edisp-bias") determine a true energy range and apply - it to reconstructed energy, effectively neglecting the energy dispersion. - - Parameters - ---------- - methods : {"aeff-default", "aeff-max", "offset-max", "bkg-peak"} - Method to use for the safe energy range. Can be a - list with a combination of those. Resulting masks - are combined with logical `and`. "aeff-default" - uses the energy ranged specified in the DL3 data - files, if available. - aeff_percent : float - Percentage of the maximal effective area to be used - as lower energy threshold for method "aeff-max". - bias_percent : float - Percentage of the energy bias to be used as lower - energy threshold for method "edisp-bias" - position : `~astropy.coordinates.SkyCoord` - Position at which the `aeff_percent` or `bias_percent` are computed. - fixed_offset : `~astropy.coordinates.Angle` - offset, calculated from the pointing position, at which - the `aeff_percent` or `bias_percent` are computed. - If neither the position nor fixed_offset is specified, - it uses the position of the center of the map by default. - offset_max : str or `~astropy.units.Quantity` - Maximum offset cut. - """ - - tag = "SafeMakerDL3" - available_methods = { - "aeff-default", - "aeff-max", - "edisp-bias", - "offset-max", - "bkg-peak", - } - - def __init__( - self, - methods=("aeff-default",), - aeff_percent=10, - bias_percent=10, - position=None, - fixed_offset=None, - offset_max="3 deg", - ): - methods = set(methods) - - if not methods.issubset(self.available_methods): - difference = methods.difference(self.available_methods) - raise ValueError(f"{difference} is not a valid method.") - - self.methods = methods - self.aeff_percent = aeff_percent / 100.0 - self.bias_percent = bias_percent / 100.0 - self.position = position - self.fixed_offset = fixed_offset - self.offset_max = Angle(offset_max) - if self.position and self.fixed_offset: - raise ValueError( - "`position` and `fixed_offset` attributes are mutually exclusive" - ) - - @staticmethod - def make_energy_aeff_default(observation): - """Select events based on from aeff default energies. - - Parameters - ---------- - observation: `~gammapy.data.Observation` - Observation to apply cuts on. - - Returns - ------- - events: `~gammapy.data.EventList` - Event list with applied cuts. - """ - - energy_max = observation.aeff.meta.get("HI_THRES", None) - - if energy_max: - energy_max = energy_max * u.TeV - else: - log.warning( - f"No default upper safe energy threshold defined for obs {observation.obs_id}" - ) - energy_max = observation.aeff.axes["energy_true"].edges[-1] - - energy_min = observation.aeff.meta.get("LO_THRES", None) - - if energy_min: - energy_min = energy_min * u.TeV - else: - log.warning( - f"No default lower safe energy threshold defined for obs {observation.obs_id}" - ) - energy_min = observation.aeff.axes["energy_true"].edges[0] - # TODO: Apply edisp - return observation.events.select_energy([energy_min, energy_max]) - - def make_energy_aeff_max(self, observation): - """Apply selection based on a minimal effective area - - TODO: Should this be offset dependent? - - Parameters - ---------- - observation: `~gammapy.data.Observation` - Observation to apply cuts on. - - Returns - ------- - events: `~gammapy.data.EventList` - Event list with applied cuts. - """ - - 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 - - values = observation.aeff.evaluate( - offset=offset, energy_true=observation.aeff.axes["energy_true"].edges - ) - valid = observation.aeff.axes["energy_true"].edges[ - values > self.aeff_percent * np.max(values) - ] - energy_min = np.min(valid) - energy_max = observation.aeff.axes["energy_true"].edges[-1] - # TODO: Apply edisp - return observation.events.select_energy([energy_min, energy_max]) - - def make_energy_bkg_peak(self, observation): - """Apply selection based the peak of the background IRF - - Parameters - ---------- - observation: `~gammapy.data.Observation` - Observation to apply cuts on. - - Returns - ------- - events: `~gammapy.data.EventList` - Event list with applied cuts. - """ - bkg = observation.bkg.to_2d() - - values = np.ravel( - bkg.integral(axis_name="offset", offset=bkg.axes["offset"].bounds[1]) - ) - energy_min = bkg.axes["energy"].center[values == np.max(values)] - energy_max = bkg.axes["energy"].center[-1] - return observation.events.select_energy([energy_min, energy_max]) - - def make_offset_max(self, observation): - """Apply selection based on maximum offset cut - - Parameters - ---------- - observation: `~gammapy.data.Observation` - Observation to apply cuts on. - - Returns - ------- - events: `~gammapy.data.EventList` - Event list with applied cuts. - """ - return observation.events.select_offset([0 * u.deg, self.offset_max]) - - def make_energy_edisp_bias(self, observation=None): - """Apply selection based on energy dispersion bias - - Parameters - ---------- - observation: `~gammapy.data.Observation` - Observation to apply cuts on. - - Returns - ------- - events: `~gammapy.data.EventList` - Event list with applied cuts. - """ - 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 - - edisp = observation.edisp.to_edisp_kernel(offset=offset) - energy_min = edisp.get_bias_energy(self.bias_percent)[0] - energy_max = observation.edisp.axes["energy_true"].center[-1] - return observation.events.select_energy([energy_min, energy_max]) - - def run(self, observation=None): - """Make safe data range mask. - - Parameters - ---------- - dataset : `~gammapy.datasets.MapDataset` or `~gammapy.datasets.SpectrumDataset` - Dataset to compute mask for. - observation: `~gammapy.data.Observation` - Observation to compute mask for. - - Returns - ------- - observation : `~gammapy.data.Observation - Observation with selected events - """ - - if "offset-max" in self.methods: - event_list = self.make_offset_max(observation) - observation = observation.copy(in_memory=True, events=event_list) - if "aeff-default" in self.methods: - event_list = self.make_energy_aeff_default(observation) - observation = observation.copy(in_memory=True, events=event_list) - if "aeff-max" in self.methods: - event_list = self.make_energy_aeff_max(observation) - observation = observation.copy(in_memory=True, events=event_list) - if "bkg-peak" in self.methods: - event_list = self.make_energy_bkg_peak(observation) - observation = observation.copy(in_memory=True, events=event_list) - if "edisp-bias" in self.methods: - event_list = self.make_energy_edisp_bias(observation) - observation = observation.copy(in_memory=True, events=event_list) - - return observation diff --git a/gammapy/makers/tests/test_safe.py b/gammapy/makers/tests/test_safe.py index 0375494912..9bffab6470 100644 --- a/gammapy/makers/tests/test_safe.py +++ b/gammapy/makers/tests/test_safe.py @@ -7,7 +7,7 @@ from regions import CircleSkyRegion from gammapy.data import DataStore from gammapy.datasets import MapDataset, SpectrumDatasetOnOff -from gammapy.makers import MapDatasetMaker, SafeMakerDL3, SafeMaskMaker +from gammapy.makers import MapDatasetMaker, SafeMaskMaker from gammapy.maps import MapAxis, RegionGeom, WcsGeom from gammapy.utils.testing import requires_data @@ -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): @@ -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): @@ -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(), 242) + @requires_data() def test_safe_mask_maker_bkg_peak(dataset, observation_cta_1dc): @@ -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): @@ -282,29 +308,3 @@ def test_safe_mask_maker_bkg_invalid(observations_hess_dl3): dataset = safe_mask_maker_nonan.run(dataset, obs) assert_allclose(dataset.mask_safe, mask_nonan) - - -@requires_data() -def test_safe_maker_dl3(observations_hess_dl3): - obs = observations_hess_dl3[0] - maker = SafeMakerDL3( - methods=["offset-max", "aeff-max", "aeff-default", "bkg-peak", "edisp-bias"] - ) - - evt = maker.make_energy_aeff_default(obs) - assert_allclose(len(evt.table), 4132) - - evt = maker.make_energy_aeff_max(obs) - assert_allclose(len(evt.table), 6903) - - evt = maker.make_energy_bkg_peak(obs) - assert_allclose(len(evt.table), 5633) - - evt = maker.make_offset_max(obs) - assert_allclose(len(evt.table), 6920) - - evt = maker.make_energy_edisp_bias(obs) - assert_allclose(len(evt.table), 3966) - - obs6 = maker.run(obs) - assert_allclose(len(obs6.events.table), 3327) From 6e94413507fe4d4814a258d05e45be889dad3da4 Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Wed, 20 Dec 2023 17:09:16 +0100 Subject: [PATCH 04/11] simplify Signed-off-by: Quentin Remy --- gammapy/makers/safe.py | 118 ++++++++++++++++------------------------- 1 file changed, 47 insertions(+), 71 deletions(-) diff --git a/gammapy/makers/safe.py b/gammapy/makers/safe.py index bfd838d035..93b7967610 100644 --- a/gammapy/makers/safe.py +++ b/gammapy/makers/safe.py @@ -153,6 +153,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. @@ -168,17 +190,16 @@ def make_mask_energy_aeff_max(self, dataset, observation=None): mask_safe : `~numpy.ndarray` Safe data range mask. """ + + if self.fixed_offset is not None and observation is None: + raise ValueError( + f"{observation} argument is mandatory with {self.fixed_offset}" + ) + geom, exposure = dataset._geom, dataset.exposure if self.irfs == "DL3": - 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 + offset = self._get_offset(observation) values = observation.aeff.evaluate( offset=offset, energy_true=observation.aeff.axes["energy_true"].edges @@ -189,22 +210,7 @@ def make_mask_energy_aeff_max(self, dataset, observation=None): energy_min = np.min(valid) else: - 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}" - ) - - elif self.position: - position = self.position - else: - position = geom.center_skydir + position = self._get_position(observation, geom) aeff = exposure.get_spectrum(position) / exposure.meta["livetime"] if not np.any(aeff.data > 0.0): @@ -246,52 +252,24 @@ def make_mask_energy_edisp_bias(self, dataset, observation=None): mask_safe : `~numpy.ndarray` Safe data range mask. """ + + if self.fixed_offset is not None and observation is None: + raise ValueError( + f"{observation} argument is mandatory with {self.fixed_offset}" + ) + edisp, geom = dataset.edisp, dataset._geom if self.irfs == "DL3": - 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 - - edisp = observation.edisp.to_edisp_kernel(offset=offset) - energy_min = edisp.get_bias_energy(self.bias_percent / 100)[-1] - + offset = self._get_offset(observation) + edisp = observation.edisp.to_edisp_kernel(offset) else: - 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 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 - ) - - energy_min = edisp.get_bias_energy(self.bias_percent / 100)[0] + 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) def make_mask_energy_bkg_peak(self, dataset, observation=None): @@ -322,15 +300,13 @@ def make_mask_energy_bkg_peak(self, dataset, observation=None): background_spectrum = np.ravel( bkg.integral(axis_name="offset", offset=bkg.axes["offset"].bounds[1]) ) - energy_min = bkg.axes["energy"].center[np.argmax(background_spectrum)] - + energy_axis = bkg.axes["energy"] else: background_spectrum = dataset.npred_background().get_spectrum() - 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) + idx = np.argmax(background_spectrum.data, axis=0) + return geom.energy_mask(energy_min=energy_axis.edges[idx]) @staticmethod def make_mask_bkg_invalid(dataset): From ef9de8f8b9ca85c7b8ea3e08fcfede6b97f79739 Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Wed, 20 Dec 2023 17:09:28 +0100 Subject: [PATCH 05/11] adapt test Signed-off-by: Quentin Remy --- gammapy/makers/tests/test_safe.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gammapy/makers/tests/test_safe.py b/gammapy/makers/tests/test_safe.py index 9bffab6470..d042fd89a2 100644 --- a/gammapy/makers/tests/test_safe.py +++ b/gammapy/makers/tests/test_safe.py @@ -233,7 +233,7 @@ def test_safe_mask_maker_edisp_bias_fixed_offset(dataset, observation_cta_1dc): mask_edisp_bias_offset = safe_mask_maker1.make_mask_energy_edisp_bias( dataset, observation_cta_1dc ) - assert_allclose(mask_edisp_bias_offset.data.sum(), 242) + assert_allclose(mask_edisp_bias_offset.data.sum(), 1694) @requires_data() From 62aba053cb6601255766fff8844b897f4252646f Mon Sep 17 00:00:00 2001 From: Quentin Remy Date: Wed, 20 Dec 2023 17:11:08 +0100 Subject: [PATCH 06/11] Update gammapy/makers/safe.py --- gammapy/makers/safe.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gammapy/makers/safe.py b/gammapy/makers/safe.py index 93b7967610..637c959851 100644 --- a/gammapy/makers/safe.py +++ b/gammapy/makers/safe.py @@ -88,7 +88,7 @@ def __init__( ) if irfs not in ["DL3", "DL4"]: - raise IOError("Invalid option fro irfs. Choose 'DL3' or 'DL4'.") + raise IOError("Invalid option for irfs. Choose 'DL3' or 'DL4'.") self.irfs = irfs def make_mask_offset_max(self, dataset, observation): From 30f44431ee8daa7ddc5972b24ac7c563af287da6 Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Fri, 19 Jan 2024 15:17:30 +0100 Subject: [PATCH 07/11] Update gammapy/makers/safe.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Régis Terrier --- gammapy/makers/safe.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gammapy/makers/safe.py b/gammapy/makers/safe.py index 637c959851..0f3b39de1e 100644 --- a/gammapy/makers/safe.py +++ b/gammapy/makers/safe.py @@ -48,7 +48,7 @@ class SafeMaskMaker(Maker): Maximum offset cut. irfs : str, "DL3" or "DL4" Whether to use reprojected ("DL4") or raw ("DL3") irfs - Default is "DL4" + Default is "DL4". """ tag = "SafeMaskMaker" From cb250d5b90342dd04d2c6131b99c108e635c9c5d Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Fri, 19 Jan 2024 15:18:03 +0100 Subject: [PATCH 08/11] Update gammapy/makers/safe.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Régis Terrier --- gammapy/makers/safe.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gammapy/makers/safe.py b/gammapy/makers/safe.py index 0f3b39de1e..d53fc6c98e 100644 --- a/gammapy/makers/safe.py +++ b/gammapy/makers/safe.py @@ -46,7 +46,7 @@ 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 : str, "DL3" or "DL4" + irfs : {"DL4", "DL3"} Whether to use reprojected ("DL4") or raw ("DL3") irfs Default is "DL4". """ From b43313e5a6aad1fe858fd63213c6828dbb1facd3 Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Fri, 19 Jan 2024 15:18:17 +0100 Subject: [PATCH 09/11] Update gammapy/makers/safe.py MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Régis Terrier --- gammapy/makers/safe.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gammapy/makers/safe.py b/gammapy/makers/safe.py index d53fc6c98e..22aaf6d4bc 100644 --- a/gammapy/makers/safe.py +++ b/gammapy/makers/safe.py @@ -47,7 +47,7 @@ class SafeMaskMaker(Maker): offset_max : str or `~astropy.units.Quantity` Maximum offset cut. irfs : {"DL4", "DL3"} - Whether to use reprojected ("DL4") or raw ("DL3") irfs + Whether to use reprojected ("DL4") or raw ("DL3") irfs. Default is "DL4". """ From e1f422aaea01a38143995b70770d2288aa502924 Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Fri, 19 Jan 2024 15:21:14 +0100 Subject: [PATCH 10/11] Update gammapy/makers/safe.py --- gammapy/makers/safe.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gammapy/makers/safe.py b/gammapy/makers/safe.py index 22aaf6d4bc..f074ca8d6d 100644 --- a/gammapy/makers/safe.py +++ b/gammapy/makers/safe.py @@ -88,7 +88,7 @@ def __init__( ) if irfs not in ["DL3", "DL4"]: - raise IOError("Invalid option for irfs. Choose 'DL3' or 'DL4'.") + ValueError("Invalid option for irfs: expected 'DL3' or 'DL4', got {irfs} instead.") self.irfs = irfs def make_mask_offset_max(self, dataset, observation): From 8b75dc5dac0c0620a560787553486cf3b3f28b76 Mon Sep 17 00:00:00 2001 From: Atreyee Sinha Date: Fri, 19 Jan 2024 16:51:38 +0100 Subject: [PATCH 11/11] black Signed-off-by: Atreyee Sinha --- gammapy/makers/safe.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/gammapy/makers/safe.py b/gammapy/makers/safe.py index f074ca8d6d..09c381c9f1 100644 --- a/gammapy/makers/safe.py +++ b/gammapy/makers/safe.py @@ -88,7 +88,9 @@ def __init__( ) if irfs not in ["DL3", "DL4"]: - ValueError("Invalid option for irfs: expected 'DL3' or 'DL4', got {irfs} instead.") + ValueError( + "Invalid option for irfs: expected 'DL3' or 'DL4', got {irfs} instead." + ) self.irfs = irfs def make_mask_offset_max(self, dataset, observation):