From 6c4206bdc1a412dbaca95f09b0a77095380c280b Mon Sep 17 00:00:00 2001 From: Axel Donath Date: Fri, 18 Oct 2019 11:46:40 +0200 Subject: [PATCH 1/6] Implement ReflectedRegionsBackgroundMaker --- gammapy/spectrum/reflected.py | 123 +++++++++++++++++++++++++++++++++- 1 file changed, 122 insertions(+), 1 deletion(-) diff --git a/gammapy/spectrum/reflected.py b/gammapy/spectrum/reflected.py index 4c75088bb1..f3e4d40072 100644 --- a/gammapy/spectrum/reflected.py +++ b/gammapy/spectrum/reflected.py @@ -7,8 +7,10 @@ from gammapy.maps import Map, WcsGeom, WcsNDMap from gammapy.maps.geom import frame_to_coordsys from .background_estimate import BackgroundEstimate +from .dataset import SpectrumDatasetOnOff -__all__ = ["ReflectedRegionsFinder", "ReflectedRegionsBackgroundEstimator"] + +__all__ = ["ReflectedRegionsFinder", "ReflectedRegionsBackgroundEstimator", "ReflectedRegionsBackgroundMaker"] log = logging.getLogger(__name__) @@ -414,3 +416,122 @@ def plot(self, fig=None, ax=None, cmap=None, idx=None, add_legend=False): ax.legend(handles=handles) return fig, ax, cbar + + +class ReflectedRegionsBackgroundMaker: + """Reflected regions background maker. + + Parameters + ---------- + region: `~regions.SkyRegion` + On region to compute off regions for. + angle_increment : `~astropy.coordinates.Angle`, optional + Rotation angle applied when a region falls in an excluded region. + min_distance : `~astropy.coordinates.Angle`, optional + Minimal distance between two consecutive reflected regions + min_distance_input : `~astropy.coordinates.Angle`, optional + Minimal distance from input region + max_region_number : int, optional + Maximum number of regions to use + exclusion_mask : `~gammapy.maps.WcsNDMap`, optional + Exclusion mask + binsz : `~astropy.coordinates.Angle` + Bin size of the reference map used for region finding. Default : 0.01 deg + """ + def __init__( + self, + region, + angle_increment="0.1 rad", + min_distance="0 rad", + min_distance_input="0.1 rad", + max_region_number=10000, + exclusion_mask=None, + binsz="0.01 deg", + ): + self.region = region + self.binsz = binsz + self.exclusion_mask = exclusion_mask + self.angle_increment = Angle(angle_increment) + self.min_distance = Angle(min_distance) + self.min_distance_input = Angle(min_distance_input) + self.max_region_number = max_region_number + + def _get_finder(self, observation): + finder = ReflectedRegionsFinder( + binsz=self.binsz, + exclusion_mask=self.exclusion_mask, + center=observation.pointing_radec, + region=self.region, + min_distance=self.min_distance, + min_distance_input=self.min_distance_input, + max_region_number=self.max_region_number, + angle_increment=self.angle_increment, + ) + return finder + + def make_counts_off(self, dataset, observation): + """Make off counts. + + Parameters + ---------- + dataset : `SpectrumDataset` + Spectrum dataset. + observation : `DatastoreObservation` + Data store observation. + + + Returns + ------- + counts_off : `CountsSpectrum` + Off counts. + """ + finder = self._get_finder(observation) + + finder.run() + regions = finder.reflected_regions + + region_union = regions[0] + + for region in regions[1:]: + region_union = region_union.union(region) + + wcs = finder.reference_map.geom.wcs + events_off = observation.events.select_region(region_union, wcs) + + counts_off = dataset.counts.copy() + counts_off.fill_events(events_off) + counts_off.regions = regions + return counts_off + + def run(self, dataset, observation): + """Run reflected regions background maker + + Parameters + ---------- + dataset : `SpectrumDataset` + Spectrum dataset. + observation : `DatastoreObservation` + Data store observation. + + Returns + ------- + dataset_on_off : `SpectrumDatasetOnOff` + On off dataset. + + """ + kwargs = {} + + kwargs["counts"] = dataset.counts + kwargs["gti"] = dataset.gti + kwargs["name"] = dataset.name + kwargs["mask_safe"] = dataset.mask_safe + kwargs["mask_fit"] = dataset.mask_fit + kwargs["aeff"] = dataset.aeff + kwargs["livetime"] = dataset.livetime + kwargs["edisp"] = dataset.edisp + + counts_off = self.make_counts_off(dataset, observation) + kwargs["acceptance"] = 1 + kwargs["acceptance_off"] = len(counts_off.regions) + kwargs["counts_off"] = counts_off + return SpectrumDatasetOnOff(**kwargs) From 9cb05209f25d2fac6b677903aaa38de0403a6a35 Mon Sep 17 00:00:00 2001 From: Axel Donath Date: Fri, 18 Oct 2019 11:47:09 +0200 Subject: [PATCH 2/6] Adapt tests in test_reflected.py --- gammapy/spectrum/tests/test_reflected.py | 54 ++++++++++++------------ 1 file changed, 28 insertions(+), 26 deletions(-) diff --git a/gammapy/spectrum/tests/test_reflected.py b/gammapy/spectrum/tests/test_reflected.py index 16c55ca259..967096962a 100644 --- a/gammapy/spectrum/tests/test_reflected.py +++ b/gammapy/spectrum/tests/test_reflected.py @@ -1,7 +1,9 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst import pytest +import numpy as np import astropy.units as u from astropy.coordinates import Angle, SkyCoord +from numpy.testing import assert_allclose from regions import ( CircleSkyRegion, EllipseAnnulusSkyRegion, @@ -10,13 +12,13 @@ ) from gammapy.data import DataStore from gammapy.maps import WcsGeom, WcsNDMap -from gammapy.spectrum import ReflectedRegionsBackgroundEstimator, ReflectedRegionsFinder +from gammapy.spectrum import ReflectedRegionsFinder, ReflectedRegionsBackgroundMaker from gammapy.utils.testing import ( assert_quantity_allclose, mpl_plot_check, requires_data, - requires_dependency, ) +from gammapy.spectrum.make import SpectrumDatasetMaker @pytest.fixture(scope="session") @@ -46,19 +48,17 @@ def observations(): return datastore.get_observations(obs_ids) -@pytest.fixture(scope="session") -def bkg_estimator(observations, exclusion_mask, on_region): - """Example background estimator for testing.""" - maker = ReflectedRegionsBackgroundEstimator( - observations=observations, - on_region=on_region, - exclusion_mask=exclusion_mask, - min_distance_input="0.2 deg", - ) - maker.run() - return maker +@pytest.fixture() +def spectrum_dataset_maker(on_region): + e_reco = np.logspace(0, 2, 5) * u.TeV + e_true = np.logspace(-0.5, 2, 11) * u.TeV + return SpectrumDatasetMaker(region=on_region, e_reco=e_reco, e_true=e_true) +@pytest.fixture() +def reflected_bkg_maker(on_region, exclusion_mask): + return ReflectedRegionsBackgroundMaker(region=on_region, exclusion_mask=exclusion_mask) + region_finder_param = [ (SkyCoord(83.2, 22.5, unit="deg"), 15, Angle("82.592 deg"), 17, 17), (SkyCoord(84.2, 22.5, unit="deg"), 17, Angle("83.636 deg"), 19, 19), @@ -161,16 +161,18 @@ def bad_on_region(exclusion_mask, on_region): @requires_data() -class TestReflectedRegionBackgroundEstimator: - def test_basic(self, bkg_estimator): - assert "ReflectedRegionsBackgroundEstimator" in str(bkg_estimator) - - def test_run(self, bkg_estimator): - assert len(bkg_estimator.result[1].off_region) == 11 - assert "Reflected" in str(bkg_estimator.result[1]) - - @requires_dependency("matplotlib") - def test_plot(self, bkg_estimator): - with mpl_plot_check(): - bkg_estimator.plot(idx=1, add_legend=True) - bkg_estimator.plot(idx=[0, 1]) +def test_reflected_bkg_maker(spectrum_dataset_maker, reflected_bkg_maker, observations): + + datasets = [] + + for obs in observations: + dataset = spectrum_dataset_maker.run(obs, selection=["counts"]) + dataset_on_off = reflected_bkg_maker.run(dataset, obs) + datasets.append(dataset_on_off) + + assert_allclose(datasets[0].counts_off.data.sum(), 76) + assert_allclose(datasets[1].counts_off.data.sum(), 60) + + assert_allclose(len(datasets[0].counts_off.regions), 11) + assert_allclose(len(datasets[1].counts_off.regions), 11) + From 1744c3aee68f0cc4e39ff2ecb0e98df449be6efc Mon Sep 17 00:00:00 2001 From: Axel Donath Date: Fri, 18 Oct 2019 11:51:40 +0200 Subject: [PATCH 3/6] Polish and add TODOs --- gammapy/spectrum/reflected.py | 11 +++++++++-- gammapy/spectrum/tests/test_reflected.py | 12 +++++++----- 2 files changed, 16 insertions(+), 7 deletions(-) diff --git a/gammapy/spectrum/reflected.py b/gammapy/spectrum/reflected.py index f3e4d40072..45ffc0da05 100644 --- a/gammapy/spectrum/reflected.py +++ b/gammapy/spectrum/reflected.py @@ -9,8 +9,11 @@ from .background_estimate import BackgroundEstimate from .dataset import SpectrumDatasetOnOff - -__all__ = ["ReflectedRegionsFinder", "ReflectedRegionsBackgroundEstimator", "ReflectedRegionsBackgroundMaker"] +__all__ = [ + "ReflectedRegionsFinder", + "ReflectedRegionsBackgroundEstimator", + "ReflectedRegionsBackgroundMaker", +] log = logging.getLogger(__name__) @@ -438,6 +441,7 @@ class ReflectedRegionsBackgroundMaker: binsz : `~astropy.coordinates.Angle` Bin size of the reference map used for region finding. Default : 0.01 deg """ + def __init__( self, region, @@ -500,6 +504,9 @@ def make_counts_off(self, dataset, observation): counts_off = dataset.counts.copy() counts_off.fill_events(events_off) + + # TODO: temporarily attach the region info to the CountsSpectrum. + # but CountsSpectrum should be modified to hold this info anyway. counts_off.regions = regions return counts_off diff --git a/gammapy/spectrum/tests/test_reflected.py b/gammapy/spectrum/tests/test_reflected.py index 967096962a..12f989128d 100644 --- a/gammapy/spectrum/tests/test_reflected.py +++ b/gammapy/spectrum/tests/test_reflected.py @@ -1,9 +1,9 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst import pytest import numpy as np +from numpy.testing import assert_allclose import astropy.units as u from astropy.coordinates import Angle, SkyCoord -from numpy.testing import assert_allclose from regions import ( CircleSkyRegion, EllipseAnnulusSkyRegion, @@ -12,13 +12,13 @@ ) from gammapy.data import DataStore from gammapy.maps import WcsGeom, WcsNDMap -from gammapy.spectrum import ReflectedRegionsFinder, ReflectedRegionsBackgroundMaker +from gammapy.spectrum import ReflectedRegionsBackgroundMaker, ReflectedRegionsFinder +from gammapy.spectrum.make import SpectrumDatasetMaker from gammapy.utils.testing import ( assert_quantity_allclose, mpl_plot_check, requires_data, ) -from gammapy.spectrum.make import SpectrumDatasetMaker @pytest.fixture(scope="session") @@ -57,7 +57,10 @@ def spectrum_dataset_maker(on_region): @pytest.fixture() def reflected_bkg_maker(on_region, exclusion_mask): - return ReflectedRegionsBackgroundMaker(region=on_region, exclusion_mask=exclusion_mask) + return ReflectedRegionsBackgroundMaker( + region=on_region, exclusion_mask=exclusion_mask + ) + region_finder_param = [ (SkyCoord(83.2, 22.5, unit="deg"), 15, Angle("82.592 deg"), 17, 17), @@ -175,4 +178,3 @@ def test_reflected_bkg_maker(spectrum_dataset_maker, reflected_bkg_maker, observ assert_allclose(len(datasets[0].counts_off.regions), 11) assert_allclose(len(datasets[1].counts_off.regions), 11) - From 76e93ea3165cb4d73310f090f79923e466ae11d8 Mon Sep 17 00:00:00 2001 From: Axel Donath Date: Fri, 18 Oct 2019 12:06:27 +0200 Subject: [PATCH 4/6] Change coding pattern in ReflectedRegionsBackgroundMaker.run() --- gammapy/spectrum/reflected.py | 28 +++++++++++++--------------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/gammapy/spectrum/reflected.py b/gammapy/spectrum/reflected.py index 45ffc0da05..86f51a7494 100644 --- a/gammapy/spectrum/reflected.py +++ b/gammapy/spectrum/reflected.py @@ -526,19 +526,17 @@ def run(self, dataset, observation): On off dataset. """ - kwargs = {} - - kwargs["counts"] = dataset.counts - kwargs["gti"] = dataset.gti - kwargs["name"] = dataset.name - kwargs["mask_safe"] = dataset.mask_safe - kwargs["mask_fit"] = dataset.mask_fit - kwargs["aeff"] = dataset.aeff - kwargs["livetime"] = dataset.livetime - kwargs["edisp"] = dataset.edisp - counts_off = self.make_counts_off(dataset, observation) - kwargs["acceptance"] = 1 - kwargs["acceptance_off"] = len(counts_off.regions) - kwargs["counts_off"] = counts_off - return SpectrumDatasetOnOff(**kwargs) + + return SpectrumDatasetOnOff( + counts=dataset.counts, + counts_off=counts_off, + gti=dataset.gti, + name=dataset.name, + livetime=dataset.livetime, + edisp=dataset.edisp, + aeff=dataset.aeff, + acceptance_off=len(counts_off.regions), + mask_safe=dataset.mask_safe, + mask_fit=dataset.mask_fit, + ) From ccf8713651897007bdbb64ba537fc7599228df08 Mon Sep 17 00:00:00 2001 From: Axel Donath Date: Fri, 18 Oct 2019 14:20:40 +0200 Subject: [PATCH 5/6] Remove default values from docstrings in reflected.py --- gammapy/spectrum/reflected.py | 18 ++++++++---------- gammapy/spectrum/tests/test_reflected.py | 1 - 2 files changed, 8 insertions(+), 11 deletions(-) diff --git a/gammapy/spectrum/reflected.py b/gammapy/spectrum/reflected.py index 86f51a7494..620ddaeb50 100644 --- a/gammapy/spectrum/reflected.py +++ b/gammapy/spectrum/reflected.py @@ -49,7 +49,7 @@ class ReflectedRegionsFinder: exclusion_mask : `~gammapy.maps.WcsNDMap`, optional Exclusion mask binsz : `~astropy.coordinates.Angle` - Bin size of the reference map used for region finding. Default : 0.01 deg + Bin size of the reference map used for region finding. Examples -------- @@ -137,9 +137,9 @@ def make_reference_map(region, center, binsz="0.01 deg", min_width="0.3 deg"): center : `~astropy.coordinates.SkyCoord` Rotation point binsz : `~astropy.coordinates.Angle` - Reference map bin size. Default : 0.01 deg + Reference map bin size. min_width : `~astropy.coordinates.Angle` - Minimal map width. Default : 0.3 deg + Minimal map width. Returns ------- @@ -270,7 +270,7 @@ class ReflectedRegionsBackgroundEstimator: observations : `~gammapy.data.Observations` Observations to process binsz : `~astropy.coordinates.Angle` - Optional, bin size of the maps used to compute the regions, Default '0.01 deg' + Optional, bin size of the maps used to compute the regions. kwargs : dict Forwarded to `~gammapy.spectrum.ReflectedRegionsFinder` """ @@ -350,9 +350,9 @@ def plot(self, fig=None, ax=None, cmap=None, idx=None, add_legend=False): cmap : `~matplotlib.colors.ListedColormap`, optional Color map to use idx : int, optional - Observations to include in the plot, default: all + Observations to include in the plot. add_legend : boolean, optional - Enable/disable legend in the plot, default: False + Enable/disable legend in the plot. """ import matplotlib.pyplot as plt @@ -439,7 +439,7 @@ class ReflectedRegionsBackgroundMaker: exclusion_mask : `~gammapy.maps.WcsNDMap`, optional Exclusion mask binsz : `~astropy.coordinates.Angle` - Bin size of the reference map used for region finding. Default : 0.01 deg + Bin size of the reference map used for region finding. """ def __init__( @@ -461,7 +461,7 @@ def __init__( self.max_region_number = max_region_number def _get_finder(self, observation): - finder = ReflectedRegionsFinder( + return ReflectedRegionsFinder( binsz=self.binsz, exclusion_mask=self.exclusion_mask, center=observation.pointing_radec, @@ -471,7 +471,6 @@ def _get_finder(self, observation): max_region_number=self.max_region_number, angle_increment=self.angle_increment, ) - return finder def make_counts_off(self, dataset, observation): """Make off counts. @@ -524,7 +523,6 @@ def run(self, dataset, observation): ------- dataset_on_off : `SpectrumDatasetOnOff` On off dataset. - """ counts_off = self.make_counts_off(dataset, observation) diff --git a/gammapy/spectrum/tests/test_reflected.py b/gammapy/spectrum/tests/test_reflected.py index 12f989128d..4dea3108ca 100644 --- a/gammapy/spectrum/tests/test_reflected.py +++ b/gammapy/spectrum/tests/test_reflected.py @@ -165,7 +165,6 @@ def bad_on_region(exclusion_mask, on_region): @requires_data() def test_reflected_bkg_maker(spectrum_dataset_maker, reflected_bkg_maker, observations): - datasets = [] for obs in observations: From c0743870fcebf89456c21d43e65d0129dfd81bcc Mon Sep 17 00:00:00 2001 From: Axel Donath Date: Fri, 18 Oct 2019 16:37:30 +0200 Subject: [PATCH 6/6] Don't use counts.copy() to create counts_off --- gammapy/spectrum/reflected.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/gammapy/spectrum/reflected.py b/gammapy/spectrum/reflected.py index 620ddaeb50..5d7a2b37e5 100644 --- a/gammapy/spectrum/reflected.py +++ b/gammapy/spectrum/reflected.py @@ -8,6 +8,7 @@ from gammapy.maps.geom import frame_to_coordsys from .background_estimate import BackgroundEstimate from .dataset import SpectrumDatasetOnOff +from .core import CountsSpectrum __all__ = [ "ReflectedRegionsFinder", @@ -501,7 +502,11 @@ def make_counts_off(self, dataset, observation): wcs = finder.reference_map.geom.wcs events_off = observation.events.select_region(region_union, wcs) - counts_off = dataset.counts.copy() + edges = dataset.counts.energy.edges + counts_off = CountsSpectrum( + energy_hi=edges[1:], + energy_lo=edges[:-1], + ) counts_off.fill_events(events_off) # TODO: temporarily attach the region info to the CountsSpectrum.