Skip to content

Commit

Permalink
Add more tests on IRF interpolations on more combinations of IRFs
Browse files Browse the repository at this point in the history
  • Loading branch information
@chaimain
chaimain committed Apr 24, 2023
1 parent a61b658 commit 78223b7
Show file tree
Hide file tree
Showing 2 changed files with 213 additions and 71 deletions.
232 changes: 162 additions & 70 deletions lstchain/high_level/tests/test_interpolate.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,12 +5,13 @@


def test_compare_irfs(
simulated_irf_file, simulated_srcdep_irf_file, simulated_dl2_file
simulated_irf_file, simulated_srcdep_irf_file, simulated_dl2_file,
):
from lstchain.high_level.interpolate import compare_irfs
from lstchain.scripts.tests.test_lstchain_scripts import run_program

# Creating a IRF with a different gammaness cut to check the comparison
# Creating IRFs with different global cuts, en-dep cuts for point-like
# IRFs, to make the comparisons
irf_file_2 = simulated_dl2_file.parent / "diff_cut_irf.fits.gz"
run_program(
"lstchain_create_irf_files",
Expand All @@ -29,14 +30,11 @@ def test_compare_irfs(
"lstchain_create_irf_files",
"--input-gamma-dl2",
simulated_dl2_file,
"--input-proton-dl2",
simulated_dl2_file,
"--input-electron-dl2",
simulated_dl2_file,
"--output-irf-file",
irf_file_3,
"--energy-dependent-gh",
"--energy-dependent-theta",
"--point-like",
"--gh-efficiency=0.7",
"--theta-containment=0.7",
)
Expand All @@ -45,14 +43,11 @@ def test_compare_irfs(
"lstchain_create_irf_files",
"--input-gamma-dl2",
simulated_dl2_file,
"--input-proton-dl2",
simulated_dl2_file,
"--input-electron-dl2",
simulated_dl2_file,
"--output-irf-file",
irf_file_4,
"--energy-dependent-gh",
"--energy-dependent-theta",
"--point-like",
"--gh-efficiency=0.8",
"--theta-containment=0.8",
)
Expand All @@ -78,72 +73,167 @@ def test_load_irf_grid(simulated_irf_file):
assert aeff_list.shape == (1, 23, 8)


def test_interp_irf(simulated_irf_file):
def test_interp_irf(simulated_irf_file, simulated_dl2_file):
from lstchain.high_level.interpolate import interpolate_irf

# Create another IRFs with different zenith and azimuth parameter
irf_file_3 = simulated_irf_file.parent / "irf_interp_0.fits.gz"
irf_file_4 = simulated_irf_file.parent / "irf_interp_1.fits.gz"
irf_file_5 = simulated_irf_file.parent / "irf_interp_final.fits.gz"

# Change the effective area for different angular pointings
aeff_1 = Table.read(simulated_irf_file, hdu="EFFECTIVE AREA")
aeff_1_meta = fits.open(simulated_irf_file)["EFFECTIVE AREA"].header
aeff_2 = aeff_1.copy()
aeff_2_meta = aeff_1_meta.copy()

zen_1 = u.Quantity(aeff_1_meta["ZEN_PNT"], "deg").to_value(u.rad)
del_1 = u.Quantity(aeff_1_meta["B_DELTA"], "deg").to_value(u.rad)

zen_2 = 2 * zen_1
del_2 = 1.2 * del_1

factor_zd = (np.cos(zen_2)) / np.cos(zen_1)
factor_del = (np.sin(del_2)) / np.sin(del_1)

aeff_1["EFFAREA"][0] *= factor_zd
aeff_2["EFFAREA"][0] *= factor_zd * factor_del

aeff_1_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
aeff_1_meta["B_DELTA"] = (del_1 * 180 / np.pi, "deg")
aeff_2_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
aeff_2_meta["B_DELTA"] = (del_2 * 180 / np.pi, "deg")

aeff_hdu_1 = fits.BinTableHDU(aeff_1, header=aeff_1_meta, name="EFFECTIVE AREA")
aeff_hdu_2 = fits.BinTableHDU(aeff_2, header=aeff_2_meta, name="EFFECTIVE AREA")

# Change the energy migration for different angular pointings
edisp_1 = Table.read(simulated_irf_file, hdu="ENERGY DISPERSION")
edisp_2 = edisp_1.copy()

edisp_1["MATRIX"][0] *= factor_zd
edisp_2["MATRIX"][0] *= factor_zd * factor_del

edisp_hdu_1 = fits.BinTableHDU(
edisp_1, header=aeff_1_meta, name="ENERGY DISPERSION"
)
edisp_hdu_2 = fits.BinTableHDU(
edisp_2, header=aeff_2_meta, name="ENERGY DISPERSION"
)
# Global cuts IRF
irf_file_g_2 = simulated_irf_file.parent / "irf_interp_0.fits.gz"
irf_file_g_3 = simulated_irf_file.parent / "irf_interp_1.fits.gz"
irf_file_g_final = simulated_irf_file.parent / "irf_interp_final.fits.gz"

hdus_2_g = [
fits.PrimaryHDU(),
]
hdus_3_g = [
fits.PrimaryHDU(),
]

# En-dep, point-like cuts IRF
irf_file_en = simulated_dl2_file.parent / "en_dep_cut_irf_1.fits.gz"

irf_file_en_2 = simulated_dl2_file.parent / "en_dep_cut_irf_2.fits.gz"
irf_file_en_3 = simulated_dl2_file.parent / "en_dep_cut_irf_3.fits.gz"
irf_file_en_final = simulated_dl2_file.parent / "interp_en_dep_cut_irf.fits.gz"

hdus_2_en = [
fits.PrimaryHDU(),
]
hdus_3_en = [
fits.PrimaryHDU(),
]

for irf in [simulated_irf_file, irf_file_en]:
# Change the effective area for different angular pointings
aeff_1 = Table.read(irf, hdu="EFFECTIVE AREA")
aeff_1_meta = fits.open(irf)["EFFECTIVE AREA"].header
aeff_2 = aeff_1.copy()
aeff_2_meta = aeff_1_meta.copy()

zen_1 = u.Quantity(aeff_1_meta["ZEN_PNT"], "deg").to_value(u.rad)
del_1 = u.Quantity(aeff_1_meta["B_DELTA"], "deg").to_value(u.rad)

zen_2 = 2 * zen_1
del_2 = 1.2 * del_1

factor_zd = (np.cos(zen_2)) / np.cos(zen_1)
factor_del = (np.sin(del_2)) / np.sin(del_1)

aeff_1["EFFAREA"][0] *= factor_zd
aeff_2["EFFAREA"][0] *= factor_zd * factor_del

aeff_1_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
aeff_1_meta["B_DELTA"] = (del_1 * 180 / np.pi, "deg")
aeff_2_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
aeff_2_meta["B_DELTA"] = (del_2 * 180 / np.pi, "deg")

aeff_hdu_1 = fits.BinTableHDU(
aeff_1, header=aeff_1_meta, name="EFFECTIVE AREA"
)
aeff_hdu_2 = fits.BinTableHDU(
aeff_2, header=aeff_2_meta, name="EFFECTIVE AREA"
)

# Change the energy migration for different angular pointings
edisp_1 = Table.read(irf, hdu="ENERGY DISPERSION")
edisp_1_meta = fits.open(irf)["ENERGY DISPERSION"].header
edisp_2 = edisp_1.copy()
edisp_2_meta = edisp_1_meta.copy()

edisp_1["MATRIX"][0] *= factor_zd
edisp_2["MATRIX"][0] *= factor_zd * factor_del

edisp_1_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
edisp_1_meta["B_DELTA"] = (del_1 * 180 / np.pi, "deg")
edisp_2_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
edisp_2_meta["B_DELTA"] = (del_2 * 180 / np.pi, "deg")

edisp_hdu_1 = fits.BinTableHDU(
edisp_1, header=edisp_1_meta, name="ENERGY DISPERSION"
)
edisp_hdu_2 = fits.BinTableHDU(
edisp_2, header=edisp_2_meta, name="ENERGY DISPERSION"
)

if "en_dep" in irf.name:
# For GH CUTS, apply the factors for cuts, lower than the max value
gh_1 = Table.read(irf, hdu="GH_CUTS")
gh_1_meta = fits.open(irf)["GH_CUTS"].header
gh_2 = gh_1.copy()
gh_2_meta = gh_1_meta.copy()

mask = gh_1["cut"] < gh_1["cut"].max()

gh_1["cut"][mask] *= factor_zd
gh_2["cut"][mask] *= factor_zd * factor_del

gh_1_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
gh_1_meta["B_DELTA"] = (del_1 * 180 / np.pi, "deg")
gh_2_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
gh_2_meta["B_DELTA"] = (del_2 * 180 / np.pi, "deg")

gh_hdu_1 = fits.BinTableHDU(gh_1, header=gh_1_meta, name="GH_CUTS")
gh_hdu_2 = fits.BinTableHDU(gh_2, header=gh_2_meta, name="GH_CUTS")

# For RAD_MAX CUTS, apply the factors for cuts, lower than the max value
th_1 = Table.read(irf, hdu="RAD_MAX")
th_1_meta = fits.open(irf)["RAD_MAX"].header
th_2 = th_1.copy()
th_2_meta = th_1_meta.copy()

mask = th_1["RAD_MAX"] < th_1["RAD_MAX"].max()

th_1["RAD_MAX"][mask] *= factor_zd
th_2["RAD_MAX"][mask] *= factor_zd * factor_del

th_1_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
th_1_meta["B_DELTA"] = (del_1 * 180 / np.pi, "deg")
th_2_meta["ZEN_PNT"] = (zen_2 * 180 / np.pi, "deg")
th_2_meta["B_DELTA"] = (del_2 * 180 / np.pi, "deg")

th_hdu_1 = fits.BinTableHDU(th_1, header=th_1_meta, name="RAD_MAX")
th_hdu_2 = fits.BinTableHDU(th_2, header=th_2_meta, name="RAD_MAX")

hdus_2_en.append(aeff_hdu_1)
hdus_2_en.append(edisp_hdu_1)
hdus_2_en.append(gh_hdu_1)
hdus_2_en.append(th_hdu_1)

hdus_3_en.append(aeff_hdu_2)
hdus_3_en.append(edisp_hdu_2)
hdus_3_en.append(gh_hdu_2)
hdus_3_en.append(th_hdu_2)
else:
hdus_2_g.append(aeff_hdu_1)
hdus_2_g.append(edisp_hdu_1)

hdus_3_g.append(aeff_hdu_2)
hdus_3_g.append(edisp_hdu_2)

fits.HDUList(hdus_2_g).writeto(irf_file_g_2, overwrite=True)
fits.HDUList(hdus_3_g).writeto(irf_file_g_3, overwrite=True)
fits.HDUList(hdus_2_en).writeto(irf_file_en_2, overwrite=True)
fits.HDUList(hdus_3_en).writeto(irf_file_en_3, overwrite=True)

irfs_g = [simulated_irf_file, irf_file_g_2, irf_file_g_3]
irfs_en = [irf_file_en, irf_file_en_2, irf_file_en_3]
data_pars = {"ZEN_PNT": 30 * u.deg, "B_DELTA": (del_1 * 0.8 * u.rad).to(u.deg)}

fits.HDUList([fits.PrimaryHDU(), aeff_hdu_1, edisp_hdu_1]).writeto(
irf_file_3, overwrite=True
)
hdu_g = interpolate_irf(irfs_g, data_pars)
hdu_g.writeto(irf_file_g_final, overwrite=True)

fits.HDUList([fits.PrimaryHDU(), aeff_hdu_2, edisp_hdu_2]).writeto(
irf_file_4, overwrite=True
)
hdu_en = interpolate_irf(irfs_en, data_pars)
hdu_en.writeto(irf_file_en_final, overwrite=True)

irfs = [simulated_irf_file, irf_file_3, irf_file_4]
data_pars = {"ZEN_PNT": 30 * u.deg, "B_DELTA": (del_1 * 0.8 * u.rad).to(u.deg)}

hdu = interpolate_irf(irfs, data_pars)
hdu.writeto(irf_file_5, overwrite=True)
assert hdu_g[1].header["ZEN_PNT"] == 30
assert irf_file_g_2.exists()
assert irf_file_g_3.exists()
assert irf_file_g_final.exists()

assert hdu[1].header["ZEN_PNT"] == 30
assert irf_file_3.exists()
assert irf_file_4.exists()
assert irf_file_5.exists()
assert hdu_en[1].header["ZEN_PNT"] == 30
assert irf_file_en_2.exists()
assert irf_file_en_3.exists()
assert irf_file_en_final.exists()


def test_check_delaunay_triangles(simulated_irf_file):
Expand All @@ -161,7 +251,9 @@ def test_check_delaunay_triangles(simulated_irf_file):

new_irfs = check_in_delaunay_triangle(irfs, data_pars)
new_irfs2 = check_in_delaunay_triangle(irfs, data_pars2)
new_irfs3 = check_in_delaunay_triangle(irfs, data_pars, use_nearest_irf_node=True)
new_irfs3 = check_in_delaunay_triangle(
irfs, data_pars, use_nearest_irf_node=True
)
new_irfs4 = check_in_delaunay_triangle([], data_pars)

t3 = Table.read(new_irfs3[0], hdu=1).meta
Expand Down
52 changes: 51 additions & 1 deletion lstchain/tools/tests/test_tools.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -253,6 +253,55 @@ def test_create_dl3(temp_dir_observed_files, observed_dl2_file, simulated_irf_fi
== 0
)

@pytest.mark.private_data
def test_create_dl3_irf_interp(
temp_dir_observed_files, observed_dl2_file, simulated_dl2_file
):
"""
Generating an DL3 file from a test DL2 files and using IRF interpolation
function, to either get the interpolated IRF, or the nearest IRF node
"""
from lstchain.tools.lstchain_create_dl3_file import DataReductionFITSWriter

assert (
run_tool(
DataReductionFITSWriter(),
argv=[
f"--input-dl2={observed_dl2_file}",
f"--output-dl3-path={temp_dir_observed_files}",
f"--input-irf-path={simulated_dl2_file.parent}",
"--irf-file-pattern=en_dep_cut*gz",
"--final-irf-file=final_irf_interp_1.fits.gz",
"--source-name=Crab",
"--source-ra=83.633deg",
"--source-dec=22.01deg",
"--overwrite",
],
cwd=temp_dir_observed_files,
)
== 0
)

assert (
run_tool(
DataReductionFITSWriter(),
argv=[
f"--input-dl2={observed_dl2_file}",
f"--output-dl3-path={temp_dir_observed_files}",
f"--input-irf-path={simulated_dl2_file.parent}",
"--irf-file-pattern=en_dep_cut*gz",
"--final-irf-file=final_irf_interp_2.fits.gz",
"--source-name=Crab",
"--source-ra=83.633deg",
"--source-dec=22.01deg",
"--overwrite",
"--use-nearest-irf-node",
],
cwd=temp_dir_observed_files,
)
== 0
)


@pytest.mark.private_data
def test_create_dl3_with_config(temp_dir_observed_files, observed_dl2_file):
Expand Down Expand Up @@ -286,7 +335,8 @@ def test_create_dl3_with_config(temp_dir_observed_files, observed_dl2_file):

@pytest.mark.private_data
def test_create_srcdep_dl3(
temp_dir_observed_srcdep_files, observed_srcdep_dl2_file, simulated_srcdep_irf_file
temp_dir_observed_srcdep_files, observed_srcdep_dl2_file,
simulated_srcdep_irf_file
):
"""
Generating a source-dependent DL3 file from a test DL2 files and test IRF file
Expand Down

0 comments on commit 78223b7

Please sign in to comment.