Skip to content

Commit

Permalink
Fixes requested by Max, added test
Browse files Browse the repository at this point in the history
  • Loading branch information
@Tobychev
Tobychev committed Nov 8, 2024
1 parent 084a2b8 commit 7fc4c74
Show file tree
Hide file tree
Showing 3 changed files with 109 additions and 11 deletions.
19 changes: 13 additions & 6 deletions pyirf/io/gadf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
):
Expand All @@ -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
Expand All @@ -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)

Expand All @@ -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)

Expand Down
10 changes: 6 additions & 4 deletions pyirf/irf/background.py
View file
Original file line number Diff line number Diff line change
@@ -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")
Expand Down Expand Up @@ -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)
91 changes: 90 additions & 1 deletion pyirf/irf/tests/test_background.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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],
)

0 comments on commit 7fc4c74

Please sign in to comment.