Repository: cta-observatory/pyirf

Issue:
Implement Background rate IRF (symmetric case)
Description here: https://open-gamma-ray-astro.github.io/gamma-astro-data-formats/irfs/full_enclosure/bkg/index.html

Discussion:

Patch:
diff --git a/docs/notebooks/comparison_with_EventDisplay.ipynb b/docs/notebooks/comparison_with_EventDisplay.ipynb
index 5df5c1a..46ccca9 100644
--- a/docs/notebooks/comparison_with_EventDisplay.ipynb
+++ b/docs/notebooks/comparison_with_EventDisplay.ipynb
@@ -593,6 +593,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "from pyirf.utils import cone_solid_angle\n",
+    "\n",
     "# Data from EventDisplay\n",
     "h = irf_eventdisplay[\"BGRate\"]\n",
     "x = 0.5 * (10**h.edges[:-1] + 10**h.edges[1:])\n",
@@ -600,15 +602,45 @@
     "y = h.values\n",
     "yerr = np.sqrt(h.variances)\n",
     "\n",
+    "# pyirf data\n",
+    "\n",
+    "bg_rate = QTable.read(pyirf_file, hdu='BACKGROUND')[0]\n",
+    "\n",
+    "reco_bins = np.append(bg_rate['ENERG_LO'], bg_rate['ENERG_HI'][-1])\n",
+    "\n",
+    "# first fov bin, [0, 1] deg\n",
+    "fov_bin = 0\n",
+    "rate_bin = bg_rate['BKG'].T[:, fov_bin]\n",
+    "\n",
+    "\n",
+    "# interpolate theta cut for given e reco bin\n",
+    "e_center_bg = 0.5 * (bg_rate['ENERG_LO'] + bg_rate['ENERG_HI'])\n",
+    "e_center_theta = 0.5 * (rad_max['ENERG_LO'] + rad_max['ENERG_HI'])\n",
+    "theta_cut = np.interp(e_center_bg, e_center_theta, rad_max['RAD_MAX'].T[:, 0])\n",
+    "\n",
+    "# undo normalization\n",
+    "rate_bin *= cone_solid_angle(theta_cut)\n",
+    "rate_bin *= np.diff(reco_bins)\n",
+    "\n",
+    "\n",
+    "# Plot function\n",
+    "plt.errorbar(x, y, xerr=xerr, yerr=yerr, ls='', label=\"EventDisplay\")\n",
+    "\n",
+    "plt.errorbar(\n",
+    "    0.5 * (bg_rate['ENERG_LO'] + bg_rate['ENERG_HI']).to_value(u.TeV)[1:-1],\n",
+    "    rate_bin.to_value(u.s**-1)[1:-1],\n",
+    "    xerr=np.diff(reco_bins).to_value(u.TeV)[1:-1] / 2,\n",
+    "    ls='',\n",
+    "    label='pyirf',\n",
+    ")\n",
+    "\n",
     "# Style settings\n",
     "plt.xscale(\"log\")\n",
     "plt.xlabel(r\"$E_\\mathrm{Reco} / \\mathrm{TeV}$\")\n",
     "plt.ylabel(\"Background rate / s⁻¹\")\n",
     "plt.grid(which=\"both\")\n",
-    "\n",
-    "# Plot function\n",
-    "plt.errorbar(x, y, xerr=xerr, yerr=yerr, fmt=\"o\", label=\"EventDisplay\")\n",
     "plt.legend(loc=\"best\")\n",
+    "plt.yscale('log')\n",
     "\n",
     "None # to remove clutter by mpl objects"
    ]
