Add basic facilities for 3D background irfs

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,57 @@
     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_offset_bins,
+    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_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()
+    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_offset_bins[np.newaxis, :].to(u.deg))
+    bkg["DETY_LO"], bkg["DETY_HI"] = binning.split_bin_lo_hi(fov_offset_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"] = "ALTAZ"
+    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

@@ -4,7 +4,7 @@
 from ..utils import cone_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 @@
             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,69 @@
     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 = np.diff(fov_offset_bins)
+    bg_rate = per_energy / t_obs / bin_solid_angle**2
+
+    return bg_rate.to(BACKGROUND_UNIT)