diff --git a/gammapy/data/data_store.py b/gammapy/data/data_store.py index 8da28ded5f..175c14f85c 100644 --- a/gammapy/data/data_store.py +++ b/gammapy/data/data_store.py @@ -12,10 +12,11 @@ from astropy.time import Time from astropy.coordinates import SkyCoord from ..utils.scripts import make_path +from ..utils.energy import Energy from .obs_table import ObservationTable from .hdu_index_table import HDUIndexTable from .utils import _earth_location_from_dict -from ..irf import EnergyDependentTablePSF +from ..irf import EnergyDependentTablePSF, IRFStacker __all__ = [ 'DataStore', @@ -654,12 +655,13 @@ def make_psf(self, position, energy=None, theta=None): if not theta: theta = self.psf.to_table_psf(theta=offset).offset - psf_value = self.psf.to_table_psf(theta=offset, offset=theta).evaluate(energy) + psf_value = self.psf.to_table_psf(theta=offset, offset=theta)\ + .evaluate(energy) arf = self.aeff.evaluate(offset=offset, energy=energy) exposure = arf * self.observation_live_time_duration - psf = EnergyDependentTablePSF(energy=energy, offset=theta, exposure=exposure, - psf_value=psf_value) + psf = EnergyDependentTablePSF(energy=energy, offset=theta, + exposure=exposure, psf_value=psf_value) return psf @@ -668,6 +670,7 @@ class ObservationList(UserList): Could be extended to hold a more generic class of observations """ + def __str__(self): s = self.__class__.__name__ + '\n' s += 'Number of observations: {}\n'.format(len(self)) @@ -676,7 +679,8 @@ def __str__(self): return s def make_psf(self, position, energy=None, theta=None): - """Make energy-dependent mean PSF for a given position and a set of observations. + """Make energy-dependent mean PSF for a given position and a set of + observations. Parameters ---------- @@ -684,10 +688,12 @@ def make_psf(self, position, energy=None, theta=None): Position at which to compute the PSF energy : `~astropy.units.Quantity` 1-dim energy array for the output PSF. - If none is given, the energy array of the PSF from the first observation is used. + If none is given, the energy array of the PSF from the first + observation is used. theta : `~astropy.coordinates.Angle` 1-dim offset array for the output PSF. - If none is given, the energy array of the PSF from the first observation is used. + If none is given, the energy array of the PSF from the first + observation is used. Returns ------- @@ -708,6 +714,57 @@ def make_psf(self, position, energy=None, theta=None): psf_value += psf.psf_value.T * psf.exposure psf_value /= exposure - psf_tot = EnergyDependentTablePSF(energy=energy, offset=theta, exposure=exposure, + psf_tot = EnergyDependentTablePSF(energy=energy, offset=theta, + exposure=exposure, psf_value=psf_value.T) return psf_tot + + def make_mean_edisp(self, position, e_true, e_reco, + low_reco_threshold=Energy(0.002, "TeV"), + high_reco_threshold=Energy(150, "TeV")): + """Make mean edisp for a given position and a set of observations. + + Compute the mean edisp of a set of observations j at a given position + + The stacking is implemented in :func:`~gammapy.irf.IRFStacker.stack_edisp` + + Parameters + ---------- + position : `~astropy.coordinates.SkyCoord` + Position at which to compute the mean EDISP + e_true : `~gammapy.utils.energy.EnergyBounds` + True energy axis + e_reco : `~gammapy.utils.energy.EnergyBounds` + Reconstructed energy axis + low_reco_threshold : `~gammapy.utils.energy.Energy` + low energy threshold in reco energy, default 0.002 TeV + high_reco_threshold : `~gammapy.utils.energy.Energy` + high energy threshold in reco energy , default 150 TeV + + Returns + ------- + stacked_edisp : `~gammapy.irf.EnergyDispersion` + Stacked EDISP for a set of observation + """ + + list_aeff = list() + list_edisp = list() + list_livetime = list() + list_low_threshold = [low_reco_threshold] * len(self) + list_high_threshold = [high_reco_threshold] * len(self) + for obs in self: + offset = position.separation(obs.pointing_radec) + list_aeff.append(obs.aeff.to_effective_area_table(offset, + energy=e_true)) + list_edisp.append(obs.edisp.to_energy_dispersion(offset, + e_reco=e_reco, + e_true=e_true)) + list_livetime.append(obs.observation_live_time_duration) + + irf_stack = IRFStacker(list_aeff=list_aeff, list_edisp=list_edisp, + list_livetime=list_livetime, + list_low_threshold=list_low_threshold, + list_high_threshold=list_high_threshold) + irf_stack.stack_edisp() + + return irf_stack.stacked_edisp diff --git a/gammapy/data/tests/test_data_store.py b/gammapy/data/tests/test_data_store.py index 5bb3ed366d..fcced6fd0c 100644 --- a/gammapy/data/tests/test_data_store.py +++ b/gammapy/data/tests/test_data_store.py @@ -1,13 +1,14 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst from __future__ import absolute_import, division, print_function, unicode_literals import numpy as np -from numpy.testing import assert_allclose +from numpy.testing import assert_allclose, assert_equal from astropy.tests.helper import pytest, assert_quantity_allclose from astropy.coordinates import Angle, SkyCoord from astropy.units import Quantity import astropy.units as u -from ...data import DataStore, DataManager +from ...data import DataStore, DataManager, ObservationList from ...utils.testing import data_manager, requires_data, requires_dependency +from ...utils.energy import Energy from ...datasets import gammapy_extra from ...utils.energy import EnergyBounds @@ -121,17 +122,26 @@ def test_data_summary(data_manager): @requires_data('gammapy-extra') @pytest.mark.parametrize("pars,result", [ (dict(energy=None, theta=None), - dict(energy_shape=18, theta_shape=300, psf_energy=2.5178505859375 * u.TeV, psf_theta=0.05 * u.deg, - psf_exposure=Quantity(6878545291473.34, "cm2 s"), psf_value=Quantity(205215.42446175334, "1/sr"))), + dict(energy_shape=18, theta_shape=300, psf_energy=2.5178505859375 * u.TeV, + psf_theta=0.05 * u.deg, + psf_exposure=Quantity(6878545291473.34, "cm2 s"), + psf_value=Quantity(205215.42446175334, "1/sr"))), (dict(energy=EnergyBounds.equal_log_spacing(1, 10, 100, "TeV"), theta=None), - dict(energy_shape=101, theta_shape=300, psf_energy=1.2589254117941673 * u.TeV, psf_theta=0.05 * u.deg, - psf_exposure=Quantity(4622187644084.735, "cm2 s"), psf_value=Quantity(119662.71915415104, "1/sr"))), + dict(energy_shape=101, theta_shape=300, + psf_energy=1.2589254117941673 * u.TeV, psf_theta=0.05 * u.deg, + psf_exposure=Quantity(4622187644084.735, "cm2 s"), + psf_value=Quantity(119662.71915415104, "1/sr"))), (dict(energy=None, theta=Angle(np.arange(0, 2, 0.002), 'deg')), - dict(energy_shape=18, theta_shape=1000, psf_energy=2.5178505859375 * u.TeV, psf_theta=0.02 * u.deg, - psf_exposure=Quantity(6878545291473.34, "cm2 s"), psf_value=Quantity(23082.369133891403, "1/sr"))), - (dict(energy=EnergyBounds.equal_log_spacing(1, 10, 100, "TeV"), theta=Angle(np.arange(0, 2, 0.002), 'deg')), - dict(energy_shape=101, theta_shape=1000, psf_energy=1.2589254117941673 * u.TeV, psf_theta=0.02 * u.deg, - psf_exposure=Quantity(4622187644084.735, "cm2 s"), psf_value=Quantity(27987.773313506143, "1/sr"))), + dict(energy_shape=18, theta_shape=1000, + psf_energy=2.5178505859375 * u.TeV, psf_theta=0.02 * u.deg, + psf_exposure=Quantity(6878545291473.34, "cm2 s"), + psf_value=Quantity(23082.369133891403, "1/sr"))), + (dict(energy=EnergyBounds.equal_log_spacing(1, 10, 100, "TeV"), + theta=Angle(np.arange(0, 2, 0.002), 'deg')), + dict(energy_shape=101, theta_shape=1000, + psf_energy=1.2589254117941673 * u.TeV, psf_theta=0.02 * u.deg, + psf_exposure=Quantity(4622187644084.735, "cm2 s"), + psf_value=Quantity(27987.773313506143, "1/sr"))), ]) def test_make_psf(pars, result): position = SkyCoord(83.63, 22.01, unit='deg') @@ -139,14 +149,52 @@ def test_make_psf(pars, result): data_store = DataStore.from_dir(store) obs1 = data_store.obs(23523) - psf = obs1.make_psf(position=position, energy=pars["energy"], theta=pars["theta"]) + psf = obs1.make_psf(position=position, energy=pars["energy"], + theta=pars["theta"]) assert_allclose(psf.offset.shape, result["theta_shape"]) assert_allclose(psf.energy.shape, result["energy_shape"]) assert_allclose(psf.exposure.shape, result["energy_shape"]) - assert_allclose(psf.psf_value.shape, (result["energy_shape"], result["theta_shape"])) + assert_allclose(psf.psf_value.shape, (result["energy_shape"], + result["theta_shape"])) assert_quantity_allclose(psf.offset[10], result["psf_theta"]) assert_quantity_allclose(psf.energy[10], result["psf_energy"]) assert_quantity_allclose(psf.exposure[10], result["psf_exposure"]) assert_quantity_allclose(psf.psf_value[10, 50], result["psf_value"]) + + +@requires_dependency('scipy') +@requires_data('gammapy-extra') +def test_make_mean_edisp(tmpdir): + position = SkyCoord(83.63, 22.01, unit='deg') + store = gammapy_extra.filename("datasets/hess-crab4-hd-hap-prod2") + data_store = DataStore.from_dir(store) + + obs1 = data_store.obs(23523) + obs2 = data_store.obs(23592) + obslist = ObservationList([obs1, obs2]) + + e_true = EnergyBounds.equal_log_spacing(0.01, 150, 80, "TeV") + e_reco = EnergyBounds.equal_log_spacing(0.5, 100, 15, "TeV") + rmf = obslist.make_mean_edisp(position=position, e_true=e_true, + e_reco=e_reco) + + assert len(rmf.e_true.nodes) == 80 + assert len(rmf.e_reco.nodes) == 15 + assert_quantity_allclose(rmf.data[53, 8], 0.0559785805550798) + + rmf2 = obslist.make_mean_edisp(position=position, e_true=e_true, + e_reco=e_reco, + low_reco_threshold=Energy(1, "TeV"), + high_reco_threshold=Energy(60, "TeV")) + i2 = np.where(rmf2.evaluate(e_reco=Energy(0.8, "TeV")) != 0)[0] + assert len(i2) == 0 + i2 = np.where(rmf2.evaluate(e_reco=Energy(61, "TeV")) != 0)[0] + assert len(i2) == 0 + i = np.where(rmf.evaluate(e_reco=Energy(1.5, "TeV")) != 0)[0] + i2 = np.where(rmf2.evaluate(e_reco=Energy(1.5, "TeV")) != 0)[0] + assert_equal(i, i2) + i = np.where(rmf.evaluate(e_reco=Energy(55, "TeV")) != 0)[0] + i2 = np.where(rmf2.evaluate(e_reco=Energy(55, "TeV")) != 0)[0] + assert_equal(i, i2) diff --git a/gammapy/irf/irf_stack.py b/gammapy/irf/irf_stack.py index ad0692f8f3..15224aa6c3 100644 --- a/gammapy/irf/irf_stack.py +++ b/gammapy/irf/irf_stack.py @@ -16,7 +16,7 @@ class IRFStacker(object): r""" Stack instrument response functions - Compute mean effective area and the mean energy dispersio for a given for a + Compute mean effective area and the mean energy dispersion for a given for a given list of instrument response functions. Results are stored as attributes. @@ -39,11 +39,11 @@ class IRFStacker(object): Parameters ---------- - list_arf: list + list_aeff: list list of `~gammapy.irf.EffectiveAreaTable` list_livetime: list list of `~astropy.units.Quantity` (livetime) - list_rmf: list + list_edisp: list list of `~gammapy.irf.EnergyDispersion` list_low_threshold: list list of low energy threshold, optional for effective area mean computation @@ -53,11 +53,11 @@ class IRFStacker(object): """ - def __init__(self, list_arf, list_livetime, list_rmf=None, + def __init__(self, list_aeff, list_livetime, list_edisp=None, list_low_threshold=None, list_high_threshold=None): - self.list_arf = list_arf + self.list_aeff = list_aeff self.list_livetime = Quantity(list_livetime) - self.list_rmf = list_rmf + self.list_edisp = list_edisp self.list_low_threshold = list_low_threshold self.list_high_threshold = list_high_threshold self.stacked_aeff = None @@ -67,40 +67,40 @@ def stack_aeff(self): """ Compute mean effective area """ - nbins = self.list_arf[0].energy.nbins + nbins = self.list_aeff[0].energy.nbins aefft = Quantity(np.zeros(nbins), 'cm2 s') livetime_tot = np.sum(self.list_livetime) - for i, arf in enumerate(self.list_arf): - aeff_data = arf.evaluate(fill_nan=True) + for i, aeff in enumerate(self.list_aeff): + aeff_data = aeff.evaluate(fill_nan=True) aefft_current = aeff_data * self.list_livetime[i] aefft += aefft_current stacked_data = aefft / livetime_tot - self.stacked_aeff = EffectiveAreaTable(energy=self.list_arf[0].energy, + self.stacked_aeff = EffectiveAreaTable(energy=self.list_aeff[0].energy, data=stacked_data.to('cm2')) - def mean_rmf(self): + def stack_edisp(self): """ Compute mean energy dispersion """ - reco_bins = self.list_rmf[0].e_reco.nbins - true_bins = self.list_rmf[0].e_true.nbins + reco_bins = self.list_edisp[0].e_reco.nbins + true_bins = self.list_edisp[0].e_true.nbins aefft = Quantity(np.zeros(true_bins), 'cm2 s') temp = np.zeros(shape=(reco_bins, true_bins)) aefftedisp = Quantity(temp, 'cm2 s') - for i, rmf in enumerate(self.list_rmf): - aeff_data = self.list_arf[i].evaluate(fill_nan=True) + for i, edisp in enumerate(self.list_edisp): + aeff_data = self.list_aeff[i].evaluate(fill_nan=True) aefft_current = aeff_data * self.list_livetime[i] aefft += aefft_current - edisp_data = rmf.pdf_in_safe_range(self.list_low_threshold[i], - self.list_high_threshold[i]) + edisp_data = edisp.pdf_in_safe_range(self.list_low_threshold[i], + self.list_high_threshold[i]) aefftedisp += edisp_data.transpose() * aefft_current stacked_edisp = np.nan_to_num(aefftedisp / aefft) - self.stacked_edisp = EnergyDispersion(e_true=self.list_rmf[0].e_true, - e_reco=self.list_rmf[0].e_reco, + self.stacked_edisp = EnergyDispersion(e_true=self.list_edisp[0].e_true, + e_reco=self.list_edisp[0].e_reco, data=stacked_edisp.transpose()) diff --git a/gammapy/spectrum/observation.py b/gammapy/spectrum/observation.py index 9845853e98..722f0abe4e 100644 --- a/gammapy/spectrum/observation.py +++ b/gammapy/spectrum/observation.py @@ -633,7 +633,7 @@ def stack_aeff(self): for o in self.obs_list: list_arf.append(o.aeff) list_livetime.append(o.livetime) - irf_stack = IRFStacker(list_arf=list_arf, list_livetime=list_livetime) + irf_stack = IRFStacker(list_aeff=list_arf, list_livetime=list_livetime) irf_stack.stack_aeff() self.stacked_aeff = irf_stack.stacked_aeff @@ -644,22 +644,22 @@ def stack_edisp(self): calls :func:`~gammapy.irf.IRFStacker.stack_edisp` """ list_arf = list() - list_rmf = list() + list_edisp = list() list_livetime = list() list_elo_threshold = list() list_ehi_threshold = list() for o in self.obs_list: list_arf.append(o.aeff) list_livetime.append(o.livetime) - list_rmf.append(o.edisp) + list_edisp.append(o.edisp) list_elo_threshold.append(o.lo_threshold) list_ehi_threshold.append(o.hi_threshold) - irf_stack = IRFStacker(list_arf=list_arf, + irf_stack = IRFStacker(list_aeff=list_arf, list_livetime=list_livetime, - list_rmf=list_rmf, + list_edisp=list_edisp, list_low_threshold=list_elo_threshold, list_high_threshold=list_ehi_threshold) - irf_stack.mean_rmf() + irf_stack.stack_edisp() self.stacked_edisp = irf_stack.stacked_edisp def stack_obs(self):