Skip to content
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

Allow for scaling true energy in IRFs to assess energy-scale systematics #1261

Merged
merged 35 commits into from
Jul 8, 2024
Merged
Show file tree
Hide file tree
Changes from 25 commits
Commits
Show all changes
35 commits
Select commit Hold shift + click to select a range
b85e538
Update event_selection.py
deborahDOR May 27, 2024
14d53d1
Update lstchain_create_irf_files.py
deborahDOR May 27, 2024
ee2ff27
Update irf_dl3_tool_config.json
deborahDOR May 27, 2024
0ac1de3
Update event_selection.py
deborahDOR May 27, 2024
e0b8f9a
Update lstchain_create_irf_files.py
deborahDOR May 27, 2024
00dbe27
Update irf_dl3_tool_config.json
deborahDOR Jun 7, 2024
de4b751
Update lstchain_create_irf_files.py
deborahDOR Jun 7, 2024
e06e0ac
Update lstchain_create_irf_files.py
deborahDOR Jun 7, 2024
eb1c4ab
Update lstchain_create_irf_files.py
deborahDOR Jun 11, 2024
eac6f10
Update lstchain_create_irf_files.py
deborahDOR Jun 12, 2024
319258a
Update lstchain_create_irf_files.py
deborahDOR Jun 12, 2024
f28019d
Update test_tools.py
deborahDOR Jul 2, 2024
d9673e1
Update test_tools.py
deborahDOR Jul 2, 2024
c239f72
Create irf_dl3_tool_config_mod.json
deborahDOR Jul 2, 2024
3d0a153
Update test_tools.py
deborahDOR Jul 3, 2024
13c085a
Merge branch 'main' into energy_systematics_in_irfs
chaimain Jul 3, 2024
b8fb8c5
Update test_tools.py
deborahDOR Jul 3, 2024
7a92409
Delete docs/examples/irf_dl3_tool_config_mod.json
deborahDOR Jul 3, 2024
be83a58
Update lstchain_create_irf_files.py
deborahDOR Jul 3, 2024
f8614a2
Update test_tools.py
deborahDOR Jul 4, 2024
a41fa4f
Update lstchain_create_irf_files.py
deborahDOR Jul 4, 2024
dc361ce
Update lstchain_create_irf_files.py
deborahDOR Jul 4, 2024
9a0605a
Update test_tools.py
deborahDOR Jul 4, 2024
221aaf7
Update test_tools.py
deborahDOR Jul 4, 2024
264f39c
[skip ci] typo
morcuended Jul 4, 2024
ce611f5
Improve the exploration of Simtel config history to handle more LST M…
gabemery Jul 3, 2024
8b7d8f6
Fix test
gabemery Jul 3, 2024
b28556c
Convert to float
gabemery Jul 3, 2024
edfaf39
Remove try block without handled exception.
gabemery Jul 4, 2024
b0820e6
Add warning log message
gabemery Jul 4, 2024
c137ade
Move comment about NSB extraction from simtel file to the function do…
gabemery Jul 5, 2024
037f900
Run black and add the missing pytest marker
chaimain Jul 8, 2024
66a99d7
Update test_tools.py by running black and adding a missing pytest marker
chaimain Jul 8, 2024
980aabe
Update test_tools.py
chaimain Jul 8, 2024
e999002
Run pyflakes to clean up rebase issues
chaimain Jul 8, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions docs/examples/irf_dl3_tool_config.json
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@
"true_energy_min": 0.005,
"true_energy_max": 500,
"true_energy_n_bins": 25,
"scale_true_energy": 1.0,
"reco_energy_min": 0.005,
"reco_energy_max": 500,
"reco_energy_n_bins": 25,
Expand Down
5 changes: 5 additions & 0 deletions lstchain/io/event_selection.py
Original file line number Diff line number Diff line change
Expand Up @@ -364,6 +364,11 @@ class DataBinning(Component):
default_value=25,
).tag(config=True)

scale_true_energy= Float(
help="Scaling value for True energy",
default_value=1.0,
).tag(config=True)

