Add basic facilities for 3D background irfs

README.rst

@@ -29,3 +29,20 @@ please cite it by using the corresponding DOI:
 
 - latest : |doilatest|
 - all versions: `Please visit Zenodo <https://zenodo.org/record/5833284>`_ and select the correct version
+
+At this point, our latest publication is the `2023 ICRC proceeding <https://doi.org/10.22323/1.444.0618>`_, which you can
+cite using the following bibtex entry, especially if using functionalities from ``pyirf.interpolation``:
+
+.. code::
+
+   @inproceedings{pyirf-icrc-2023,
+     author = {Dominik, Rune Michael and Linhoff, Maximilian and Sitarek, Julian},
+     title = {Interpolation of Instrument Response Functions for the Cherenkov Telescope Array in the Context of pyirf},
+     usera = {for the CTA Consortium},
+     doi = {10.22323/1.444.0703},
+     booktitle = {Proceedings, 38th International Cosmic Ray Conference},
+     year=2023,
+     volume={444},
+     number={618},
+     location={Nagoya, Japan},
+   }

pyirf/io/__init__.py

@@ -5,6 +5,7 @@
     create_psf_table_hdu,
     create_rad_max_hdu,
     create_background_2d_hdu,
+    create_background_3d_hdu,
 )
 
 
@@ -16,4 +17,5 @@
     "create_psf_table_hdu",
     "create_rad_max_hdu",
     "create_background_2d_hdu",
+    "create_background_3d_hdu",
 ]

pyirf/io/gadf.py