diff --git a/examples/calculate_eventdisplay_irfs.py b/examples/calculate_eventdisplay_irfs.py
index 94208ce..8bea478 100644
--- a/examples/calculate_eventdisplay_irfs.py
+++ b/examples/calculate_eventdisplay_irfs.py
@@ -36,6 +36,7 @@ from pyirf.irf import (
     point_like_effective_area,
     energy_dispersion,
     psf_table,
+    background_2d,
 )
 
 from pyirf.io import (
@@ -43,6 +44,7 @@ from pyirf.io import (
     create_psf_table_hdu,
     create_energy_dispersion_hdu,
     create_rad_max_hdu,
+    create_background_2d_hdu,
 )
 
 
@@ -220,6 +222,9 @@ def main():
     true_energy_bins = add_overflow_bins(
         create_bins_per_decade(10 ** -1.9 * u.TeV, 10 ** 2.31 * u.TeV, 10,)
     )
+    reco_energy_bins = add_overflow_bins(
+        create_bins_per_decade(10 ** -1.9 * u.TeV, 10 ** 2.31 * u.TeV, 10,)
+    )
     fov_offset_bins = [0, 0.5] * u.deg
     source_offset_bins = np.arange(0, 1 + 1e-4, 1e-3) * u.deg
     energy_migration_bins = np.geomspace(0.2, 5, 200)
@@ -258,7 +263,6 @@ def main():
         gammas[gammas["selected"]], true_energy_bins,
     )
     ang_res = angular_resolution(gammas[gammas["selected_gh"]], true_energy_bins,)
-
     psf = psf_table(
         gammas[gammas["selected_gh"]],
         true_energy_bins,
@@ -266,16 +270,24 @@ def main():
         source_offset_bins=source_offset_bins,
     )
 
