Skip to content

Commit

Permalink
Fix the shape of the interpolated Effective Area while creating the HDU
Browse files Browse the repository at this point in the history
  • Loading branch information
chaimain committed Jun 25, 2024
1 parent 609a977 commit b9220bd
Show file tree
Hide file tree
Showing 2 changed files with 19 additions and 15 deletions.
24 changes: 12 additions & 12 deletions lstchain/high_level/interpolate.py
Original file line number Diff line number Diff line change
Expand Up @@ -422,7 +422,7 @@ def interpolate_irf(irfs, data_pars, interp_method="linear"):
main_headers = fits.open(irfs[0])[1].header

point_like = main_headers["HDUCLAS3"] == "POINT-LIKE"

# Update headers to be added to the final IRFs
extra_headers = dict()

Expand Down Expand Up @@ -471,19 +471,19 @@ def interpolate_irf(irfs, data_pars, interp_method="linear"):
}

for hdu_name in hdu_names:

if hdu_name == "BACKGROUND":
log.warning("The interpolation of BACKGROUND is not yet supported")
continue

if interp_col_keys[hdu_name] == 'cut':
gadf_irf = False
else:
gadf_irf = True
gadf_irf = True

irf_list = load_irf_grid(
irfs,
extname=hdu_name,
irfs,
extname=hdu_name,
interp_col=interp_col_keys[hdu_name],
gadf_irf=gadf_irf
)
Expand All @@ -492,17 +492,17 @@ def interpolate_irf(irfs, data_pars, interp_method="linear"):
if hdu_name in ['EFFECTIVE AREA', 'ENERGY DISPERSION', 'PSF']:
e_true = np.append(temp_irf["ENERG_LO"][0], temp_irf["ENERG_HI"][0][-1])
fov_off = np.append(temp_irf["THETA_LO"][0], temp_irf["THETA_HI"][0][-1])

if hdu_name == "EFFECTIVE AREA":
aeff_estimator = EffectiveAreaEstimator(
grid_points=irf_pars_sel,
effective_area=irf_list,
interpolator_kwargs={"method": interp_method},
)
aeff_interp = aeff_estimator(interp_pars_sel)

aeff_hdu_interp = create_aeff2d_hdu(
effective_area=aeff_interp.T[0],
effective_area=aeff_interp[0],
true_energy_bins=e_true,
fov_offset_bins=fov_off,
point_like=point_like,
Expand All @@ -519,7 +519,7 @@ def interpolate_irf(irfs, data_pars, interp_method="linear"):
energy_dispersion=irf_list,
)
edisp_interp = edisp_estimator(interp_pars_sel)

edisp_hdu_interp = create_energy_dispersion_hdu(
energy_dispersion=edisp_interp[0],
true_energy_bins=e_true,
Expand Down Expand Up @@ -556,7 +556,7 @@ def interpolate_irf(irfs, data_pars, interp_method="linear"):
cut_header = fits.Header()
cut_header["CREATOR"] = f"lstchain v{__version__}"
cut_header["DATE"] = Time.now().utc.iso

for k, v in extra_headers.items():
cut_header[k] = v

Expand Down
10 changes: 7 additions & 3 deletions lstchain/high_level/tests/test_interpolate.py
Original file line number Diff line number Diff line change
Expand Up @@ -257,17 +257,17 @@ def test_interp_irf(
psf_1_meta = fits.open(irf)["PSF"].header
psf_2 = psf_1.copy()
psf_2_meta = psf_1_meta.copy()

psf_1["RPSF"][0] *= factor_czd
psf_2["RPSF"][0] *= factor_czd * factor_sd

psf_1_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
psf_1_meta["AZ_PNT"] = (az_1 * 180 / np.pi, "deg")
psf_1_meta["B_DELTA"] = (del_1 * 180 / np.pi, "deg")
psf_2_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
psf_2_meta["AZ_PNT"] = (az_2 * 180 / np.pi, "deg")
psf_2_meta["B_DELTA"] = (del_2 * 180 / np.pi, "deg")

psf_hdu_1 = fits.BinTableHDU(
psf_1, header=psf_1_meta, name="PSF"
)
Expand Down Expand Up @@ -333,10 +333,14 @@ def test_interp_irf(
hdu_en_srcdep = interpolate_irf(irfs_en_srcdep, data_pars)
hdu_en_srcdep.writeto(irf_file_en_srcdep_final, overwrite=True)

aeff_shape_final = Table.read(irf_file_g_final, hdu=1)["EFFAREA"].shape
aeff_shape_2 = Table.read(irf_file_g_2, hdu=1)["EFFAREA"].shape

assert hdu_g[1].header["ZEN_PNT"] == zen_t
assert irf_file_g_2.exists()
assert irf_file_g_3.exists()
assert irf_file_g_final.exists()
assert aeff_shape_final == aeff_shape_2

assert hdu_en[1].header["ZEN_PNT"] == zen_t
assert irf_file_en_2.exists()
Expand Down

0 comments on commit b9220bd

Please sign in to comment.