@@ -260,6 +260,64 @@ def create_background_2d_hdu(
     return BinTableHDU(bkg, header=header, name=extname)
 
 
+@u.quantity_input(
+    background=GADF_BACKGROUND_UNIT, reco_energy_bins=u.TeV, fov_offset_bins=u.deg,
+)
+def create_background_3d_hdu(
+    background_3d,
+    reco_energy_bins,
+    fov_lon_bins,
+    fov_lat_bins,
+    alignment = "ALTAZ",
+    extname="BACKGROUND",
+    **header_cards,
+):
+    """
+    Create a fits binary table HDU in GADF format for the background 3d table. Assumes ALTAZ coordinates.
+    See the specification at
+    https://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/full_enclosure/bkg/index.html#bkg-2d
+
+    Parameters
+    ----------
+    background_3d: astropy.units.Quantity[(MeV s sr)^-1]
+        Background rate, must have shape
+        (n_energy_bins, n_fov_offset_bins, n_fov_offset_bins)
+    reco_energy_bins: astropy.units.Quantity[energy]
+        Bin edges in reconstructed energy
+    fov_lon_bins: astropy.units.Quantity[angle]
+        Bin edges in the field of view system, becomes the DETX values
+    fov_lat_bins: astropy.units.Quantity[angle]
+        Bin edges in the field of view system, becomes the DETY values
+    alignment: str
+        Wheter the FOV coordinates are aligned with the ALTAZ or RADEC system, more details at
+        https://gamma-astro-data-formats.readthedocs.io/en/latest/general/coordinates.html
+    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()
+    bkg["ENERG_LO"], bkg["ENERG_HI"] = binning.split_bin_lo_hi(reco_energy_bins[np.newaxis, :].to(u.TeV))
+    bkg["DETX_LO"], bkg["DETX_HI"] = binning.split_bin_lo_hi(fov_lon_bins[np.newaxis, :].to(u.deg))
+    bkg["DETY_LO"], bkg["DETY_HI"] = binning.split_bin_lo_hi(fov_lat_bins[np.newaxis, :].to(u.deg))
+    # transpose as FITS uses opposite dimension order
+    bkg["BKG"] = background_3d.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["FOVALIGN"] = alignment
+    header["DATE"] = Time.now().utc.iso
+    _add_header_cards(header, **header_cards)
+
+    return BinTableHDU(bkg, header=header, name=extname)
+
+
 @u.quantity_input(
     rad_max=u.deg,
     reco_energy_bins=u.TeV,

pyirf/irf/__init__.py

@@ -5,7 +5,7 @@
 )
 from .energy_dispersion import energy_dispersion
 from .psf import psf_table
-from .background import background_2d
+from .background import background_2d, background_3d
 
 __all__ = [
     "effective_area",
@@ -14,4 +14,5 @@
     "energy_dispersion",
     "psf_table",
     "background_2d",
+    "background_3d",
 ]

pyirf/irf/background.py

@@ -1,10 +1,10 @@
 import astropy.units as u
 import numpy as np
 
-from ..utils import cone_solid_angle
+from ..utils import cone_solid_angle, rectangle_solid_angle
 
 #: Unit of the background rate IRF
-BACKGROUND_UNIT = u.Unit('s-1 TeV-1 sr-1')
+BACKGROUND_UNIT = u.Unit("s-1 TeV-1 sr-1")
 
 
 def background_2d(events, reco_energy_bins, fov_offset_bins, t_obs):
@@ -43,7 +43,7 @@ def background_2d(events, reco_energy_bins, fov_offset_bins, t_obs):
             reco_energy_bins.to_value(u.TeV),
             fov_offset_bins.to_value(u.deg),
         ],
-        weights=events['weight'],
+        weights=events["weight"],
     )
 
     # divide all energy bins by their width
@@ -56,3 +56,71 @@ def background_2d(events, reco_energy_bins, fov_offset_bins, t_obs):
     bg_rate = per_energy / t_obs / bin_solid_angle
 
     return bg_rate.to(BACKGROUND_UNIT)
+
+
+def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs):
+    """
+    Calculate background rates in square 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-3d
+
+    Parameters
+    ----------
+    events: astropy.table.QTable
+        DL2 events table of the selected background events.
+        Needed columns for this function: `reco_fov_lon`, `reco_fov_lat`,
+        `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, either a (N,) or (1,N) or a (2,N) array
+    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, len(fov_offset_bins) - 1)
+    """
+    if (fov_offset_bins.shape[0] == 1) or (len(fov_offset_bins.shape) == 1):
+        fov_x_offset_bins = fov_offset_bins
+        fov_y_offset_bins = fov_offset_bins
+    elif fov_offset_bins.shape[0] == 2:
+        fov_x_offset_bins = fov_offset_bins[0, :]
+        fov_y_offset_bins = fov_offset_bins[1, :]
+    else:
+        raise ValueError(
+            f"fov_offset_bins must be eiher (N,) or (1,N) or (2,N), found {fov_offset_bins.shape}"
+        )
+
+    hist, _ = np.histogramdd(
+        [
+            events["reco_energy"].to_value(u.TeV),
+            (events["reco_fov_lon"]).to_value(u.deg),
+            (events["reco_fov_lat"]).to_value(u.deg),
+        ],
+        bins=[
+            reco_energy_bins.to_value(u.TeV),
+            fov_x_offset_bins.to_value(u.deg),
+            fov_y_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  = rectangle_solid_angle(fov_x_offset_bins[:-1],
+                                             fov_x_offset_bins[1:],
+                                             fov_y_offset_bins[:-1],
+                                             fov_y_offset_bins[1:])
+    bg_rate = per_energy / t_obs / bin_solid_angle
+    return bg_rate.to(BACKGROUND_UNIT)

pyirf/irf/tests/test_background.py

@@ -36,4 +36,93 @@ def test_background():
     # 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)
+    assert np.allclose(
+        np.sum((bg.T * e_width).T * bin_solid_angle, axis=1), [1000, 100] / u.s
+    )
+
+
+def test_background_3d():
+    from pyirf.irf import background_3d
+    from pyirf.utils import rectangle_solid_angle
+    from pyirf.irf.background import BACKGROUND_UNIT
+
+    reco_energy_bins = [0.1, 1.1, 11.1, 111.1] * u.TeV
+    fov_lon_bins = [-1.0, 0, 1.0] * u.deg
+    fov_lat_bins = [-1.0, 0, 1.0] * u.deg
+
+    N_low = 4000
+    N_high = 40
+    N_tot = N_low + N_high
+
+    # Fill values
+    E_low, E_hig = 0.5, 5
+    Lon_low, Lon_hig = (-0.5, 0.5) * u.deg
+    Lat_low, Lat_hig = (-0.5, 0.5) * u.deg
+
+    t_obs = 100 * u.s
+    bin_width_energy = np.diff(reco_energy_bins)
+    bin_solid_angle = rectangle_solid_angle(
+        fov_lon_bins[:-1], fov_lon_bins[1:], fov_lat_bins[:-1], fov_lat_bins[1:]
+    )
+
+    # Toy events with two energies and four different sky positions
+    selected_events = QTable(
+        {
+            "reco_energy": np.concatenate(
+                [
+                    np.full(N_low // 4, E_low),
+                    np.full(N_high // 4, E_hig),
+                    np.full(N_low // 4, E_low),
+                    np.full(N_high // 4, E_hig),
+                    np.full(N_low // 4, E_low),
+                    np.full(N_high // 4, E_hig),
+                    np.full(N_low // 4, E_low),
+                    np.full(N_high // 4, E_hig),
+                ]
+            )
+            * u.TeV,
+            "reco_fov_lon": np.concatenate(
+                [
+                    np.full(N_low // 4, Lon_low),
+                    np.full(N_high // 4, Lon_hig),
+                    np.full(N_low // 4, Lon_low),
+                    np.full(N_high // 4, Lon_hig),
+                    np.full(N_low // 4, Lon_low),
+                    np.full(N_high // 4, Lon_hig),
+                    np.full(N_low // 4, Lon_low),
+                    np.full(N_high // 4, Lon_hig),
+                ]
+            )
+            * u.deg,
+            "reco_fov_lat": np.append(
+                np.full(N_tot // 2, Lat_low),
+                np.full(N_tot // 2, Lat_hig)
+            )
+            * u.deg,
+            "weight": np.full(N_tot, 1.0),
+        }
+    )
+
+    bkg_rate = background_3d(
+        selected_events,
+        reco_energy_bins=reco_energy_bins,
+        fov_offset_bins=np.vstack((fov_lat_bins, fov_lon_bins)),
+        t_obs=t_obs,
+    )
+    assert bkg_rate.shape == (
+        len(reco_energy_bins) - 1,
+        len(fov_lon_bins) - 1,
+        len(fov_lat_bins) - 1,
+    )
+    assert bkg_rate.unit == BACKGROUND_UNIT
+
+    # Convert to counts, project to energy axis, and check counts round-trip correctly
+    assert np.allclose(
+        (bin_solid_angle * (bkg_rate.T * bin_width_energy).T).sum(axis=(1, 2)) * t_obs,
+        [N_low, N_high, 0],
+    )
+    # Convert to counts, project to latitude axis, and check counts round-trip correctly
+    assert np.allclose(
+        (bin_solid_angle * (bkg_rate.T * bin_width_energy).T).sum(axis=(0, 1)) * t_obs,
+        2 * [N_tot // 2],
+    )

setup.py

@@ -16,20 +16,17 @@
         "notebook",
         "tables",
         "towncrier",
-        "astropy~=5.3",  # gammapy doesn't yet support 6.0 but doesn't pin it
         gammapy,
     ],
     "tests": [
         "pytest",
         "pytest-cov",
-        "astropy~=5.3",  # gammapy doesn't yet support 6.0 but doesn't pin it
         gammapy,
         "ogadf-schema~=0.2.3",
         "uproot~=4.0",
         "awkward~=1.0",
     ],
     "gammapy": [
-        "astropy~=5.3",  # gammapy doesn't yet support 6.0 but doesn't pin it
         gammapy,
     ],
 }
@@ -48,7 +45,7 @@
     install_requires=[
         "astropy>=5.3,<7.0.0a0",
         "numpy>=1.21",
-        "scipy<1.12",
+        "scipy",
         "tqdm",
     ],
     include_package_data=True,