From 2227729af44e58ab5055d0f8ff862a10e762a605 Mon Sep 17 00:00:00 2001 From: AtreyeeS Date: Fri, 18 Oct 2019 20:01:39 +0200 Subject: [PATCH 01/11] Adapt MapDataset --- gammapy/cube/fit.py | 81 ++++++++++++++++++---------------- gammapy/cube/tests/test_fit.py | 28 +++++++----- 2 files changed, 61 insertions(+), 48 deletions(-) diff --git a/gammapy/cube/fit.py b/gammapy/cube/fit.py index 77f29b3d53..c135fabc98 100644 --- a/gammapy/cube/fit.py +++ b/gammapy/cube/fit.py @@ -12,7 +12,7 @@ from gammapy.cube.psf_map import PSFMap from gammapy.data import GTI from gammapy.irf import EffectiveAreaTable, EnergyDispersion, apply_containment_fraction -from gammapy.maps import Map, MapAxis +from gammapy.maps import Map, MapAxis, WcsGeom from gammapy.modeling import Dataset, Parameters from gammapy.modeling.models import BackgroundModel, SkyModel, SkyModels from gammapy.spectrum import SpectrumDataset @@ -33,7 +33,12 @@ MIGRA_AXIS_DEFAULT = MapAxis.from_bounds( 0.2, 5, nbin=48, node_type="edges", name="migra" ) +ENERGY_AXIS_DEFAULT = MapAxis.from_edges( + np.logspace(-2.0, 2.5, 109), name="energy", unit=u.TeV, interp="log" +) BINSZ_IRF = 0.2 +BINSZ_IRF_DEFAULT = BINSZ_IRF +MARGIN_IRF_DEFAULT = 0.5 # TODO: Choose optimal binnings depending on IRFs @@ -363,9 +368,11 @@ def from_geoms( def create( cls, geom, - geom_irf=None, + energy_axis_true=None, migra_axis=None, rad_axis=None, + binsz_irf=None, + margin_irf=None, reference_time="2000-01-01", name="", **kwargs, @@ -376,54 +383,54 @@ def create( ---------- geom: `~gammapy.maps.WcsGeom` Reference target geometry in reco energy, used for counts and background maps - geom_irf: `~gammapy.maps.WcsGeom` - Reference image geometry in true energy, used for IRF maps. + energy_axis_true: `~gammapy.maps.MapAxis` + True energy axis used for IRF maps migra_axis: `~gammapy.maps.MapAxis` Migration axis for the energy dispersion map rad_axis: `~gammapy.maps.MapAxis` Rad axis for the psf map + binsz_irf_maps: float + IRF Map pixel size in degrees. + margin_irf: float + IRF map margin size in degrees reference_time: `~astropy.time.Time` the reference time to use in GTI definition name : str Name of the dataset. + + Returns + -------- + empty_maps: `MapDataset` + A MapDataset containing zero filled maps """ - geom_irf = geom_irf or geom.to_binsz(BINSZ_IRF) + migra_axis = migra_axis or MIGRA_AXIS_DEFAULT rad_axis = rad_axis or RAD_AXIS_DEFAULT + energy_axis_true = energy_axis_true or ENERGY_AXIS_DEFAULT + binsz_irf = binsz_irf or BINSZ_IRF_DEFAULT + margin_irf = margin_irf or MARGIN_IRF_DEFAULT + margin_irf = margin_irf * u.deg + + wcs = geom.to_image() + geom_exposure = wcs.to_cube([energy_axis_true]) + + wcs_irf = WcsGeom.create( + binsz=binsz_irf, + width=wcs.width + margin_irf, + skydir=wcs.center_skydir, + proj=wcs.projection, + coordsys=wcs.coordsys, + ) - counts = Map.from_geom(geom, unit="") - - background = Map.from_geom(geom, unit="") - background_model = BackgroundModel(background) - - energy_axis = geom_irf.get_axis_by_name("ENERGY") - - exposure_geom = geom.to_image().to_cube([energy_axis]) - exposure = Map.from_geom(exposure_geom, unit="m2 s") - exposure_irf = Map.from_geom(geom_irf, unit="m2 s") - - mask_safe = np.zeros(geom.data_shape, dtype=bool) - - gti = GTI.create([] * u.s, [] * u.s, reference_time=reference_time) - - geom_migra = geom_irf.to_image().to_cube([migra_axis, energy_axis]) - edisp_map = Map.from_geom(geom_migra, unit="") - loc = migra_axis.edges.searchsorted(1.0) - edisp_map.data[:, loc, :, :] = 1.0 - edisp = EDispMap(edisp_map, exposure_irf) - - geom_rad = geom_irf.to_image().to_cube([rad_axis, energy_axis]) - psf_map = Map.from_geom(geom_rad, unit="sr-1") - psf = PSFMap(psf_map, exposure_irf) + geom_psf = wcs_irf.to_cube([rad_axis, energy_axis_true]) + geom_edisp = wcs_irf.to_cube([migra_axis, energy_axis_true]) - return cls( - counts=counts, - exposure=exposure, - psf=psf, - edisp=edisp, - background_model=background_model, - gti=gti, - mask_safe=mask_safe, + return cls.from_geoms( + geom, + geom_exposure, + geom_psf, + geom_edisp, + reference_time=reference_time, name=name, **kwargs, ) diff --git a/gammapy/cube/tests/test_fit.py b/gammapy/cube/tests/test_fit.py index f827801cd2..2cb6635e93 100644 --- a/gammapy/cube/tests/test_fit.py +++ b/gammapy/cube/tests/test_fit.py @@ -360,22 +360,29 @@ def test_map_fit_one_energy_bin(sky_model, geom_image): def test_create(geom, geom_etrue): # tests empty datasets created - migra_axis = MapAxis(nodes=np.linspace(0.0, 3.0, 51), unit="", name="migra") rad_axis = MapAxis(nodes=np.linspace(0.0, 1.0, 51), unit="deg", name="theta") + e_reco = MapAxis.from_edges( + np.logspace(-1.0, 1.0, 3), name="energy", unit=u.TeV, interp="log" + ) + e_true = MapAxis.from_edges( + np.logspace(-1.0, 1.0, 4), name="energy", unit=u.TeV, interp="log" + ) + geom = WcsGeom.create(binsz=0.02, width=(2, 2), axes=[e_reco]) + empty_dataset = MapDataset.create( + geom=geom, energy_axis_true=e_true, migra_axis=migra_axis, rad_axis=rad_axis + ) - empty_dataset = MapDataset.create(geom, geom_etrue, migra_axis, rad_axis) - - assert_allclose(empty_dataset.counts.data.sum(), 0.0) - assert_allclose(empty_dataset.background_model.map.data.sum(), 0.0) + assert empty_dataset.counts.data.shape == (2, 100, 100) - assert empty_dataset.psf.psf_map.data.shape == (3, 50, 100, 100) - assert empty_dataset.psf.exposure_map.data.shape == (3, 1, 100, 100) + assert empty_dataset.exposure.data.shape == (3, 100, 100) - assert empty_dataset.edisp.edisp_map.data.shape == (3, 50, 100, 100) - assert empty_dataset.edisp.exposure_map.data.shape == (3, 1, 100, 100) + assert empty_dataset.psf.psf_map.data.shape == (3, 50, 12, 12) + assert empty_dataset.psf.exposure_map.data.shape == (3, 1, 12, 12) - assert_allclose(empty_dataset.edisp.edisp_map.data.sum(), 30000) + assert empty_dataset.edisp.edisp_map.data.shape == (3, 50, 12, 12) + assert empty_dataset.edisp.exposure_map.data.shape == (3, 1, 12, 12) + assert_allclose(empty_dataset.edisp.edisp_map.data.sum(), 432) assert_allclose(empty_dataset.gti.time_delta, 0.0 * u.s) @@ -394,7 +401,6 @@ def test_from_geoms(): wcs_irf = WcsGeom.create(binsz=0.1, width=(2.5, 2.5)) geom = wcs.to_cube([e_reco]) geom_exposure = wcs.to_cube([e_true]) - geom_irf = wcs_irf.to_cube([e_true]) geom_psf = wcs_irf.to_cube([rad_axis, e_true]) geom_edisp = wcs_irf.to_cube([migra_axis, e_true]) From 18447e7c60f7ac1dcc1ede650f52278c00be8950 Mon Sep 17 00:00:00 2001 From: AtreyeeS Date: Sat, 19 Oct 2019 01:45:57 +0200 Subject: [PATCH 02/11] Adapt MapDatasetMaker --- gammapy/cube/fit.py | 4 +- gammapy/cube/make.py | 65 +++++++++++++++++++++++++-------- gammapy/cube/tests/test_make.py | 49 +++++++++++++++++-------- 3 files changed, 84 insertions(+), 34 deletions(-) diff --git a/gammapy/cube/fit.py b/gammapy/cube/fit.py index c135fabc98..e247e208bd 100644 --- a/gammapy/cube/fit.py +++ b/gammapy/cube/fit.py @@ -36,8 +36,8 @@ ENERGY_AXIS_DEFAULT = MapAxis.from_edges( np.logspace(-2.0, 2.5, 109), name="energy", unit=u.TeV, interp="log" ) -BINSZ_IRF = 0.2 -BINSZ_IRF_DEFAULT = BINSZ_IRF + +BINSZ_IRF_DEFAULT = 0.2 MARGIN_IRF_DEFAULT = 0.5 # TODO: Choose optimal binnings depending on IRFs diff --git a/gammapy/cube/make.py b/gammapy/cube/make.py index 49cb491735..77b8677a46 100644 --- a/gammapy/cube/make.py +++ b/gammapy/cube/make.py @@ -2,16 +2,24 @@ import logging from functools import lru_cache import numpy as np +import astropy.units as u from astropy.coordinates import Angle from astropy.nddata.utils import NoOverlapError from astropy.utils import lazyproperty from gammapy.irf import EnergyDependentMultiGaussPSF -from gammapy.maps import Map +from gammapy.maps import Map, WcsGeom from gammapy.modeling.models import BackgroundModel from .background import make_map_background_irf from .edisp_map import make_edisp_map from .exposure import _map_spectrum_weight, make_map_exposure_true_energy -from .fit import BINSZ_IRF, MIGRA_AXIS_DEFAULT, RAD_AXIS_DEFAULT, MapDataset +from .fit import ( + BINSZ_IRF_DEFAULT, + MARGIN_IRF_DEFAULT, + MIGRA_AXIS_DEFAULT, + RAD_AXIS_DEFAULT, + ENERGY_AXIS_DEFAULT, + MapDataset +) from .psf_map import make_psf_map __all__ = ["MapDatasetMaker", "MapMakerRing"] @@ -28,14 +36,16 @@ class MapDatasetMaker: Reference image geometry in reco energy, used for counts and background maps offset_max : `~astropy.coordinates.Angle` Maximum offset angle - geom_true : `~gammapy.maps.WcsGeom` - Reference image geometry in true energy, used for IRF maps. It can have a coarser - spatial bins than the counts geom. - If none, the same as geom is assumed + energy_axis_true: `~gammapy.maps.MapAxis` + True energy axis used for IRF maps migra_axis : `~gammapy.maps.MapAxis` Migration axis for edisp map rad_axis : `~gammapy.maps.MapAxis` Radial axis for psf map. + binsz_irf: float + IRF Map pixel size in degrees. + margin_irf: float + IRF map margin size in degrees cutout : bool Whether to cutout the observation. cutout_mode : {'trim', 'partial', 'strict'} @@ -47,19 +57,27 @@ def __init__( self, geom, offset_max, - geom_true=None, background_oversampling=None, + energy_axis_true=None, migra_axis=None, rad_axis=None, + binsz_irf=None, + margin_irf=None, cutout_mode="trim", cutout=True, ): + self.geom = geom - self.geom_true = geom_true if geom_true else geom.to_binsz(BINSZ_IRF) self.offset_max = Angle(offset_max) self.background_oversampling = background_oversampling self.migra_axis = migra_axis if migra_axis else MIGRA_AXIS_DEFAULT self.rad_axis = rad_axis if rad_axis else RAD_AXIS_DEFAULT + self.energy_axis_true = energy_axis_true or ENERGY_AXIS_DEFAULT + self.binsz_irf = binsz_irf or BINSZ_IRF_DEFAULT + + self.margin_irf = margin_irf or MARGIN_IRF_DEFAULT + self.margin_irf = self.margin_irf * u.deg + self.cutout_mode = cutout_mode self.cutout_width = 2 * self.offset_max self.cutout = cutout @@ -75,25 +93,39 @@ def _cutout_geom(self, geom, observation): else: return geom + @lazyproperty + def wcs_irf(self): + """Spatial geometry of IRF Maps (`Geom`)""" + wcs = self.geom.to_image() + return WcsGeom.create( + binsz=self.binsz_irf, + width=wcs.width + self.margin_irf, + skydir=wcs.center_skydir, + proj=wcs.projection, + coordsys=wcs.coordsys, + ) + + @lazyproperty + def geom_irf(self): + """Geom of Exposure map associated with IRFs (`Geom`)""" + return self.wcs_irf.to_cube([self.energy_axis_true]) + @lazyproperty def geom_exposure(self): """Exposure map geom (`Geom`)""" - energy_axis = self.geom_true.get_axis_by_name("energy") - geom_exposure = self.geom.to_image().to_cube([energy_axis]) + geom_exposure = self.geom.to_image().to_cube([self.energy_axis_true]) return geom_exposure @lazyproperty def geom_psf(self): """PSFMap geom (`Geom`)""" - energy_axis = self.geom_true.get_axis_by_name("ENERGY") - geom_psf = self.geom_true.to_image().to_cube([self.rad_axis, energy_axis]) + geom_psf = self.wcs_irf.to_cube([self.rad_axis, self.energy_axis_true]) return geom_psf @lazyproperty def geom_edisp(self): """EdispMap geom (`Geom`)""" - energy_axis = self.geom_true.get_axis_by_name("ENERGY") - geom_edisp = self.geom_true.to_image().to_cube([self.migra_axis, energy_axis]) + geom_edisp = self.wcs_irf.to_cube([self.migra_axis, self.energy_axis_true]) return geom_edisp def make_counts(self, observation): @@ -150,7 +182,8 @@ def make_exposure_irf(self, observation): exposure : `Map` Exposure map. """ - geom = self._cutout_geom(self.geom_true, observation) + + geom = self._cutout_geom(self.geom_irf, observation) exposure = make_map_exposure_true_energy( pointing=observation.pointing_radec, livetime=observation.observation_live_time_duration, @@ -281,7 +314,7 @@ def make_mask_safe_irf(self, observation): mask : `~numpy.ndarray` Mask """ - geom = self._cutout_geom(self.geom_true, observation) + geom = self._cutout_geom(self.geom_irf, observation) offset = geom.separation(observation.pointing_radec) return offset >= self.offset_max diff --git a/gammapy/cube/tests/test_make.py b/gammapy/cube/tests/test_make.py index 5e3d48d141..f43fd5f1b6 100644 --- a/gammapy/cube/tests/test_make.py +++ b/gammapy/cube/tests/test_make.py @@ -36,64 +36,76 @@ def geom(ebounds, binsz=0.5): "pars", [ { - # Default, normal test case + # Default, same e_true and reco "geom": geom(ebounds=[0.1, 1, 10]), - "geom_true": None, + "e_true": MapAxis.from_edges([0.1, 1, 10], name="energy", unit="TeV", interp="log"), "counts": 34366, "exposure": 9.995376e08, "exposure_image": 7.921993e10, "background": 27989.05, + "binsz_irf": 0.5, + "margin_irf": 0.0, }, { # Test single energy bin "geom": geom(ebounds=[0.1, 10]), - "geom_true": None, + "e_true": MapAxis.from_edges([0.1, 10], name="energy", unit="TeV", interp="log"), "counts": 34366, "exposure": 5.843302e08, "exposure_image": 1.16866e11, "background": 30424.451, + "binsz_irf": 0.5, + "margin_irf": 0.0, }, { # Test single energy bin with exclusion mask "geom": geom(ebounds=[0.1, 10]), - "geom_true": None, + "e_true": MapAxis.from_edges([0.1, 10], name="energy", unit="TeV", interp="log"), "exclusion_mask": Map.from_geom(geom(ebounds=[0.1, 10])), "counts": 34366, "exposure": 5.843302e08, "exposure_image": 1.16866e11, "background": 30424.451, + "binsz_irf": 0.5, + "margin_irf": 0.0, }, { # Test for different e_true and e_reco bins "geom": geom(ebounds=[0.1, 1, 10]), - "geom_true": geom(ebounds=[0.1, 0.5, 2.5, 10.0]), + "e_true": MapAxis.from_edges([0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log"), "counts": 34366, "exposure": 9.951827e08, "exposure_image": 6.492968e10, "background": 28760.283, "background_oversampling": 2, + "binsz_irf": 0.5, + "margin_irf": 0.0, }, { # Test for different e_true and e_reco bins "geom": geom(ebounds=[0.1, 1, 10]), - "geom_true": geom(ebounds=[0.1, 0.5, 2.5, 10.0], binsz=1), + "e_true": MapAxis.from_edges([0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log"), "counts": 34366, "exposure": 9.951827e08, "exposure_image": 6.492968e10, "background": 28760.283, "background_oversampling": 2, + "binsz_irf": 1.0, + "margin_irf": 0.0, }, ], ) -@pytest.mark.parametrize("keepdims", [True, False]) -def test_map_maker(pars, observations, keepdims): +def test_map_maker(pars, observations): - stacked = MapDataset.create(geom=pars["geom"], geom_irf=pars["geom_true"]) + stacked = MapDataset.create(geom=pars["geom"], + energy_axis_true=pars["e_true"], + binsz_irf=pars["binsz_irf"], + margin_irf=pars["margin_irf"]) for obs in observations: maker = MapDatasetMaker( geom=pars["geom"], - geom_true=pars["geom_true"], + energy_axis_true=pars["e_true"], offset_max="2 deg", background_oversampling=pars.get("background_oversampling"), ) @@ -157,19 +169,24 @@ def test_map_maker_obs(observations): # Test for different spatial geoms and etrue, ereco bins geom_reco = geom(ebounds=[0.1, 1, 10]) - geom_true = geom(ebounds=[0.1, 0.5, 2.5, 10.0], binsz=1.0) + e_true = MapAxis.from_edges([0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log") geom_exp = geom(ebounds=[0.1, 0.5, 2.5, 10.0]) maker_obs = MapDatasetMaker( - geom=geom_reco, geom_true=geom_true, offset_max=2.0 * u.deg, cutout=False + geom=geom_reco, + energy_axis_true=e_true, + binsz_irf=1.0, + margin_irf=1.0, + offset_max=2.0 * u.deg, + cutout=False ) map_dataset = maker_obs.run(observations[0]) assert map_dataset.counts.geom == geom_reco assert map_dataset.background_model.map.geom == geom_reco assert map_dataset.exposure.geom == geom_exp - assert map_dataset.edisp.edisp_map.data.shape == (3, 48, 5, 10) - assert map_dataset.edisp.exposure_map.data.shape == (3, 1, 5, 10) - assert map_dataset.psf.psf_map.data.shape == (3, 66, 5, 10) - assert map_dataset.psf.exposure_map.data.shape == (3, 1, 5, 10) + assert map_dataset.edisp.edisp_map.data.shape == (3, 48, 6, 11) + assert map_dataset.edisp.exposure_map.data.shape == (3, 1, 6, 11) + assert map_dataset.psf.psf_map.data.shape == (3, 66, 6, 11) + assert map_dataset.psf.exposure_map.data.shape == (3, 1, 6, 11) assert_allclose(map_dataset.gti.time_delta, 1800.0 * u.s) assert map_dataset.name == "obs_110380" From 618e759e9dc14cefd585d9445b9bc6ec5e5582e2 Mon Sep 17 00:00:00 2001 From: AtreyeeS Date: Sat, 19 Oct 2019 01:55:22 +0200 Subject: [PATCH 03/11] change BINSZ to BINSZ_DEFAULT --- gammapy/cube/make.py | 2 +- gammapy/cube/tests/test_make.py | 36 +++++++++++++++++++++++---------- gammapy/scripts/analysis.py | 4 ++-- 3 files changed, 28 insertions(+), 14 deletions(-) diff --git a/gammapy/cube/make.py b/gammapy/cube/make.py index 77b8677a46..e699b6f2c1 100644 --- a/gammapy/cube/make.py +++ b/gammapy/cube/make.py @@ -18,7 +18,7 @@ MIGRA_AXIS_DEFAULT, RAD_AXIS_DEFAULT, ENERGY_AXIS_DEFAULT, - MapDataset + MapDataset, ) from .psf_map import make_psf_map diff --git a/gammapy/cube/tests/test_make.py b/gammapy/cube/tests/test_make.py index f43fd5f1b6..a69887b195 100644 --- a/gammapy/cube/tests/test_make.py +++ b/gammapy/cube/tests/test_make.py @@ -38,7 +38,9 @@ def geom(ebounds, binsz=0.5): { # Default, same e_true and reco "geom": geom(ebounds=[0.1, 1, 10]), - "e_true": MapAxis.from_edges([0.1, 1, 10], name="energy", unit="TeV", interp="log"), + "e_true": MapAxis.from_edges( + [0.1, 1, 10], name="energy", unit="TeV", interp="log" + ), "counts": 34366, "exposure": 9.995376e08, "exposure_image": 7.921993e10, @@ -49,7 +51,9 @@ def geom(ebounds, binsz=0.5): { # Test single energy bin "geom": geom(ebounds=[0.1, 10]), - "e_true": MapAxis.from_edges([0.1, 10], name="energy", unit="TeV", interp="log"), + "e_true": MapAxis.from_edges( + [0.1, 10], name="energy", unit="TeV", interp="log" + ), "counts": 34366, "exposure": 5.843302e08, "exposure_image": 1.16866e11, @@ -60,7 +64,9 @@ def geom(ebounds, binsz=0.5): { # Test single energy bin with exclusion mask "geom": geom(ebounds=[0.1, 10]), - "e_true": MapAxis.from_edges([0.1, 10], name="energy", unit="TeV", interp="log"), + "e_true": MapAxis.from_edges( + [0.1, 10], name="energy", unit="TeV", interp="log" + ), "exclusion_mask": Map.from_geom(geom(ebounds=[0.1, 10])), "counts": 34366, "exposure": 5.843302e08, @@ -72,7 +78,9 @@ def geom(ebounds, binsz=0.5): { # Test for different e_true and e_reco bins "geom": geom(ebounds=[0.1, 1, 10]), - "e_true": MapAxis.from_edges([0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log"), + "e_true": MapAxis.from_edges( + [0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log" + ), "counts": 34366, "exposure": 9.951827e08, "exposure_image": 6.492968e10, @@ -84,7 +92,9 @@ def geom(ebounds, binsz=0.5): { # Test for different e_true and e_reco bins "geom": geom(ebounds=[0.1, 1, 10]), - "e_true": MapAxis.from_edges([0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log"), + "e_true": MapAxis.from_edges( + [0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log" + ), "counts": 34366, "exposure": 9.951827e08, "exposure_image": 6.492968e10, @@ -97,10 +107,12 @@ def geom(ebounds, binsz=0.5): ) def test_map_maker(pars, observations): - stacked = MapDataset.create(geom=pars["geom"], - energy_axis_true=pars["e_true"], - binsz_irf=pars["binsz_irf"], - margin_irf=pars["margin_irf"]) + stacked = MapDataset.create( + geom=pars["geom"], + energy_axis_true=pars["e_true"], + binsz_irf=pars["binsz_irf"], + margin_irf=pars["margin_irf"], + ) for obs in observations: maker = MapDatasetMaker( @@ -169,7 +181,9 @@ def test_map_maker_obs(observations): # Test for different spatial geoms and etrue, ereco bins geom_reco = geom(ebounds=[0.1, 1, 10]) - e_true = MapAxis.from_edges([0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log") + e_true = MapAxis.from_edges( + [0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log" + ) geom_exp = geom(ebounds=[0.1, 0.5, 2.5, 10.0]) maker_obs = MapDatasetMaker( geom=geom_reco, @@ -177,7 +191,7 @@ def test_map_maker_obs(observations): binsz_irf=1.0, margin_irf=1.0, offset_max=2.0 * u.deg, - cutout=False + cutout=False, ) map_dataset = maker_obs.run(observations[0]) diff --git a/gammapy/scripts/analysis.py b/gammapy/scripts/analysis.py index 741adbb649..981c475354 100644 --- a/gammapy/scripts/analysis.py +++ b/gammapy/scripts/analysis.py @@ -10,7 +10,7 @@ import jsonschema import yaml from gammapy.cube import MapDataset, MapDatasetMaker -from gammapy.cube.fit import BINSZ_IRF +from gammapy.cube.fit import BINSZ_IRF_DEFAULT from gammapy.data import DataStore, ObservationTable from gammapy.maps import Map, MapAxis, WcsGeom from gammapy.modeling import Datasets, Fit @@ -245,7 +245,7 @@ def _map_making(self): if "geom-irf" in self.settings["datasets"]: geom_irf = self._create_geometry(self.settings["datasets"]["geom-irf"]) else: - geom_irf = geom.to_binsz(binsz=BINSZ_IRF) + geom_irf = geom.to_binsz(binsz=BINSZ_IRF_DEFAULT) offset_max = Angle(self.settings["datasets"]["offset-max"]) stack_datasets = self.settings["datasets"]["stack-datasets"] From acb4cc07ac9bd149ee50d733b4289d99d388d09a Mon Sep 17 00:00:00 2001 From: AtreyeeS Date: Fri, 18 Oct 2019 20:01:39 +0200 Subject: [PATCH 04/11] Adapt MapDataset --- gammapy/cube/fit.py | 81 ++++++++++++++++++---------------- gammapy/cube/tests/test_fit.py | 28 +++++++----- 2 files changed, 61 insertions(+), 48 deletions(-) diff --git a/gammapy/cube/fit.py b/gammapy/cube/fit.py index 77f29b3d53..c135fabc98 100644 --- a/gammapy/cube/fit.py +++ b/gammapy/cube/fit.py @@ -12,7 +12,7 @@ from gammapy.cube.psf_map import PSFMap from gammapy.data import GTI from gammapy.irf import EffectiveAreaTable, EnergyDispersion, apply_containment_fraction -from gammapy.maps import Map, MapAxis +from gammapy.maps import Map, MapAxis, WcsGeom from gammapy.modeling import Dataset, Parameters from gammapy.modeling.models import BackgroundModel, SkyModel, SkyModels from gammapy.spectrum import SpectrumDataset @@ -33,7 +33,12 @@ MIGRA_AXIS_DEFAULT = MapAxis.from_bounds( 0.2, 5, nbin=48, node_type="edges", name="migra" ) +ENERGY_AXIS_DEFAULT = MapAxis.from_edges( + np.logspace(-2.0, 2.5, 109), name="energy", unit=u.TeV, interp="log" +) BINSZ_IRF = 0.2 +BINSZ_IRF_DEFAULT = BINSZ_IRF +MARGIN_IRF_DEFAULT = 0.5 # TODO: Choose optimal binnings depending on IRFs @@ -363,9 +368,11 @@ def from_geoms( def create( cls, geom, - geom_irf=None, + energy_axis_true=None, migra_axis=None, rad_axis=None, + binsz_irf=None, + margin_irf=None, reference_time="2000-01-01", name="", **kwargs, @@ -376,54 +383,54 @@ def create( ---------- geom: `~gammapy.maps.WcsGeom` Reference target geometry in reco energy, used for counts and background maps - geom_irf: `~gammapy.maps.WcsGeom` - Reference image geometry in true energy, used for IRF maps. + energy_axis_true: `~gammapy.maps.MapAxis` + True energy axis used for IRF maps migra_axis: `~gammapy.maps.MapAxis` Migration axis for the energy dispersion map rad_axis: `~gammapy.maps.MapAxis` Rad axis for the psf map + binsz_irf_maps: float + IRF Map pixel size in degrees. + margin_irf: float + IRF map margin size in degrees reference_time: `~astropy.time.Time` the reference time to use in GTI definition name : str Name of the dataset. + + Returns + -------- + empty_maps: `MapDataset` + A MapDataset containing zero filled maps """ - geom_irf = geom_irf or geom.to_binsz(BINSZ_IRF) + migra_axis = migra_axis or MIGRA_AXIS_DEFAULT rad_axis = rad_axis or RAD_AXIS_DEFAULT + energy_axis_true = energy_axis_true or ENERGY_AXIS_DEFAULT + binsz_irf = binsz_irf or BINSZ_IRF_DEFAULT + margin_irf = margin_irf or MARGIN_IRF_DEFAULT + margin_irf = margin_irf * u.deg + + wcs = geom.to_image() + geom_exposure = wcs.to_cube([energy_axis_true]) + + wcs_irf = WcsGeom.create( + binsz=binsz_irf, + width=wcs.width + margin_irf, + skydir=wcs.center_skydir, + proj=wcs.projection, + coordsys=wcs.coordsys, + ) - counts = Map.from_geom(geom, unit="") - - background = Map.from_geom(geom, unit="") - background_model = BackgroundModel(background) - - energy_axis = geom_irf.get_axis_by_name("ENERGY") - - exposure_geom = geom.to_image().to_cube([energy_axis]) - exposure = Map.from_geom(exposure_geom, unit="m2 s") - exposure_irf = Map.from_geom(geom_irf, unit="m2 s") - - mask_safe = np.zeros(geom.data_shape, dtype=bool) - - gti = GTI.create([] * u.s, [] * u.s, reference_time=reference_time) - - geom_migra = geom_irf.to_image().to_cube([migra_axis, energy_axis]) - edisp_map = Map.from_geom(geom_migra, unit="") - loc = migra_axis.edges.searchsorted(1.0) - edisp_map.data[:, loc, :, :] = 1.0 - edisp = EDispMap(edisp_map, exposure_irf) - - geom_rad = geom_irf.to_image().to_cube([rad_axis, energy_axis]) - psf_map = Map.from_geom(geom_rad, unit="sr-1") - psf = PSFMap(psf_map, exposure_irf) + geom_psf = wcs_irf.to_cube([rad_axis, energy_axis_true]) + geom_edisp = wcs_irf.to_cube([migra_axis, energy_axis_true]) - return cls( - counts=counts, - exposure=exposure, - psf=psf, - edisp=edisp, - background_model=background_model, - gti=gti, - mask_safe=mask_safe, + return cls.from_geoms( + geom, + geom_exposure, + geom_psf, + geom_edisp, + reference_time=reference_time, name=name, **kwargs, ) diff --git a/gammapy/cube/tests/test_fit.py b/gammapy/cube/tests/test_fit.py index f827801cd2..2cb6635e93 100644 --- a/gammapy/cube/tests/test_fit.py +++ b/gammapy/cube/tests/test_fit.py @@ -360,22 +360,29 @@ def test_map_fit_one_energy_bin(sky_model, geom_image): def test_create(geom, geom_etrue): # tests empty datasets created - migra_axis = MapAxis(nodes=np.linspace(0.0, 3.0, 51), unit="", name="migra") rad_axis = MapAxis(nodes=np.linspace(0.0, 1.0, 51), unit="deg", name="theta") + e_reco = MapAxis.from_edges( + np.logspace(-1.0, 1.0, 3), name="energy", unit=u.TeV, interp="log" + ) + e_true = MapAxis.from_edges( + np.logspace(-1.0, 1.0, 4), name="energy", unit=u.TeV, interp="log" + ) + geom = WcsGeom.create(binsz=0.02, width=(2, 2), axes=[e_reco]) + empty_dataset = MapDataset.create( + geom=geom, energy_axis_true=e_true, migra_axis=migra_axis, rad_axis=rad_axis + ) - empty_dataset = MapDataset.create(geom, geom_etrue, migra_axis, rad_axis) - - assert_allclose(empty_dataset.counts.data.sum(), 0.0) - assert_allclose(empty_dataset.background_model.map.data.sum(), 0.0) + assert empty_dataset.counts.data.shape == (2, 100, 100) - assert empty_dataset.psf.psf_map.data.shape == (3, 50, 100, 100) - assert empty_dataset.psf.exposure_map.data.shape == (3, 1, 100, 100) + assert empty_dataset.exposure.data.shape == (3, 100, 100) - assert empty_dataset.edisp.edisp_map.data.shape == (3, 50, 100, 100) - assert empty_dataset.edisp.exposure_map.data.shape == (3, 1, 100, 100) + assert empty_dataset.psf.psf_map.data.shape == (3, 50, 12, 12) + assert empty_dataset.psf.exposure_map.data.shape == (3, 1, 12, 12) - assert_allclose(empty_dataset.edisp.edisp_map.data.sum(), 30000) + assert empty_dataset.edisp.edisp_map.data.shape == (3, 50, 12, 12) + assert empty_dataset.edisp.exposure_map.data.shape == (3, 1, 12, 12) + assert_allclose(empty_dataset.edisp.edisp_map.data.sum(), 432) assert_allclose(empty_dataset.gti.time_delta, 0.0 * u.s) @@ -394,7 +401,6 @@ def test_from_geoms(): wcs_irf = WcsGeom.create(binsz=0.1, width=(2.5, 2.5)) geom = wcs.to_cube([e_reco]) geom_exposure = wcs.to_cube([e_true]) - geom_irf = wcs_irf.to_cube([e_true]) geom_psf = wcs_irf.to_cube([rad_axis, e_true]) geom_edisp = wcs_irf.to_cube([migra_axis, e_true]) From 51e7297fe4765432e9f11bb9901022294460d1af Mon Sep 17 00:00:00 2001 From: AtreyeeS Date: Sat, 19 Oct 2019 01:45:57 +0200 Subject: [PATCH 05/11] Adapt MapDatasetMaker --- gammapy/cube/fit.py | 4 +- gammapy/cube/make.py | 65 +++++++++++++++++++++++++-------- gammapy/cube/tests/test_make.py | 49 +++++++++++++++++-------- 3 files changed, 84 insertions(+), 34 deletions(-) diff --git a/gammapy/cube/fit.py b/gammapy/cube/fit.py index c135fabc98..e247e208bd 100644 --- a/gammapy/cube/fit.py +++ b/gammapy/cube/fit.py @@ -36,8 +36,8 @@ ENERGY_AXIS_DEFAULT = MapAxis.from_edges( np.logspace(-2.0, 2.5, 109), name="energy", unit=u.TeV, interp="log" ) -BINSZ_IRF = 0.2 -BINSZ_IRF_DEFAULT = BINSZ_IRF + +BINSZ_IRF_DEFAULT = 0.2 MARGIN_IRF_DEFAULT = 0.5 # TODO: Choose optimal binnings depending on IRFs diff --git a/gammapy/cube/make.py b/gammapy/cube/make.py index 49cb491735..77b8677a46 100644 --- a/gammapy/cube/make.py +++ b/gammapy/cube/make.py @@ -2,16 +2,24 @@ import logging from functools import lru_cache import numpy as np +import astropy.units as u from astropy.coordinates import Angle from astropy.nddata.utils import NoOverlapError from astropy.utils import lazyproperty from gammapy.irf import EnergyDependentMultiGaussPSF -from gammapy.maps import Map +from gammapy.maps import Map, WcsGeom from gammapy.modeling.models import BackgroundModel from .background import make_map_background_irf from .edisp_map import make_edisp_map from .exposure import _map_spectrum_weight, make_map_exposure_true_energy -from .fit import BINSZ_IRF, MIGRA_AXIS_DEFAULT, RAD_AXIS_DEFAULT, MapDataset +from .fit import ( + BINSZ_IRF_DEFAULT, + MARGIN_IRF_DEFAULT, + MIGRA_AXIS_DEFAULT, + RAD_AXIS_DEFAULT, + ENERGY_AXIS_DEFAULT, + MapDataset +) from .psf_map import make_psf_map __all__ = ["MapDatasetMaker", "MapMakerRing"] @@ -28,14 +36,16 @@ class MapDatasetMaker: Reference image geometry in reco energy, used for counts and background maps offset_max : `~astropy.coordinates.Angle` Maximum offset angle - geom_true : `~gammapy.maps.WcsGeom` - Reference image geometry in true energy, used for IRF maps. It can have a coarser - spatial bins than the counts geom. - If none, the same as geom is assumed + energy_axis_true: `~gammapy.maps.MapAxis` + True energy axis used for IRF maps migra_axis : `~gammapy.maps.MapAxis` Migration axis for edisp map rad_axis : `~gammapy.maps.MapAxis` Radial axis for psf map. + binsz_irf: float + IRF Map pixel size in degrees. + margin_irf: float + IRF map margin size in degrees cutout : bool Whether to cutout the observation. cutout_mode : {'trim', 'partial', 'strict'} @@ -47,19 +57,27 @@ def __init__( self, geom, offset_max, - geom_true=None, background_oversampling=None, + energy_axis_true=None, migra_axis=None, rad_axis=None, + binsz_irf=None, + margin_irf=None, cutout_mode="trim", cutout=True, ): + self.geom = geom - self.geom_true = geom_true if geom_true else geom.to_binsz(BINSZ_IRF) self.offset_max = Angle(offset_max) self.background_oversampling = background_oversampling self.migra_axis = migra_axis if migra_axis else MIGRA_AXIS_DEFAULT self.rad_axis = rad_axis if rad_axis else RAD_AXIS_DEFAULT + self.energy_axis_true = energy_axis_true or ENERGY_AXIS_DEFAULT + self.binsz_irf = binsz_irf or BINSZ_IRF_DEFAULT + + self.margin_irf = margin_irf or MARGIN_IRF_DEFAULT + self.margin_irf = self.margin_irf * u.deg + self.cutout_mode = cutout_mode self.cutout_width = 2 * self.offset_max self.cutout = cutout @@ -75,25 +93,39 @@ def _cutout_geom(self, geom, observation): else: return geom + @lazyproperty + def wcs_irf(self): + """Spatial geometry of IRF Maps (`Geom`)""" + wcs = self.geom.to_image() + return WcsGeom.create( + binsz=self.binsz_irf, + width=wcs.width + self.margin_irf, + skydir=wcs.center_skydir, + proj=wcs.projection, + coordsys=wcs.coordsys, + ) + + @lazyproperty + def geom_irf(self): + """Geom of Exposure map associated with IRFs (`Geom`)""" + return self.wcs_irf.to_cube([self.energy_axis_true]) + @lazyproperty def geom_exposure(self): """Exposure map geom (`Geom`)""" - energy_axis = self.geom_true.get_axis_by_name("energy") - geom_exposure = self.geom.to_image().to_cube([energy_axis]) + geom_exposure = self.geom.to_image().to_cube([self.energy_axis_true]) return geom_exposure @lazyproperty def geom_psf(self): """PSFMap geom (`Geom`)""" - energy_axis = self.geom_true.get_axis_by_name("ENERGY") - geom_psf = self.geom_true.to_image().to_cube([self.rad_axis, energy_axis]) + geom_psf = self.wcs_irf.to_cube([self.rad_axis, self.energy_axis_true]) return geom_psf @lazyproperty def geom_edisp(self): """EdispMap geom (`Geom`)""" - energy_axis = self.geom_true.get_axis_by_name("ENERGY") - geom_edisp = self.geom_true.to_image().to_cube([self.migra_axis, energy_axis]) + geom_edisp = self.wcs_irf.to_cube([self.migra_axis, self.energy_axis_true]) return geom_edisp def make_counts(self, observation): @@ -150,7 +182,8 @@ def make_exposure_irf(self, observation): exposure : `Map` Exposure map. """ - geom = self._cutout_geom(self.geom_true, observation) + + geom = self._cutout_geom(self.geom_irf, observation) exposure = make_map_exposure_true_energy( pointing=observation.pointing_radec, livetime=observation.observation_live_time_duration, @@ -281,7 +314,7 @@ def make_mask_safe_irf(self, observation): mask : `~numpy.ndarray` Mask """ - geom = self._cutout_geom(self.geom_true, observation) + geom = self._cutout_geom(self.geom_irf, observation) offset = geom.separation(observation.pointing_radec) return offset >= self.offset_max diff --git a/gammapy/cube/tests/test_make.py b/gammapy/cube/tests/test_make.py index 5e3d48d141..f43fd5f1b6 100644 --- a/gammapy/cube/tests/test_make.py +++ b/gammapy/cube/tests/test_make.py @@ -36,64 +36,76 @@ def geom(ebounds, binsz=0.5): "pars", [ { - # Default, normal test case + # Default, same e_true and reco "geom": geom(ebounds=[0.1, 1, 10]), - "geom_true": None, + "e_true": MapAxis.from_edges([0.1, 1, 10], name="energy", unit="TeV", interp="log"), "counts": 34366, "exposure": 9.995376e08, "exposure_image": 7.921993e10, "background": 27989.05, + "binsz_irf": 0.5, + "margin_irf": 0.0, }, { # Test single energy bin "geom": geom(ebounds=[0.1, 10]), - "geom_true": None, + "e_true": MapAxis.from_edges([0.1, 10], name="energy", unit="TeV", interp="log"), "counts": 34366, "exposure": 5.843302e08, "exposure_image": 1.16866e11, "background": 30424.451, + "binsz_irf": 0.5, + "margin_irf": 0.0, }, { # Test single energy bin with exclusion mask "geom": geom(ebounds=[0.1, 10]), - "geom_true": None, + "e_true": MapAxis.from_edges([0.1, 10], name="energy", unit="TeV", interp="log"), "exclusion_mask": Map.from_geom(geom(ebounds=[0.1, 10])), "counts": 34366, "exposure": 5.843302e08, "exposure_image": 1.16866e11, "background": 30424.451, + "binsz_irf": 0.5, + "margin_irf": 0.0, }, { # Test for different e_true and e_reco bins "geom": geom(ebounds=[0.1, 1, 10]), - "geom_true": geom(ebounds=[0.1, 0.5, 2.5, 10.0]), + "e_true": MapAxis.from_edges([0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log"), "counts": 34366, "exposure": 9.951827e08, "exposure_image": 6.492968e10, "background": 28760.283, "background_oversampling": 2, + "binsz_irf": 0.5, + "margin_irf": 0.0, }, { # Test for different e_true and e_reco bins "geom": geom(ebounds=[0.1, 1, 10]), - "geom_true": geom(ebounds=[0.1, 0.5, 2.5, 10.0], binsz=1), + "e_true": MapAxis.from_edges([0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log"), "counts": 34366, "exposure": 9.951827e08, "exposure_image": 6.492968e10, "background": 28760.283, "background_oversampling": 2, + "binsz_irf": 1.0, + "margin_irf": 0.0, }, ], ) -@pytest.mark.parametrize("keepdims", [True, False]) -def test_map_maker(pars, observations, keepdims): +def test_map_maker(pars, observations): - stacked = MapDataset.create(geom=pars["geom"], geom_irf=pars["geom_true"]) + stacked = MapDataset.create(geom=pars["geom"], + energy_axis_true=pars["e_true"], + binsz_irf=pars["binsz_irf"], + margin_irf=pars["margin_irf"]) for obs in observations: maker = MapDatasetMaker( geom=pars["geom"], - geom_true=pars["geom_true"], + energy_axis_true=pars["e_true"], offset_max="2 deg", background_oversampling=pars.get("background_oversampling"), ) @@ -157,19 +169,24 @@ def test_map_maker_obs(observations): # Test for different spatial geoms and etrue, ereco bins geom_reco = geom(ebounds=[0.1, 1, 10]) - geom_true = geom(ebounds=[0.1, 0.5, 2.5, 10.0], binsz=1.0) + e_true = MapAxis.from_edges([0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log") geom_exp = geom(ebounds=[0.1, 0.5, 2.5, 10.0]) maker_obs = MapDatasetMaker( - geom=geom_reco, geom_true=geom_true, offset_max=2.0 * u.deg, cutout=False + geom=geom_reco, + energy_axis_true=e_true, + binsz_irf=1.0, + margin_irf=1.0, + offset_max=2.0 * u.deg, + cutout=False ) map_dataset = maker_obs.run(observations[0]) assert map_dataset.counts.geom == geom_reco assert map_dataset.background_model.map.geom == geom_reco assert map_dataset.exposure.geom == geom_exp - assert map_dataset.edisp.edisp_map.data.shape == (3, 48, 5, 10) - assert map_dataset.edisp.exposure_map.data.shape == (3, 1, 5, 10) - assert map_dataset.psf.psf_map.data.shape == (3, 66, 5, 10) - assert map_dataset.psf.exposure_map.data.shape == (3, 1, 5, 10) + assert map_dataset.edisp.edisp_map.data.shape == (3, 48, 6, 11) + assert map_dataset.edisp.exposure_map.data.shape == (3, 1, 6, 11) + assert map_dataset.psf.psf_map.data.shape == (3, 66, 6, 11) + assert map_dataset.psf.exposure_map.data.shape == (3, 1, 6, 11) assert_allclose(map_dataset.gti.time_delta, 1800.0 * u.s) assert map_dataset.name == "obs_110380" From 7df9907fae5dc92d3273fbfaa7d2cd5e1b48c7cf Mon Sep 17 00:00:00 2001 From: AtreyeeS Date: Sat, 19 Oct 2019 01:55:22 +0200 Subject: [PATCH 06/11] change BINSZ to BINSZ_DEFAULT --- gammapy/cube/make.py | 2 +- gammapy/cube/tests/test_make.py | 36 +++++++++++++++++++++++---------- gammapy/scripts/analysis.py | 4 ++-- 3 files changed, 28 insertions(+), 14 deletions(-) diff --git a/gammapy/cube/make.py b/gammapy/cube/make.py index 77b8677a46..e699b6f2c1 100644 --- a/gammapy/cube/make.py +++ b/gammapy/cube/make.py @@ -18,7 +18,7 @@ MIGRA_AXIS_DEFAULT, RAD_AXIS_DEFAULT, ENERGY_AXIS_DEFAULT, - MapDataset + MapDataset, ) from .psf_map import make_psf_map diff --git a/gammapy/cube/tests/test_make.py b/gammapy/cube/tests/test_make.py index f43fd5f1b6..a69887b195 100644 --- a/gammapy/cube/tests/test_make.py +++ b/gammapy/cube/tests/test_make.py @@ -38,7 +38,9 @@ def geom(ebounds, binsz=0.5): { # Default, same e_true and reco "geom": geom(ebounds=[0.1, 1, 10]), - "e_true": MapAxis.from_edges([0.1, 1, 10], name="energy", unit="TeV", interp="log"), + "e_true": MapAxis.from_edges( + [0.1, 1, 10], name="energy", unit="TeV", interp="log" + ), "counts": 34366, "exposure": 9.995376e08, "exposure_image": 7.921993e10, @@ -49,7 +51,9 @@ def geom(ebounds, binsz=0.5): { # Test single energy bin "geom": geom(ebounds=[0.1, 10]), - "e_true": MapAxis.from_edges([0.1, 10], name="energy", unit="TeV", interp="log"), + "e_true": MapAxis.from_edges( + [0.1, 10], name="energy", unit="TeV", interp="log" + ), "counts": 34366, "exposure": 5.843302e08, "exposure_image": 1.16866e11, @@ -60,7 +64,9 @@ def geom(ebounds, binsz=0.5): { # Test single energy bin with exclusion mask "geom": geom(ebounds=[0.1, 10]), - "e_true": MapAxis.from_edges([0.1, 10], name="energy", unit="TeV", interp="log"), + "e_true": MapAxis.from_edges( + [0.1, 10], name="energy", unit="TeV", interp="log" + ), "exclusion_mask": Map.from_geom(geom(ebounds=[0.1, 10])), "counts": 34366, "exposure": 5.843302e08, @@ -72,7 +78,9 @@ def geom(ebounds, binsz=0.5): { # Test for different e_true and e_reco bins "geom": geom(ebounds=[0.1, 1, 10]), - "e_true": MapAxis.from_edges([0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log"), + "e_true": MapAxis.from_edges( + [0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log" + ), "counts": 34366, "exposure": 9.951827e08, "exposure_image": 6.492968e10, @@ -84,7 +92,9 @@ def geom(ebounds, binsz=0.5): { # Test for different e_true and e_reco bins "geom": geom(ebounds=[0.1, 1, 10]), - "e_true": MapAxis.from_edges([0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log"), + "e_true": MapAxis.from_edges( + [0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log" + ), "counts": 34366, "exposure": 9.951827e08, "exposure_image": 6.492968e10, @@ -97,10 +107,12 @@ def geom(ebounds, binsz=0.5): ) def test_map_maker(pars, observations): - stacked = MapDataset.create(geom=pars["geom"], - energy_axis_true=pars["e_true"], - binsz_irf=pars["binsz_irf"], - margin_irf=pars["margin_irf"]) + stacked = MapDataset.create( + geom=pars["geom"], + energy_axis_true=pars["e_true"], + binsz_irf=pars["binsz_irf"], + margin_irf=pars["margin_irf"], + ) for obs in observations: maker = MapDatasetMaker( @@ -169,7 +181,9 @@ def test_map_maker_obs(observations): # Test for different spatial geoms and etrue, ereco bins geom_reco = geom(ebounds=[0.1, 1, 10]) - e_true = MapAxis.from_edges([0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log") + e_true = MapAxis.from_edges( + [0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log" + ) geom_exp = geom(ebounds=[0.1, 0.5, 2.5, 10.0]) maker_obs = MapDatasetMaker( geom=geom_reco, @@ -177,7 +191,7 @@ def test_map_maker_obs(observations): binsz_irf=1.0, margin_irf=1.0, offset_max=2.0 * u.deg, - cutout=False + cutout=False, ) map_dataset = maker_obs.run(observations[0]) diff --git a/gammapy/scripts/analysis.py b/gammapy/scripts/analysis.py index 741adbb649..981c475354 100644 --- a/gammapy/scripts/analysis.py +++ b/gammapy/scripts/analysis.py @@ -10,7 +10,7 @@ import jsonschema import yaml from gammapy.cube import MapDataset, MapDatasetMaker -from gammapy.cube.fit import BINSZ_IRF +from gammapy.cube.fit import BINSZ_IRF_DEFAULT from gammapy.data import DataStore, ObservationTable from gammapy.maps import Map, MapAxis, WcsGeom from gammapy.modeling import Datasets, Fit @@ -245,7 +245,7 @@ def _map_making(self): if "geom-irf" in self.settings["datasets"]: geom_irf = self._create_geometry(self.settings["datasets"]["geom-irf"]) else: - geom_irf = geom.to_binsz(binsz=BINSZ_IRF) + geom_irf = geom.to_binsz(binsz=BINSZ_IRF_DEFAULT) offset_max = Angle(self.settings["datasets"]["offset-max"]) stack_datasets = self.settings["datasets"]["stack-datasets"] From f116b93ce885358d460afda965b256ddf5a9c184 Mon Sep 17 00:00:00 2001 From: Jose Enrique Ruiz Date: Sat, 19 Oct 2019 21:19:07 +0200 Subject: [PATCH 07/11] Fix high-level interface --- gammapy/scripts/analysis.py | 28 +++++++++++++++---------- gammapy/scripts/config/docs.yaml | 12 ++++++++++- gammapy/scripts/config/schema.yaml | 8 ++++++- gammapy/scripts/config/template-3d.yaml | 5 +---- tutorials/hess.ipynb | 11 ++++++++++ 5 files changed, 47 insertions(+), 17 deletions(-) diff --git a/gammapy/scripts/analysis.py b/gammapy/scripts/analysis.py index 981c475354..115871f254 100644 --- a/gammapy/scripts/analysis.py +++ b/gammapy/scripts/analysis.py @@ -10,7 +10,6 @@ import jsonschema import yaml from gammapy.cube import MapDataset, MapDatasetMaker -from gammapy.cube.fit import BINSZ_IRF_DEFAULT from gammapy.data import DataStore, ObservationTable from gammapy.maps import Map, MapAxis, WcsGeom from gammapy.modeling import Datasets, Fit @@ -225,35 +224,42 @@ def get_flux_points(self, source="source"): @staticmethod def _create_geometry(params): """Create the geometry.""" - # TODO: handled in jsonschema validation class geom_params = copy.deepcopy(params) - axes = [] - for axis_params in geom_params.get("axes", []): + for axis_params in params.get("axes", []): ax = MapAxis.from_bounds(**axis_params) axes.append(ax) - geom_params["axes"] = axes - geom_params["skydir"] = tuple(geom_params["skydir"]) + if "skydir" in geom_params: + geom_params["skydir"] = tuple(geom_params["skydir"]) return WcsGeom.create(**geom_params) def _map_making(self): """Make maps and datasets for 3d analysis.""" log.info("Creating geometry.") + geom = self._create_geometry(self.settings["datasets"]["geom"]) + geom_irf = dict( + energy_axis_true=None, + binsz_irf=None, + margin_irf=None + ) if "geom-irf" in self.settings["datasets"]: - geom_irf = self._create_geometry(self.settings["datasets"]["geom-irf"]) - else: - geom_irf = geom.to_binsz(binsz=BINSZ_IRF_DEFAULT) + settings_irf = self.settings["datasets"]["geom-irf"] + for axis_params in settings_irf.get("axes", []): + if axis_params["name"] == "energy": + geom_irf["energy_axis_true"] = MapAxis.from_bounds(**axis_params) + geom_irf["binsz_irf"] = settings_irf.get("binsz", None) + geom_irf["margin_irf"] = settings_irf.get("margin", None) offset_max = Angle(self.settings["datasets"]["offset-max"]) stack_datasets = self.settings["datasets"]["stack-datasets"] log.info("Creating datasets.") - maker = MapDatasetMaker(geom=geom, geom_true=geom_irf, offset_max=offset_max) + maker = MapDatasetMaker(geom=geom, offset_max=offset_max, **geom_irf) if stack_datasets: - stacked = MapDataset.create(geom=geom, geom_irf=geom_irf, name="stacked") + stacked = MapDataset.create(geom=geom, name="stacked", **geom_irf) for obs in self.observations: dataset = maker.run(obs) stacked.stack(dataset) diff --git a/gammapy/scripts/config/docs.yaml b/gammapy/scripts/config/docs.yaml index 9207ed4ce7..0f62261f64 100644 --- a/gammapy/scripts/config/docs.yaml +++ b/gammapy/scripts/config/docs.yaml @@ -54,7 +54,6 @@ datasets: hdu: IMAGE containment_correction: true # geometry settings - # note that there may be different values for geom and geom-irf geom: # on region in spectrum extraction region: @@ -76,6 +75,17 @@ datasets: nbin: 73 unit: TeV interp: log + geom-irf: + # geom values for the IRF + margin: 0.05 # IRF map pixel size in degrees + binsz: 0.02 # IRF map pixel size in degrees + axes: + - name: energy + lo_bnd: 0.01 + hi_bnd: 100 + nbin: 73 + unit: TeV + interp: log # Section: fit # Fitting process diff --git a/gammapy/scripts/config/schema.yaml b/gammapy/scripts/config/schema.yaml index 3c553f5364..0f52debb47 100644 --- a/gammapy/scripts/config/schema.yaml +++ b/gammapy/scripts/config/schema.yaml @@ -100,7 +100,13 @@ properties: geom: "$ref": "#/definitions/map_geom" geom-irf: - "$ref": "#/definitions/map_geom" + margin: {type: number} + binsz: {type: number} + axes: + type: array + additionalProperties: false + items: + "$ref": "#/definitions/map_axis" required: [dataset-type, geom, stack-datasets] allOf: - if: diff --git a/gammapy/scripts/config/template-3d.yaml b/gammapy/scripts/config/template-3d.yaml index fca3c5bb60..97f3a520bb 100644 --- a/gammapy/scripts/config/template-3d.yaml +++ b/gammapy/scripts/config/template-3d.yaml @@ -31,11 +31,8 @@ datasets: node_type: edges unit: TeV geom-irf: - skydir: [83.633, 22.014] - width: [5, 5] + margin: 0.5 binsz: 0.2 - coordsys: CEL - proj: TAN axes: - name: energy hi_bnd: 10 diff --git a/tutorials/hess.ipynb b/tutorials/hess.ipynb index 5af63b09c8..9311344530 100644 --- a/tutorials/hess.ipynb +++ b/tutorials/hess.ipynb @@ -102,6 +102,17 @@ " interp: log\n", " node_type: edges\n", " unit: TeV\n", + " geom-irf:\n", + " margin: 0.5\n", + " binsz: 0.2\n", + " axes:\n", + " - name: energy\n", + " hi_bnd: 10\n", + " lo_bnd: 1\n", + " nbin: 11\n", + " interp: log\n", + " node_type: edges\n", + " unit: TeV \n", "\n", "fit:\n", " fit_range:\n", From d936eb93f4a7817b4ea1cee387a7095386d86c8c Mon Sep 17 00:00:00 2001 From: Jose Enrique Ruiz Date: Tue, 22 Oct 2019 15:39:16 +0200 Subject: [PATCH 08/11] address review comments --- gammapy/scripts/analysis.py | 12 +++++------- gammapy/scripts/config/docs.yaml | 20 +++++++++----------- gammapy/scripts/config/schema.yaml | 13 +++++-------- gammapy/scripts/config/template-3d.yaml | 22 +++++++++++----------- gammapy/scripts/tests/test_analysis.py | 1 - tutorials/hess.ipynb | 25 ++++++++++++------------- 6 files changed, 42 insertions(+), 51 deletions(-) diff --git a/gammapy/scripts/analysis.py b/gammapy/scripts/analysis.py index 115871f254..f593ff87b6 100644 --- a/gammapy/scripts/analysis.py +++ b/gammapy/scripts/analysis.py @@ -245,13 +245,11 @@ def _map_making(self): binsz_irf=None, margin_irf=None ) - if "geom-irf" in self.settings["datasets"]: - settings_irf = self.settings["datasets"]["geom-irf"] - for axis_params in settings_irf.get("axes", []): - if axis_params["name"] == "energy": - geom_irf["energy_axis_true"] = MapAxis.from_bounds(**axis_params) - geom_irf["binsz_irf"] = settings_irf.get("binsz", None) - geom_irf["margin_irf"] = settings_irf.get("margin", None) + if "energy-axis-true" in self.settings["datasets"]: + axis_params = self.settings["datasets"]["energy-axis-true"] + geom_irf["energy_axis_true"] = MapAxis.from_bounds(**axis_params) + geom_irf["binsz_irf"] = self.settings["datasets"].get("binsz", None) + geom_irf["margin_irf"] = self.settings["datasets"].get("margin", None) offset_max = Angle(self.settings["datasets"]["offset-max"]) stack_datasets = self.settings["datasets"]["stack-datasets"] diff --git a/gammapy/scripts/config/docs.yaml b/gammapy/scripts/config/docs.yaml index 0f62261f64..0b9bd07f72 100644 --- a/gammapy/scripts/config/docs.yaml +++ b/gammapy/scripts/config/docs.yaml @@ -75,17 +75,15 @@ datasets: nbin: 73 unit: TeV interp: log - geom-irf: - # geom values for the IRF - margin: 0.05 # IRF map pixel size in degrees - binsz: 0.02 # IRF map pixel size in degrees - axes: - - name: energy - lo_bnd: 0.01 - hi_bnd: 100 - nbin: 73 - unit: TeV - interp: log + energy-axis-true: + name: energy + lo_bnd: 0.01 + hi_bnd: 100 + nbin: 73 + unit: TeV + interp: log + binsz-irf: 0.02 # IRF map pixel size in degrees + margin-irf: 0.05 # IRF map pixel size in degrees # Section: fit # Fitting process diff --git a/gammapy/scripts/config/schema.yaml b/gammapy/scripts/config/schema.yaml index 0f52debb47..2e708d03ca 100644 --- a/gammapy/scripts/config/schema.yaml +++ b/gammapy/scripts/config/schema.yaml @@ -99,14 +99,11 @@ properties: psf-kernel-radius: {type: number} geom: "$ref": "#/definitions/map_geom" - geom-irf: - margin: {type: number} - binsz: {type: number} - axes: - type: array - additionalProperties: false - items: - "$ref": "#/definitions/map_axis" + energy-axis-true: + "$ref": "#/definitions/map_axis" + binsz-irf: {type: number} + margin-irf: {type: number} + required: [dataset-type, geom, stack-datasets] allOf: - if: diff --git a/gammapy/scripts/config/template-3d.yaml b/gammapy/scripts/config/template-3d.yaml index 97f3a520bb..e49e038c61 100644 --- a/gammapy/scripts/config/template-3d.yaml +++ b/gammapy/scripts/config/template-3d.yaml @@ -30,17 +30,17 @@ datasets: interp: log node_type: edges unit: TeV - geom-irf: - margin: 0.5 - binsz: 0.2 - axes: - - name: energy - hi_bnd: 10 - lo_bnd: 1 - nbin: 11 - interp: log - node_type: center - unit: TeV + energy-axis-true: + name: energy + hi_bnd: 10 + lo_bnd: 1 + nbin: 5 + interp: log + node_type: edges + unit: TeV + binsz-irf: 0.2 + margin-irf: 0.5 + fit: {} diff --git a/gammapy/scripts/tests/test_analysis.py b/gammapy/scripts/tests/test_analysis.py index 2037ecf5f1..e24137aded 100644 --- a/gammapy/scripts/tests/test_analysis.py +++ b/gammapy/scripts/tests/test_analysis.py @@ -240,7 +240,6 @@ def test_help(): def test_analysis_3d_no_geom_irf(): config = AnalysisConfig.from_template("3d") analysis = Analysis(config) - del analysis.settings["datasets"]["geom-irf"] analysis.get_observations() analysis.get_datasets() diff --git a/tutorials/hess.ipynb b/tutorials/hess.ipynb index 9311344530..2271e740ff 100644 --- a/tutorials/hess.ipynb +++ b/tutorials/hess.ipynb @@ -101,19 +101,18 @@ " nbin: 5\n", " interp: log\n", " node_type: edges\n", - " unit: TeV\n", - " geom-irf:\n", - " margin: 0.5\n", - " binsz: 0.2\n", - " axes:\n", - " - name: energy\n", - " hi_bnd: 10\n", - " lo_bnd: 1\n", - " nbin: 11\n", - " interp: log\n", - " node_type: edges\n", - " unit: TeV \n", - "\n", + " unit: TeV \n", + " energy-axis-true:\n", + " name: energy\n", + " hi_bnd: 10\n", + " lo_bnd: 1\n", + " nbin: 11\n", + " interp: log\n", + " node_type: edges\n", + " unit: TeV\n", + " binsz-irf: 0.2\n", + " margin-irf: 0.5\n", + " \n", "fit:\n", " fit_range:\n", " max: 30 TeV\n", From 4224305cb925dc846df4e0ba185570d7cb2f43be Mon Sep 17 00:00:00 2001 From: AtreyeeS Date: Wed, 23 Oct 2019 17:53:56 +0200 Subject: [PATCH 09/11] Make e_true default to e_ereco --- gammapy/cube/fit.py | 5 +---- gammapy/cube/make.py | 19 ++++++++++--------- gammapy/cube/tests/test_make.py | 14 ++++---------- 3 files changed, 15 insertions(+), 23 deletions(-) diff --git a/gammapy/cube/fit.py b/gammapy/cube/fit.py index e247e208bd..ce55ec1211 100644 --- a/gammapy/cube/fit.py +++ b/gammapy/cube/fit.py @@ -33,9 +33,6 @@ MIGRA_AXIS_DEFAULT = MapAxis.from_bounds( 0.2, 5, nbin=48, node_type="edges", name="migra" ) -ENERGY_AXIS_DEFAULT = MapAxis.from_edges( - np.logspace(-2.0, 2.5, 109), name="energy", unit=u.TeV, interp="log" -) BINSZ_IRF_DEFAULT = 0.2 MARGIN_IRF_DEFAULT = 0.5 @@ -406,7 +403,7 @@ def create( migra_axis = migra_axis or MIGRA_AXIS_DEFAULT rad_axis = rad_axis or RAD_AXIS_DEFAULT - energy_axis_true = energy_axis_true or ENERGY_AXIS_DEFAULT + energy_axis_true = energy_axis_true or geom.get_axis_by_name("energy") binsz_irf = binsz_irf or BINSZ_IRF_DEFAULT margin_irf = margin_irf or MARGIN_IRF_DEFAULT margin_irf = margin_irf * u.deg diff --git a/gammapy/cube/make.py b/gammapy/cube/make.py index e699b6f2c1..3ac5c600ed 100644 --- a/gammapy/cube/make.py +++ b/gammapy/cube/make.py @@ -17,7 +17,6 @@ MARGIN_IRF_DEFAULT, MIGRA_AXIS_DEFAULT, RAD_AXIS_DEFAULT, - ENERGY_AXIS_DEFAULT, MapDataset, ) from .psf_map import make_psf_map @@ -72,7 +71,7 @@ def __init__( self.background_oversampling = background_oversampling self.migra_axis = migra_axis if migra_axis else MIGRA_AXIS_DEFAULT self.rad_axis = rad_axis if rad_axis else RAD_AXIS_DEFAULT - self.energy_axis_true = energy_axis_true or ENERGY_AXIS_DEFAULT + self.energy_axis_true = energy_axis_true or geom.get_axis_by_name("energy") self.binsz_irf = binsz_irf or BINSZ_IRF_DEFAULT self.margin_irf = margin_irf or MARGIN_IRF_DEFAULT @@ -94,7 +93,7 @@ def _cutout_geom(self, geom, observation): return geom @lazyproperty - def wcs_irf(self): + def geom_image_irf(self): """Spatial geometry of IRF Maps (`Geom`)""" wcs = self.geom.to_image() return WcsGeom.create( @@ -106,9 +105,9 @@ def wcs_irf(self): ) @lazyproperty - def geom_irf(self): + def geom_exposure_irf(self): """Geom of Exposure map associated with IRFs (`Geom`)""" - return self.wcs_irf.to_cube([self.energy_axis_true]) + return self.geom_image_irf.to_cube([self.energy_axis_true]) @lazyproperty def geom_exposure(self): @@ -119,13 +118,15 @@ def geom_exposure(self): @lazyproperty def geom_psf(self): """PSFMap geom (`Geom`)""" - geom_psf = self.wcs_irf.to_cube([self.rad_axis, self.energy_axis_true]) + geom_psf = self.geom_image_irf.to_cube([self.rad_axis, self.energy_axis_true]) return geom_psf @lazyproperty def geom_edisp(self): """EdispMap geom (`Geom`)""" - geom_edisp = self.wcs_irf.to_cube([self.migra_axis, self.energy_axis_true]) + geom_edisp = self.geom_image_irf.to_cube( + [self.migra_axis, self.energy_axis_true] + ) return geom_edisp def make_counts(self, observation): @@ -183,7 +184,7 @@ def make_exposure_irf(self, observation): Exposure map. """ - geom = self._cutout_geom(self.geom_irf, observation) + geom = self._cutout_geom(self.geom_exposure_irf, observation) exposure = make_map_exposure_true_energy( pointing=observation.pointing_radec, livetime=observation.observation_live_time_duration, @@ -314,7 +315,7 @@ def make_mask_safe_irf(self, observation): mask : `~numpy.ndarray` Mask """ - geom = self._cutout_geom(self.geom_irf, observation) + geom = self._cutout_geom(self.geom_exposure_irf, observation) offset = geom.separation(observation.pointing_radec) return offset >= self.offset_max diff --git a/gammapy/cube/tests/test_make.py b/gammapy/cube/tests/test_make.py index a69887b195..0ceba507e2 100644 --- a/gammapy/cube/tests/test_make.py +++ b/gammapy/cube/tests/test_make.py @@ -38,9 +38,7 @@ def geom(ebounds, binsz=0.5): { # Default, same e_true and reco "geom": geom(ebounds=[0.1, 1, 10]), - "e_true": MapAxis.from_edges( - [0.1, 1, 10], name="energy", unit="TeV", interp="log" - ), + "e_true": None, "counts": 34366, "exposure": 9.995376e08, "exposure_image": 7.921993e10, @@ -51,9 +49,7 @@ def geom(ebounds, binsz=0.5): { # Test single energy bin "geom": geom(ebounds=[0.1, 10]), - "e_true": MapAxis.from_edges( - [0.1, 10], name="energy", unit="TeV", interp="log" - ), + "e_true": None, "counts": 34366, "exposure": 5.843302e08, "exposure_image": 1.16866e11, @@ -64,9 +60,7 @@ def geom(ebounds, binsz=0.5): { # Test single energy bin with exclusion mask "geom": geom(ebounds=[0.1, 10]), - "e_true": MapAxis.from_edges( - [0.1, 10], name="energy", unit="TeV", interp="log" - ), + "e_true": None, "exclusion_mask": Map.from_geom(geom(ebounds=[0.1, 10])), "counts": 34366, "exposure": 5.843302e08, @@ -90,7 +84,7 @@ def geom(ebounds, binsz=0.5): "margin_irf": 0.0, }, { - # Test for different e_true and e_reco bins + # Test for different e_true and e_reco and spatial bins "geom": geom(ebounds=[0.1, 1, 10]), "e_true": MapAxis.from_edges( [0.1, 0.5, 2.5, 10.0], name="energy", unit="TeV", interp="log" From d662561da82c2a56af03363ac1d10bf6e3de3ad9 Mon Sep 17 00:00:00 2001 From: Axel Donath Date: Thu, 24 Oct 2019 10:26:20 +0200 Subject: [PATCH 10/11] Minor cleanup and fix of lightcurve notebook --- gammapy/cube/fit.py | 18 +++++++++--------- gammapy/scripts/analysis.py | 1 - tutorials/hess.ipynb | 10 ---------- tutorials/light_curve.ipynb | 15 +++------------ 4 files changed, 12 insertions(+), 32 deletions(-) diff --git a/gammapy/cube/fit.py b/gammapy/cube/fit.py index ce55ec1211..bea4807589 100644 --- a/gammapy/cube/fit.py +++ b/gammapy/cube/fit.py @@ -408,19 +408,19 @@ def create( margin_irf = margin_irf or MARGIN_IRF_DEFAULT margin_irf = margin_irf * u.deg - wcs = geom.to_image() - geom_exposure = wcs.to_cube([energy_axis_true]) + geom_image = geom.to_image() + geom_exposure = geom_image.to_cube([energy_axis_true]) - wcs_irf = WcsGeom.create( + geom_irf = WcsGeom.create( binsz=binsz_irf, - width=wcs.width + margin_irf, - skydir=wcs.center_skydir, - proj=wcs.projection, - coordsys=wcs.coordsys, + width=geom_image.width + margin_irf, + skydir=geom_image.center_skydir, + proj=geom_image.projection, + coordsys=geom_image.coordsys, ) - geom_psf = wcs_irf.to_cube([rad_axis, energy_axis_true]) - geom_edisp = wcs_irf.to_cube([migra_axis, energy_axis_true]) + geom_psf = geom_irf.to_cube([rad_axis, energy_axis_true]) + geom_edisp = geom_irf.to_cube([migra_axis, energy_axis_true]) return cls.from_geoms( geom, diff --git a/gammapy/scripts/analysis.py b/gammapy/scripts/analysis.py index 23d2417270..f593ff87b6 100644 --- a/gammapy/scripts/analysis.py +++ b/gammapy/scripts/analysis.py @@ -251,7 +251,6 @@ def _map_making(self): geom_irf["binsz_irf"] = self.settings["datasets"].get("binsz", None) geom_irf["margin_irf"] = self.settings["datasets"].get("margin", None) - offset_max = Angle(self.settings["datasets"]["offset-max"]) stack_datasets = self.settings["datasets"]["stack-datasets"] log.info("Creating datasets.") diff --git a/tutorials/hess.ipynb b/tutorials/hess.ipynb index 2271e740ff..5154ff78c1 100644 --- a/tutorials/hess.ipynb +++ b/tutorials/hess.ipynb @@ -102,16 +102,6 @@ " interp: log\n", " node_type: edges\n", " unit: TeV \n", - " energy-axis-true:\n", - " name: energy\n", - " hi_bnd: 10\n", - " lo_bnd: 1\n", - " nbin: 11\n", - " interp: log\n", - " node_type: edges\n", - " unit: TeV\n", - " binsz-irf: 0.2\n", - " margin-irf: 0.5\n", " \n", "fit:\n", " fit_range:\n", diff --git a/tutorials/light_curve.ipynb b/tutorials/light_curve.ipynb index eda89300a0..6e16599171 100644 --- a/tutorials/light_curve.ipynb +++ b/tutorials/light_curve.ipynb @@ -136,19 +136,10 @@ " axes=[energy_axis],\n", ")\n", "\n", - "etrue_axis = MapAxis.from_bounds(\n", + "energy_axis_true = MapAxis.from_bounds(\n", " 0.1, 20, 20, unit=\"TeV\", name=\"energy\", interp=\"log\"\n", ")\n", "\n", - "geom_true = WcsGeom.create(\n", - " skydir=target_position,\n", - " binsz=0.2,\n", - " width=(2, 2),\n", - " coordsys=\"CEL\",\n", - " proj=\"CAR\",\n", - " axes=[etrue_axis],\n", - ")\n", - "\n", "offset_max = 2 * u.deg" ] }, @@ -213,7 +204,7 @@ "\n", "datasets = []\n", "\n", - "maker = MapDatasetMaker(geom=geom, geom_true=geom_true, offset_max=offset_max)\n", + "maker = MapDatasetMaker(geom=geom, energy_axis_true=energy_axis_true, offset_max=offset_max)\n", "\n", "for time_interval in time_intervals:\n", " # get filtered observation lists in time interval\n", @@ -230,7 +221,7 @@ " )\n", " continue\n", "\n", - " stacked = MapDataset.create(geom=geom, geom_irf=geom_true)\n", + " stacked = MapDataset.create(geom=geom, energy_axis_true=energy_axis_true)\n", "\n", " for obs in observations:\n", " dataset = maker.run(obs)\n", From 24c8c0d170e45c66710c59f634f6e9cbb80d94c0 Mon Sep 17 00:00:00 2001 From: Axel Donath Date: Thu, 24 Oct 2019 12:14:35 +0200 Subject: [PATCH 11/11] Fix light-curve notebook --- tutorials/light_curve.ipynb | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/tutorials/light_curve.ipynb b/tutorials/light_curve.ipynb index 6e16599171..a39e9f3449 100644 --- a/tutorials/light_curve.ipynb +++ b/tutorials/light_curve.ipynb @@ -232,9 +232,8 @@ " position=target_position, e_reco=energy_axis.edges\n", " )\n", "\n", - " geom_psf = geom_true.to_binsz(binsz=0.02)\n", " stacked.psf = stacked.psf.get_psf_kernel(\n", - " position=target_position, geom=geom_psf, max_radius=\"0.3 deg\"\n", + " position=target_position, geom=stacked.exposure.geom, max_radius=\"0.3 deg\"\n", " )\n", "\n", " stacked.counts.meta[\"t_start\"] = time_interval[0]\n", @@ -420,7 +419,7 @@ " bkg_estimate=bkg_estimator.result,\n", " containment_correction=True,\n", " e_reco=energy_axis.edges,\n", - " e_true=etrue_axis.edges,\n", + " e_true=energy_axis_true.edges,\n", ")\n", "extraction.run()\n", "datasets_1d = extraction.spectrum_observations\n",