diff --git a/pyirf/io/gadf.py b/pyirf/io/gadf.py index 824d00839..963b6de78 100644 --- a/pyirf/io/gadf.py +++ b/pyirf/io/gadf.py @@ -266,7 +266,9 @@ def create_background_2d_hdu( def create_background_3d_hdu( background_3d, reco_energy_bins, - fov_offset_bins, + fov_lon_bins, + fov_lat_bins, + alignment = "ALTAZ", extname="BACKGROUND", **header_cards, ): @@ -282,8 +284,13 @@ def create_background_3d_hdu( (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. + 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 @@ -293,8 +300,8 @@ def create_background_3d_hdu( 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)) + 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) @@ -304,7 +311,7 @@ def create_background_3d_hdu( header["HDUCLAS2"] = "BKG" header["HDUCLAS3"] = "FULL-ENCLOSURE" header["HDUCLAS4"] = "BKG_2D" - header["FOVALIGN"] = "ALTAZ" + header["FOVALIGN"] = alignment header["DATE"] = Time.now().utc.iso _add_header_cards(header, **header_cards) diff --git a/pyirf/irf/background.py b/pyirf/irf/background.py index 1c4242d12..a5f9b3b0d 100644 --- a/pyirf/irf/background.py +++ b/pyirf/irf/background.py @@ -1,7 +1,7 @@ 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") @@ -118,7 +118,9 @@ def background_3d(events, reco_energy_bins, fov_offset_bins, t_obs): 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 - + 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) diff --git a/pyirf/irf/tests/test_background.py b/pyirf/irf/tests/test_background.py index 23f760eb2..5d3e009ae 100644 --- a/pyirf/irf/tests/test_background.py +++ b/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], + )