Skip to content

Commit

Permalink
Merge pull request #5186 from Astro-Kirsty/mapaxis_edisp
Browse files Browse the repository at this point in the history
Adjust `to_edisp_kernel` input to `MapAxis`
  • Loading branch information
registerrier committed Jun 5, 2024
2 parents a8de5e6 + e8e4a7c commit 619f6f4
Show file tree
Hide file tree
Showing 4 changed files with 39 additions and 19 deletions.
6 changes: 3 additions & 3 deletions gammapy/datasets/tests/test_map.py
Original file line number Diff line number Diff line change
Expand Up @@ -164,10 +164,10 @@ def get_psf():
def get_edisp(geom, geom_etrue):
filename = "$GAMMAPY_DATA/hess-dl3-dr1/data/hess_dl3_dr1_obs_id_020136.fits.gz"
edisp2d = EnergyDispersion2D.read(filename, hdu="EDISP")
energy = geom.axes["energy"].edges
energy_true = geom_etrue.axes["energy_true"].edges
energy_axis = geom.axes["energy"]
energy_axis_true = geom_etrue.axes["energy_true"]
edisp_kernel = edisp2d.to_edisp_kernel(
offset="1.2 deg", energy=energy, energy_true=energy_true
offset="1.2 deg", energy_axis=energy_axis, energy_axis_true=energy_axis_true
)
edisp = EDispKernelMap.from_edisp_kernel(edisp_kernel)
return edisp
Expand Down
37 changes: 25 additions & 12 deletions gammapy/irf/edisp/core.py
Original file line number Diff line number Diff line change
@@ -1,18 +1,24 @@
# Licensed under a 3-clause BSD style license - see LICENSE.rst
import logging
import numpy as np
import scipy.special
from astropy import units as u
from astropy.coordinates import Angle, SkyCoord
from astropy.units import Quantity
from astropy.visualization import quantity_support
import matplotlib.pyplot as plt
from matplotlib.colors import PowerNorm
from gammapy.maps import MapAxes, MapAxis, RegionGeom
from gammapy.utils.deprecation import deprecated_renamed_argument
from gammapy.visualization.utils import add_colorbar
from ..core import IRF

__all__ = ["EnergyDispersion2D"]


log = logging.getLogger(__name__)


class EnergyDispersion2D(IRF):
"""Offset-dependent energy dispersion matrix.
Expand Down Expand Up @@ -41,8 +47,9 @@ class EnergyDispersion2D(IRF):
Create energy dispersion matrix (`~gammapy.irf.EnergyDispersion`)
for a given field of view offset and energy binning:
>>> energy = MapAxis.from_bounds(0.1, 20, nbin=60, unit="TeV", interp="log").edges
>>> edisp = edisp2d.to_edisp_kernel(offset='1.2 deg', energy=energy, energy_true=energy)
>>> energy_axis = MapAxis.from_bounds(0.1, 20, nbin=60, unit="TeV", interp="log", name='energy')
>>> edisp = edisp2d.to_edisp_kernel(offset='1.2 deg', energy_axis=energy_axis,
... energy_axis_true=energy_axis.copy(name='energy_true'))
See Also
--------
Expand Down Expand Up @@ -117,7 +124,13 @@ def from_gauss(
data=data.value,
)

def to_edisp_kernel(self, offset, energy_true=None, energy=None):
@deprecated_renamed_argument(
["energy_true", "energy"],
["energy_axis_true", "energy_axis"],
["v1.3", "v1.3"],
arg_in_kwargs=True,
)
def to_edisp_kernel(self, offset, energy_axis_true=None, energy_axis=None):
"""Detector response R(Delta E_reco, Delta E_true).
Probability to reconstruct an energy in a given true energy band
Expand All @@ -127,9 +140,9 @@ def to_edisp_kernel(self, offset, energy_true=None, energy=None):
----------
offset : `~astropy.coordinates.Angle`
Offset.
energy_true : `~astropy.units.Quantity`, optional
energy_axis_true : `~gammapy.maps.MapAxis`, optional
True energy axis. Default is None.
energy : `~astropy.units.Quantity`, optional
energy_axis : `~gammapy.maps.MapAxis`, optional
Reconstructed energy axis. Default is None.
Returns
Expand All @@ -141,18 +154,18 @@ def to_edisp_kernel(self, offset, energy_true=None, energy=None):

offset = Angle(offset)

if energy is None:
if isinstance(energy_axis, Quantity):
energy_axis = MapAxis.from_energy_edges(energy_axis)
if energy_axis is None:
energy_axis = self.axes["energy_true"].copy(name="energy")
else:
energy_axis = MapAxis.from_energy_edges(energy)

if energy_true is None:
energy_axis_true = self.axes["energy_true"]
else:
if isinstance(energy_axis_true, Quantity):
energy_axis_true = MapAxis.from_energy_edges(
energy_true,
energy_axis_true,
name="energy_true",
)
if energy_axis_true is None:
energy_axis_true = self.axes["energy_true"]

pointing = SkyCoord("0d", "0d")

Expand Down
2 changes: 1 addition & 1 deletion gammapy/irf/edisp/kernel.py
Original file line number Diff line number Diff line change
Expand Up @@ -150,7 +150,7 @@ def from_gauss(cls, energy_axis_true, energy_axis, sigma, bias, pdf_threshold=1e
pdf_threshold=pdf_threshold,
)
return edisp.to_edisp_kernel(
offset=offset_axis.center[0], energy=energy_axis.edges
offset=offset_axis.center[0], energy_axis=energy_axis
)

@classmethod
Expand Down
13 changes: 10 additions & 3 deletions gammapy/irf/edisp/tests/test_core.py
Original file line number Diff line number Diff line change
Expand Up @@ -66,10 +66,17 @@ def test_evaluation(self):
def test_exporter(self):
# Check RMF exporter
offset = Angle(0.612, "deg")
e_reco = MapAxis.from_energy_bounds(1, 10, 7, "TeV").edges
e_true = MapAxis.from_energy_bounds(0.8, 5, 5, "TeV").edges
rmf = self.edisp.to_edisp_kernel(offset, energy_true=e_true, energy=e_reco)
e_reco_axis = MapAxis.from_energy_bounds(1, 10, 7, "TeV")
e_true_axis = MapAxis.from_energy_bounds(0.8, 5, 5, "TeV", name="energy_true")
rmf = self.edisp.to_edisp_kernel(
offset, energy_axis_true=e_true_axis, energy_axis=e_reco_axis
)
assert_allclose(rmf.data.data[2, 3], 0.08, atol=5e-2) # same tolerance as above
# Test for v1.3
rmf2 = self.edisp.to_edisp_kernel(
offset, energy_axis_true=e_true_axis.edges, energy_axis=e_reco_axis.edges
)
assert_allclose(rmf2.data.data[2, 3], 0.08, atol=5e-2)

def test_write(self):
energy_axis_true = MapAxis.from_energy_bounds(
Expand Down

0 comments on commit 619f6f4

Please sign in to comment.