reco_energy_min = Float(
help="Minimum value for Reco Energy bins in TeV units",
default_value=0.005,
Expand Down
28 changes: 26 additions & 2 deletions lstchain/tools/lstchain_create_irf_files.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,15 @@
If you want to generate source-dependent IRFs, source-dep flag should be activated.
The global alpha cut used to generate IRFs is stored as AL_CUT in the HDU header.

Modified IRFs with true energy scaled by a given factor can be created to evaluate
the systematic uncertainty in the light collection efficiency. This can be done by
setting a value different from one for the "scale_true_energy" argument present in
the DataBinning Component of the configuration file of the IRF creation Tool.
(The true energy of the MC events will be scaled before filling the IRFs histograms
when pyirf commands are used. The effects expected are a non-diagonal energy dispersion
matrix and a different spectrum).


"""

from astropy import table
Expand Down Expand Up @@ -136,7 +145,12 @@ class IRFFITSWriter(Tool):
--global-alpha-cut 10
--source-dep

"""
To build modified IRFs by specifying a scaling factor applying to the true energy (without using a config file):
> lstchain_create_irf_files
-g /path/to/DL2_MC_gamma_file.h5
-o /path/to/irf.fits.gz
--scale-true-energy 1.15
"""

input_gamma_dl2 = traits.Path(
help="Input MC gamma DL2 file",
Expand Down Expand Up @@ -221,6 +235,7 @@ class IRFFITSWriter(Tool):
"global-alpha-cut": "DL3Cuts.global_alpha_cut",
"allowed-tels": "DL3Cuts.allowed_tels",
"overwrite": "IRFFITSWriter.overwrite",
"scale-true-energy": "DataBinning.scale_true_energy"
}

flags = {
Expand Down Expand Up @@ -317,6 +332,12 @@ def start(self):
p["geomag_params"],
) = read_mc_dl2_to_QTable(p["file"])


if self.data_bin.scale_true_energy != 1.0:
p["events"]["true_energy"] *= self.data_bin.scale_true_energy
p["simulation_info"].energy_min *= self.data_bin.scale_true_energy
p["simulation_info"].energy_max *= self.data_bin.scale_true_energy

p["mc_type"] = check_mc_type(p["file"])

self.log.debug(
Expand Down Expand Up @@ -527,7 +548,10 @@ def start(self):
geomag_params["GEOMAG_DELTA"].to_value(u.deg),
"deg",
)

extra_headers["ETRUE_SCALE"]= (
self.data_bin.scale_true_energy
)

if self.point_like:
self.log.info("Generating point_like IRF HDUs")
else:
Expand Down
101 changes: 100 additions & 1 deletion lstchain/tools/tests/test_tools.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@
from astropy.io import fits
import numpy as np


def test_create_irf_full_enclosure(temp_dir_observed_files, simulated_dl2_file):
"""
Generating full enclosure IRF file from a test DL2 files
Expand Down Expand Up @@ -415,3 +414,103 @@ def test_index_srcdep_dl3_files(temp_dir_observed_srcdep_files):
'EFFECTIVE AREA', 'ENERGY DISPERSION'
]:
assert hdu_name in data.hdu_table['HDU_NAME']

def test_add_scale_true_energy_in_irfs(temp_dir_observed_files, simulated_dl2_file):
"""
Checking the validity of modified IRFs after scaling the True Energy by a factor.
"""

from lstchain.tools.lstchain_create_irf_files import IRFFITSWriter
from gammapy.irf import EffectiveAreaTable2D, EnergyDispersion2D
import astropy.units as u

irf_file=temp_dir_observed_files / "fe_irf.fits.gz"
irf_file_mod = temp_dir_observed_files / "mod_irf.fits.gz"
config_file = os.path.join(os.getcwd(), "docs/examples/irf_dl3_tool_config.json")

assert (
run_tool(
IRFFITSWriter(),
argv=[
f"--input-gamma-dl2={simulated_dl2_file}",
f"--input-proton-dl2={simulated_dl2_file}",
f"--input-electron-dl2={simulated_dl2_file}",
f"--output-irf-file={irf_file}",
f"--config={config_file}",
"--overwrite",
"--DataBinning.true_energy_n_bins=2",
"--DataBinning.reco_energy_n_bins=2",
"--DataBinning.true_energy_min: 0.2",
"--DataBinning.true_energy_max: 0.3",
"--DL3Cuts.min_event_p_en_bin=2",
],
cwd=temp_dir_observed_files,
)
== 0
)
assert (
run_tool(
IRFFITSWriter(),
argv=[
f"--input-gamma-dl2={simulated_dl2_file}",
f"--input-proton-dl2={simulated_dl2_file}",
f"--input-electron-dl2={simulated_dl2_file}",
f"--output-irf-file={irf_file_mod}",
f"--config={config_file}",
"--overwrite",
"--DataBinning.true_energy_n_bins=2",
"--DataBinning.reco_energy_n_bins=2",
"--DataBinning.true_energy_min: 0.2",
"--DataBinning.true_energy_max: 0.3",
"--DL3Cuts.min_event_p_en_bin=2",
"--DataBinning.scale_true_energy=1.15",
],
cwd=temp_dir_observed_files,
)
== 0
)

aeff_hdu=EffectiveAreaTable2D.read(irf_file, hdu='EFFECTIVE AREA')
aeff_mod_hdu=EffectiveAreaTable2D.read(irf_file_mod, hdu='EFFECTIVE AREA')

edisp_hdu=EnergyDispersion2D.read(irf_file, hdu='ENERGY DISPERSION')
edisp_mod_hdu=EnergyDispersion2D.read(irf_file_mod, hdu='ENERGY DISPERSION')

assert aeff_mod_hdu.data.shape==aeff_hdu.data.shape
assert edisp_mod_hdu.data.shape==edisp_hdu.data.shape

edisp=EnergyDispersion2D.read(irf_file)
edisp_mod=EnergyDispersion2D.read(irf_file_mod)

e_migra = edisp.axes["migra"].center
e_migra_mod=edisp_mod.axes["migra"].center

e_true_list=[0.2,2,20]
e_migra_prob=[]
e_migra_prob_mod=[]

for i in e_true_list:
e_true = i* u.TeV
e_migra_prob.append(edisp.evaluate(
offset=0.4*u.deg,
energy_true=e_true,
migra=e_migra
))
e_migra_prob_mod.append(edisp_mod.evaluate(
offset=0.4*u.deg,
energy_true=e_true,
migra=e_migra_mod
))

#check that the maximum of the density probability of the migration has shifted
order_max = []
order_max_mod = []
for idx, _ in enumerate(e_true_list):
for j in range(len(e_migra)):
if e_migra_prob[idx][j] > e_migra_prob[idx][j-1]:
order_max.append(j)
if e_migra_prob_mod[idx][j] > e_migra_prob_mod[idx][j-1]:
order_max_mod.append(j)

for i in range(len(order_max)):
assert order_max[i]!=order_max_mod[i]