Add asymmetric energy dispersion/migration support

docs/changes/300.feature.rst

@@ -0,0 +1,11 @@
+Add functions to calculate asymmetric energy dispersion and energy migration matrices to ``pyirf.irf.energy_dispersion``, supporting both binning in ``true_source_fov_offset`` / ``true_source_fov_position_angle`` and ``true_source_fov_lon`` / ``true_source_fov_lat``.
+These functions are
+* ``energy_dispersion_asymmetric_polar``
+* ``energy_dispersion_asymmetric_lonlat``
+* ``energy_migration_matrix_asymmetric_polar``
+* ``energy_migration_matrix_asymmetric_lonlat``
+* ``energy_dispersion_to_migration_asymmetric``
+
+Also adds functions to create fits table HDUs in a format to the existing GADF definition for the energy dispersion to ``pyirf.io.gadf``.
+* ``create_energy_dispersion_asymmetric_polar_hdu``
+* ``create_energy_dispersion_asymmetric_lonlat_hdu``

pyirf/io/__init__.py

@@ -2,6 +2,8 @@
 from .gadf import (
     create_aeff2d_hdu,
     create_energy_dispersion_hdu,
+    create_energy_dispersion_3d_polar_hdu,
+    create_energy_dispersion_3d_lonlat_hdu,
     create_psf_table_hdu,
     create_rad_max_hdu,
     create_background_2d_hdu,
@@ -13,6 +15,8 @@
     "create_psf_table_hdu",
     "create_aeff2d_hdu",
     "create_energy_dispersion_hdu",
+    "create_energy_dispersion_3d_polar_hdu",
+    "create_energy_dispersion_3d_lonlat_hdu",
     "create_psf_table_hdu",
     "create_rad_max_hdu",
     "create_background_2d_hdu",

pyirf/io/gadf.py

@@ -11,6 +11,8 @@
 __all__ = [
     "create_aeff2d_hdu",
     "create_energy_dispersion_hdu",
+    "create_energy_dispersion_3d_polar_hdu",
+    "create_energy_dispersion_3d_lonlat_hdu",
     "create_psf_table_hdu",
     "create_rad_max_hdu",
 ]
@@ -202,6 +204,130 @@ def create_energy_dispersion_hdu(
     return BinTableHDU(edisp, header=header, name=extname)
 
 
+def create_energy_dispersion_3d_polar_hdu(
+    energy_dispersion,
+    true_energy_bins,
+    fov_offset_bins,
+    fov_position_angle_bins,
+    migration_bins,
+    point_like=True,
+    extname="EDISP",
+    **header_cards,
+):
+    """
+    Create a fits binary table HDU in a format similar to the GADF definition of the energy dispersion, but using two-axis polar FoV coordinates.
+    See the similar GADF specification at
+    https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/edisp/index.html
+
+    Parameters
+    ----------
+    energy_dispersion: numpy.ndarray
+        Energy dispersion array, must have shape
+        (n_energy_bins, n_migra_bins, n_source_offset_bins)
+    true_energy_bins: astropy.units.Quantity[energy]
+        Bin edges in true energy
+    fov_offset_bins: astropy.units.Quantity[angle]
+        Bin edges in the field of view offset.
+        For Point-Like IRFs, only giving a single bin is appropriate.
+    fov_position_angle_bins: astropy.units.Quantity[angle]
+        Bin edges in the field of view position angle.
+        For Point-Like IRFs, only giving a single bin is appropriate.
+    migration_bins: numpy.ndarray
+        Bin edges for the relative energy migration (``reco_energy / true_energy``)
+    point_like: bool
+        If the provided effective area was calculated after applying a direction cut,
+        pass ``True``, else ``False`` for a full-enclosure effective area.
+    extname: str
+        Name for BinTableHDU
+    **header_cards
+        Additional metadata to add to the header, use this to set e.g. TELESCOP or
+        INSTRUME.
+    """
+
+    edisp = QTable()
+    edisp["ENERG_LO"], edisp["ENERG_HI"] = binning.split_bin_lo_hi(true_energy_bins[np.newaxis, :].to(u.TeV))
+    edisp["THETA_LO"], edisp["THETA_HI"] = binning.split_bin_lo_hi(fov_offset_bins[np.newaxis, :].to(u.deg))
+    edisp["PHI_LO"], edisp["PHI_HI"] = binning.split_bin_lo_hi(fov_position_angle_bins[np.newaxis, :].to(u.deg))
+    edisp["MIGRA_LO"], edisp["MIGRA_HI"] = binning.split_bin_lo_hi(migration_bins[np.newaxis, :])
+    # transpose as FITS uses opposite dimension order
+    edisp["MATRIX"] = u.Quantity(energy_dispersion.T[np.newaxis, ...]).to(u.one)
+
+    # required header keywords
+    header = DEFAULT_HEADER.copy()
+    header["HDUCLAS1"] = "RESPONSE"
+    header["HDUCLAS2"] = "EDISP"
+    header["HDUCLAS3"] = "POINT-LIKE" if point_like else "FULL-ENCLOSURE"
+    header["HDUCLAS4"] = "EDISP_3D"
+    header["DATE"] = Time.now().utc.iso
+    idx = edisp.colnames.index("MATRIX") + 1
+    header[f"CREF{idx}"] = "(ENERG_LO:ENERG_HI,THETA_LO:THETA_HI,PHI_LO:PHI_HI,MIGRA_LO:MIGRA_HI)"
+    _add_header_cards(header, **header_cards)
+
+    return BinTableHDU(edisp, header=header, name=extname)
+
+
+def create_energy_dispersion_3d_lonlat_hdu(
+    energy_dispersion,
+    true_energy_bins,
+    fov_longitude_bins,
+    fov_latitude_bins,
+    migration_bins,
+    point_like=True,
+    extname="EDISP",
+    **header_cards,
+):
+    """
+    Create a fits binary table HDU in a format similar to the GADF definition of the energy dispersion, but using two-axis longitude-latitude FoV coordinates.
+    See the similar GADF specification at
+    https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/edisp/index.html
+
+    Parameters
+    ----------
+    energy_dispersion: numpy.ndarray
+        Energy dispersion array, must have shape
+        (n_energy_bins, n_migra_bins, n_source_offset_bins)
+    true_energy_bins: astropy.units.Quantity[energy]
+        Bin edges in true energy
+    fov_longitude_bins: astropy.units.Quantity[angle]
+        Bin edges in the field of view longitude.
+        For Point-Like IRFs, only giving a single bin is appropriate.
+    fov_latitiude_bins: astropy.units.Quantity[angle]
+        Bin edges in the field of view latitude.
+        For Point-Like IRFs, only giving a single bin is appropriate.
+    migration_bins: numpy.ndarray
+        Bin edges for the relative energy migration (``reco_energy / true_energy``)
+    point_like: bool
+        If the provided effective area was calculated after applying a direction cut,
+        pass ``True``, else ``False`` for a full-enclosure effective area.
+    extname: str
+        Name for BinTableHDU
+    **header_cards
+        Additional metadata to add to the header, use this to set e.g. TELESCOP or
+        INSTRUME.
+    """
+
+    edisp = QTable()
+    edisp["ENERG_LO"], edisp["ENERG_HI"] = binning.split_bin_lo_hi(true_energy_bins[np.newaxis, :].to(u.TeV))
+    edisp["DETX_LO"], edisp["DETX_HI"] = binning.split_bin_lo_hi(fov_longitude_bins[np.newaxis, :].to(u.deg))
+    edisp["DETY_LO"], edisp["DETY_HI"] = binning.split_bin_lo_hi(fov_latitude_bins[np.newaxis, :].to(u.deg))
+    edisp["MIGRA_LO"], edisp["MIGRA_HI"] = binning.split_bin_lo_hi(migration_bins[np.newaxis, :])
+    # transpose as FITS uses opposite dimension order
+    edisp["MATRIX"] = u.Quantity(energy_dispersion.T[np.newaxis, ...]).to(u.one)
+
+    # required header keywords
+    header = DEFAULT_HEADER.copy()
+    header["HDUCLAS1"] = "RESPONSE"
+    header["HDUCLAS2"] = "EDISP"
+    header["HDUCLAS3"] = "POINT-LIKE" if point_like else "FULL-ENCLOSURE"
+    header["HDUCLAS4"] = "EDISP_3D"
+    header["DATE"] = Time.now().utc.iso
+    idx = edisp.colnames.index("MATRIX") + 1
+    header[f"CREF{idx}"] = "(ENERG_LO:ENERG_HI,DETX_LO:DETX_HI,DETY_LO:DETY_HI,MIGRA_LO:MIGRA_HI)"
+    _add_header_cards(header, **header_cards)
+
+    return BinTableHDU(edisp, header=header, name=extname)
+
+
 #: Unit to store background rate in GADF format
 #:
 #: see https://github.com/open-gamma-ray-astro/gamma-astro-data-formats/issues/153

pyirf/io/tests/test_gadf.py

@@ -45,6 +45,44 @@ def edisp_hdus():
     return edisp, hdus
 
 
+@pytest.fixture
+def edisp_3d_polar_hdus():
+    from pyirf.io import create_energy_dispersion_3d_polar_hdu
+
+    edisp = np.zeros(
+        (len(e_bins) - 1, len(fov_bins) - 1, len(fov_bins) - 1, len(migra_bins) - 1)
+    )
+    edisp[..., 50] = 1.0
+
+    hdus = [
+        create_energy_dispersion_3d_polar_hdu(
+            edisp, e_bins, fov_bins, fov_bins, migra_bins, point_like=point_like
+        )
+        for point_like in [True, False]
+    ]
+
+    return edisp, hdus
+
+
+@pytest.fixture
+def edisp_3d_lonlat_hdus():
+    from pyirf.io import create_energy_dispersion_3d_lonlat_hdu
+
+    edisp = np.zeros(
+        (len(e_bins) - 1, len(fov_bins) - 1, len(fov_bins) - 1, len(migra_bins) - 1)
+    )
+    edisp[..., 50] = 1.0
+
+    hdus = [
+        create_energy_dispersion_3d_lonlat_hdu(
+            edisp, e_bins, fov_bins, fov_bins, migra_bins, point_like=point_like
+        )
+        for point_like in [True, False]
+    ]
+
+    return edisp, hdus
+
+
 @pytest.fixture
 def psf_hdu():
     from pyirf.io import create_psf_table_hdu
@@ -133,6 +171,26 @@ def test_energy_dispersion_schema(edisp_hdus):
         EDISP_2D.validate_hdu(hdu)
 
 
+@pytest.mark.xfail(reason="EDISP_3D not implemented yet", raises=ImportError)
+def test_energy_dispersion_schema_3d_polar(edisp_3d_polar_hdus):
+    from ogadf_schema.irfs import EDISP_3D
+
+    _, hdus = edisp_3d_polar_hdus
+
+    for hdu in hdus:
+        EDISP_3D.validate_hdu(hdu)
+
+
+@pytest.mark.xfail(reason="EDISP_3D not implemented yet", raises=ImportError)
+def test_energy_dispersion_schema_3d_lonlat(edisp_3d_lonlat_hdus):
+    from ogadf_schema.irfs import EDISP_3D
+
+    _, hdus = edisp_3d_lonlat_hdus
+
+    for hdu in hdus:
+        EDISP_3D.validate_hdu(hdu)
+
+
 def test_psf_table_gammapy(psf_hdu):
     '''Test our psf is readable by gammapy'''
     from gammapy.irf import PSF3D

pyirf/irf/__init__.py

@@ -5,7 +5,11 @@
     effective_area_3d_polar,
     effective_area_3d_lonlat,
 )
-from .energy_dispersion import energy_dispersion
+from .energy_dispersion import (
+    energy_dispersion,
+    energy_dispersion_3d_polar,
+    energy_dispersion_3d_lonlat,
+)
 from .psf import psf_table
 from .background import background_2d
 
@@ -16,6 +20,8 @@
     "effective_area_3d_polar",
     "effective_area_3d_lonlat",
     "energy_dispersion",
+    "energy_dispersion_3d_polar",
+    "energy_dispersion_3d_lonlat",
     "psf_table",
     "background_2d",
 ]