-
Notifications
You must be signed in to change notification settings - Fork 25
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add basic facilities for 3D background irfs #276
base: main
Are you sure you want to change the base?
Changes from all commits
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -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", | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. needs to define the alignment |
||
**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" | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. This should be an argument, not hard-coded |
||
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, | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -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) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. |
||
bg_rate = per_energy / t_obs / bin_solid_angle**2 | ||
|
||
return bg_rate.to(BACKGROUND_UNIT) | ||
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
should have
fov_lon_bins
/fov_lat_bins