-    hdus.append(
-        create_psf_table_hdu(
-            psf, true_energy_bins, source_offset_bins, fov_offset_bins,
-        )
+    background_rate = background_2d(
+        background[background['selected_gh']],
+        reco_energy_bins,
+        fov_offset_bins=np.arange(0, 11) * u.deg,
+        t_obs=T_OBS,
     )
-    hdus.append(
-        create_rad_max_hdu(
+
+    hdus.append(create_background_2d_hdu(
+        background_rate,
+        reco_energy_bins,
+        fov_offset_bins=np.arange(0, 11) * u.deg,
+    ))
+    hdus.append(create_psf_table_hdu(
+            psf, true_energy_bins, source_offset_bins, fov_offset_bins,
+    ))
+    hdus.append(create_rad_max_hdu(
             theta_bins, fov_offset_bins, rad_max=theta_cuts_opt["cut"][:, np.newaxis]
-        )
-    )
+    ))
     hdus.append(fits.BinTableHDU(ang_res, name="ANGULAR_RESOLUTION"))
     hdus.append(fits.BinTableHDU(bias_resolution, name="ENERGY_BIAS_RESOLUTION"))
     fits.HDUList(hdus).writeto("pyirf_eventdisplay.fits.gz", overwrite=True)
diff --git a/pyirf/io/__init__.py b/pyirf/io/__init__.py
index 48dc609..aeefad5 100644
--- a/pyirf/io/__init__.py
+++ b/pyirf/io/__init__.py
@@ -4,6 +4,7 @@ from .gadf import (
     create_energy_dispersion_hdu,
     create_psf_table_hdu,
     create_rad_max_hdu,
+    create_background_2d_hdu,
 )
 
 
@@ -14,4 +15,5 @@ __all__ = [
     "create_energy_dispersion_hdu",
     "create_psf_table_hdu",
     "create_rad_max_hdu",
+    "create_background_2d_hdu",
 ]
diff --git a/pyirf/io/gadf.py b/pyirf/io/gadf.py
index 16c4a78..5563e1f 100644
--- a/pyirf/io/gadf.py
+++ b/pyirf/io/gadf.py
@@ -212,10 +212,73 @@ def create_energy_dispersion_hdu(
     return BinTableHDU(psf, 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
+#: for a discussion on why this is MeV not TeV as everywhere else
+GADF_BACKGROUND_UNIT = u.Unit('MeV-1 s-1 sr-1')
+
+
 @u.quantity_input(
-    psf=u.sr ** -1,
-    true_energy_bins=u.TeV,
+    background=GADF_BACKGROUND_UNIT,
+    reco_energy_bins=u.TeV,
+    fov_offset_bins=u.deg,
+)
+def create_background_2d_hdu(
+    background_2d,
+    reco_energy_bins,
+    fov_offset_bins,
+    extname="BACKGROUND",
+    **header_cards,
+):
+    """
+    Create a fits binary table HDU in GADF format for the background 2d table.
+    See the specification at
+    https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/bkg/index.html#bkg-2d
+
+    Parameters
+    ----------
+    background_2d: astropy.units.Quantity[(MeV s sr)^-1]
+        Background rate, must have shape
+        (n_energy_bins, n_fov_offset_bins)
+    reco_energy_bins: astropy.units.Quantity[energy]
+        Bin edges in reconstructed energy
+    fov_offset_bins: astropy.units.Quantity[angle]
+        Bin edges in the field of view offset.
+    extname: str
+        Name for BinTableHDU
+    **header_cards
+        Additional metadata to add to the header, use this to set e.g. TELESCOP or
+        INSTRUME.
+    """
+
+    bkg = QTable(
+        {
+            "ENERG_LO": u.Quantity(reco_energy_bins[:-1], ndmin=2).to(u.TeV),
+            "ENERG_HI": u.Quantity(reco_energy_bins[1:], ndmin=2).to(u.TeV),
+            "THETA_LO": u.Quantity(fov_offset_bins[:-1], ndmin=2).to(u.deg),
+            "THETA_HI": u.Quantity(fov_offset_bins[1:], ndmin=2).to(u.deg),
+            # transpose as FITS uses opposite dimension order
+            "BKG": background_2d.T[np.newaxis, ...].to(GADF_BACKGROUND_UNIT),
+        }
+    )
+
+    # required header keywords
+    header = DEFAULT_HEADER.copy()
+    header["HDUCLAS1"] = "RESPONSE"
+    header["HDUCLAS2"] = "BKG"
+    header["HDUCLAS3"] = "FULL-ENCLOSURE"
+    header["HDUCLAS4"] = "BKG_2D"
+    header["DATE"] = Time.now().utc.iso
+    _add_header_cards(header, **header_cards)
+
+    return BinTableHDU(bkg, header=header, name=extname)
+
+
+@u.quantity_input(
+    reco_energy_bins=u.TeV,
     fov_offset_bins=u.deg,
+    rad_max=u.deg,
     source_offset_bins=u.deg,
 )
 def create_rad_max_hdu(
diff --git a/pyirf/irf/__init__.py b/pyirf/irf/__init__.py
index 6a5e7df..a350667 100644
--- a/pyirf/irf/__init__.py
+++ b/pyirf/irf/__init__.py
@@ -1,10 +1,12 @@
 from .effective_area import effective_area, point_like_effective_area
 from .energy_dispersion import energy_dispersion
 from .psf import psf_table
+from .background import background_2d
 
 __all__ = [
     "effective_area",
     "point_like_effective_area",
     "energy_dispersion",
     "psf_table",
+    "background_2d",
 ]
diff --git a/pyirf/irf/background.py b/pyirf/irf/background.py
new file mode 100644
index 0000000..27209f0
--- /dev/null
+++ b/pyirf/irf/background.py
@@ -0,0 +1,58 @@
+import astropy.units as u
+import numpy as np
+
+from ..utils import cone_solid_angle
+
+#: Unit of the background rate IRF
+BACKGROUND_UNIT = u.Unit('s-1 TeV-1 sr-1')
+
+
+def background_2d(events, reco_energy_bins, fov_offset_bins, t_obs):
+    """
+    Calculate background rates in radially symmetric bins in the field of view.
+
+    GADF documentation here:
+    https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/bkg/index.html#bkg-2d
+
+    Parameters
+    ----------
+    events: astropy.table.QTable
+        DL2 events table of the selected background events.
+        Needed columns for this function: `source_fov_offset`, `reco_energy`, `weight`
+    reco_energy: astropy.units.Quantity[energy]
+        The bins in reconstructed energy to be used for the IRF
+    fov_offset_bins: astropy.units.Quantity[angle]
+        The bins in the field of view offset to be used for the IRF
+    t_obs: astropy.units.Quantity[time]
+        Observation time. This must match with how the individual event
+        weights are calculated.
+
+    Returns
+    -------
+    bg_rate: astropy.units.Quantity
+        The background rate as particles per energy, time and solid angle
+        in the specified bins.
+
+        Shape: (len(reco_energy_bins) - 1, len(fov_offset_bins) - 1)
+    """
+
+    hist, _, _ = np.histogram2d(
+        events["reco_energy"].to_value(u.TeV),
+        events["source_fov_offset"].to_value(u.deg),
+        bins=[
+            reco_energy_bins.to_value(u.TeV),
+            fov_offset_bins.to_value(u.deg),
+        ],
+        weights=events['weight'],
+    )
+
+    # divide all energy bins by their width
+    # hist has shape (n_energy, n_fov_offset) so we need to transpose and then back
+    bin_width_energy = np.diff(reco_energy_bins)
+    per_energy = (hist.T / bin_width_energy).T
+
+    # divide by solid angle in each fov bin and the observation time
+    bin_solid_angle = np.diff(cone_solid_angle(fov_offset_bins))
+    bg_rate = per_energy / t_obs / bin_solid_angle
+
+    return bg_rate.to(BACKGROUND_UNIT)


Test patch:
diff --git a/pyirf/io/tests/test_gadf.py b/pyirf/io/tests/test_gadf.py
index 0ae43eb..a049947 100644
--- a/pyirf/io/tests/test_gadf.py
+++ b/pyirf/io/tests/test_gadf.py
@@ -69,7 +69,6 @@ def test_psf_table():
     psf[:, 0, :] = 1
     psf = psf / cone_solid_angle(source_bins[1])
 
-
     for point_like in [True, False]:
         with tempfile.NamedTemporaryFile(suffix='.fits') as f:
             hdu = create_psf_table_hdu(
@@ -85,3 +84,33 @@ def test_psf_table():
             # unlike how it handles effective area and edisp
             # see https://github.com/gammapy/gammapy/issues/3025
             assert u.allclose(psf, psf3d.psf_value.T, atol=1e-16 / u.sr)
+
+
+# gammapy uses inconsistent axis order, should be fixed before gammapy 1.0
+# see https://github.com/gammapy/gammapy/issues/2067
+# TODO: remove xfail when this is fixed in gammapy and bump the required gammapy version
+@pytest.mark.xfail
+def test_background_2d():
+    '''Test our background hdu is readable by gammapy'''
+
+    pytest.importorskip('gammapy')
+    from pyirf.io import create_background_2d_hdu
+    from gammapy.irf import Background2D
+
+    e_bins = np.geomspace(0.1, 100, 31) * u.TeV
+    fov_bins = [0, 1, 2, 3] * u.deg
+    background = np.column_stack([
+        np.geomspace(1e9, 1e3, 3),
+        np.geomspace(0.5e9, 0.5e3, 3),
+        np.geomspace(1e8, 1e2, 3),
+    ]) * u.Unit('TeV-1 s-1 sr-1')
+
+    with tempfile.NamedTemporaryFile(suffix='.fits') as f:
+        hdu = create_background_2d_hdu(background, e_bins, fov_bins)
+
+        fits.HDUList([fits.PrimaryHDU(), hdu]).writeto(f.name)
+
+        # test reading with gammapy works
+        bg2d = Background2D.read(f.name, 'BACKGROUND')
+
+        assert u.allclose(background, bg2d.data.data, atol=1e-16 * u.Unit('TeV-1 s-1 sr-1'))
diff --git a/pyirf/irf/tests/test_background.py b/pyirf/irf/tests/test_background.py
new file mode 100644
index 0000000..c476850
--- /dev/null
+++ b/pyirf/irf/tests/test_background.py
@@ -0,0 +1,39 @@
+import astropy.units as u
+from astropy.table import QTable
+import numpy as np
+
+
+def test_background():
+    from pyirf.irf import background_2d
+    from pyirf.utils import cone_solid_angle
+
+    np.random.seed(0)
+
+    N1 = 1000
+    N2 = 100
+    N = N1 + N2
+
+    # toy event data set with just two energies
+    # and a psf per energy bin, point-like
+    events = QTable(
+        {
+            "reco_energy": np.append(np.full(N1, 1), np.full(N2, 2)) * u.TeV,
+            "source_fov_offset": np.zeros(N) * u.deg,
+            "weight": np.ones(N),
+        }
+    )
+
+    energy_bins = [0, 1.5, 3] * u.TeV
+    fov_bins = [0, 1] * u.deg
+
+    # We return a table with one row as needed for gadf
+    bg = background_2d(events, energy_bins, fov_bins, t_obs=1 * u.s)
+
+    # 2 energy bins, 1 fov bin, 200 source distance bins
+    assert bg.shape == (2, 1)
+    assert bg.unit == u.Unit("TeV-1 s-1 sr-1")
+
+    # check that psf is normalized
+    bin_solid_angle = np.diff(cone_solid_angle(fov_bins))
+    e_width = np.diff(energy_bins)
+    assert np.allclose(np.sum((bg.T * e_width).T * bin_solid_angle, axis=1), [1000, 100] / u.s)