Fixes

docs/changes/276.feature.rst

@@ -0,0 +1 @@
+Add function to make 3d background for lon/lat coordinates.

pyirf/io/gadf.py

@@ -268,7 +268,7 @@ def create_background_3d_hdu(
     reco_energy_bins,
     fov_lon_bins,
     fov_lat_bins,
-    alignment = "ALTAZ",
+    alignment="ALTAZ",
     extname="BACKGROUND",
     **header_cards,
 ):
@@ -289,7 +289,7 @@ def create_background_3d_hdu(
     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
+        Whether 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
@@ -310,8 +310,11 @@ def create_background_3d_hdu(
     header["HDUCLAS1"] = "RESPONSE"
     header["HDUCLAS2"] = "BKG"
     header["HDUCLAS3"] = "FULL-ENCLOSURE"
-    header["HDUCLAS4"] = "BKG_2D"
-    header["FOVALIGN"] = alignment
+    header["HDUCLAS4"] = "BKG_3D"
+    if alignment in ["ALTAZ", "RADEC"]:
+        header["FOVALIGN"] = alignment
+    else:
+        raise ValueError(f"'alignment' must be one of 'ALTAZ' or 'RADEC', got {alignment}")
     header["DATE"] = Time.now().utc.iso
     _add_header_cards(header, **header_cards)
 

pyirf/irf/__init__.py

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

pyirf/irf/background.py

@@ -58,7 +58,7 @@ def background_2d(events, reco_energy_bins, fov_offset_bins, t_obs):
     return bg_rate.to(BACKGROUND_UNIT)
 
 
-def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs):
+def background_3d_lonlat(events, reco_energy_bins, fov_lon_bins, fov_lat_bins, t_obs):
     """
     Calculate background rates in square bins in the field of view.
 
@@ -73,8 +73,12 @@ def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs):
         `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
+    fov_lon_bins: astropy.units.Quantity[angle]
+        The bins in the field of view longitudinal direction to be used for the IRF.
+        Will become DETX bins.
+    fov_lat_bins: astropy.units.Quantity[angle]
+        The bins in the field of view latitudinal direction to be used for the IRF.
+        Will become DETY bins.
     t_obs: astropy.units.Quantity[time]
         Observation time. This must match with how the individual event
         weights are calculated.
@@ -85,24 +89,16 @@ def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs):
         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)
+        Shape: (len(reco_energy_bins) - 1,  len(fov_lon_bins) - 1, len(fov_lon_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}"
-        )
+    fov_x_offset_bins = fov_lon_bins
+    fov_y_offset_bins = fov_lat_bins
 
     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),
+            events["reco_fov_lon"].to_value(u.deg),
+            events["reco_fov_lat"].to_value(u.deg),
         ],
         bins=[
             reco_energy_bins.to_value(u.TeV),
@@ -115,12 +111,14 @@ def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs):
     # 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
+    per_energy = (hist / bin_width_energy[:, np.newaxis, np.newaxis])
 
     # 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:])
+    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

@@ -41,8 +41,8 @@ def test_background():
     )
 
 
-def test_background_3d():
-    from pyirf.irf import background_3d
+def test_background_3d_lonlat():
+    from pyirf.irf import background_3d_lonlat
     from pyirf.utils import rectangle_solid_angle
     from pyirf.irf.background import BACKGROUND_UNIT
 
@@ -103,10 +103,11 @@ def test_background_3d():
         }
     )
 
-    bkg_rate = background_3d(
+    bkg_rate = background_3d_lonlat(
         selected_events,
         reco_energy_bins=reco_energy_bins,
-        fov_offset_bins=np.vstack((fov_lat_bins, fov_lon_bins)),
+        fov_lon_bins=fov_lon_bins,
+        fov_lat_bins=fov_lat_bins,
         t_obs=t_obs,
     )
     assert bkg_rate.shape == (
@@ -118,11 +119,11 @@ def test_background_3d():
 
     # 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,
+        (bin_solid_angle * bkg_rate * bin_width_energy[:, np.newaxis, np.newaxis]).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,
+        (bin_solid_angle * bkg_rate * bin_width_energy[:, np.newaxis, np.newaxis]).sum(axis=(0, 1)) * t_obs,
         2 * [N_tot // 2],
     )