diff --git a/docs/examples/irf_dl3_tool_config.json b/docs/examples/irf_dl3_tool_config.json index 61de3ef1d6..a9c75a9e39 100644 --- a/docs/examples/irf_dl3_tool_config.json +++ b/docs/examples/irf_dl3_tool_config.json @@ -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, diff --git a/lstchain/io/event_selection.py b/lstchain/io/event_selection.py index 533e45e1e2..c13f87c69c 100644 --- a/lstchain/io/event_selection.py +++ b/lstchain/io/event_selection.py @@ -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, diff --git a/lstchain/io/io.py b/lstchain/io/io.py index 04cd0c9ffe..79a82ad556 100644 --- a/lstchain/io/io.py +++ b/lstchain/io/io.py @@ -1,6 +1,5 @@ import logging import os -import re import warnings from multiprocessing import Pool from contextlib import ExitStack @@ -22,9 +21,8 @@ from ctapipe.instrument import SubarrayDescription from ctapipe.io import HDF5TableReader, HDF5TableWriter -from eventio import Histograms, EventIOFile -from eventio.search_utils import yield_toplevel_of_type, yield_all_subobjects -from eventio.simtel.objects import History, HistoryConfig +from eventio import Histograms, SimTelFile +from eventio.search_utils import yield_toplevel_of_type from pyirf.simulations import SimulatedEventsInfo @@ -59,7 +57,6 @@ 'global_metadata', 'merge_dl2_runs', 'merging_check', - 'parse_cfg_bytestring', 'read_data_dl2_to_QTable', 'read_dl2_params', 'read_mc_dl2_to_QTable', @@ -1264,37 +1261,34 @@ def remove_duplicated_events(data): data.remove_rows(remove_row_list) -def parse_cfg_bytestring(bytestring): - """ - Parse configuration as read by eventio - :param bytes bytestring: A ``Bytes`` object with configuration data for one parameter - :return: Tuple in form ``('parameter_name', 'value')`` +def extract_simulation_nsb(filename): """ - line_decoded = bytestring.decode('utf-8').rstrip() - if 'ECHO' in line_decoded or '#' in line_decoded: - return None - line_list = line_decoded.split('%', 1)[0] # drop comment - res = re.sub(' +', ' ', line_list).strip().split(' ', 1) # remove extra whitespaces and split - return res[0].upper(), res[1] + Get current run NSB from configuration in simtel file. + WARNING : In current MC, correct NSB are logged after 'STORE_PHOTOELECTRONS' entries + In any new production, behaviour needs to be verified. + New version of simtel will allow to use better metadata. -def extract_simulation_nsb(filename): - """ - Get current run NSB from configuration in simtel file :param str filename: Input file name - :return array of `float` by tel_id: NSB rate - """ - nsb = [] - with EventIOFile(filename) as f: - for o in yield_all_subobjects(f, [History, HistoryConfig]): - if hasattr(o, 'parse'): - try: - cfg_element = parse_cfg_bytestring(o.parse()[1]) - if cfg_element is not None: - if cfg_element[0] == 'NIGHTSKY_BACKGROUND': - nsb.append(float(cfg_element[1].strip('all:'))) - except Exception as e: - print('Unexpected end of %s,\n caught exception %s', filename, e) + :return dict of `float` by tel_id: NSB rate + """ + nsb = {} + next_nsb = False + tel_id = 1 + with SimTelFile(filename) as f: + for _, line in f.history: + line = line.decode('utf-8').strip().split(' ') + if next_nsb and line[0] == 'NIGHTSKY_BACKGROUND': + nsb[tel_id] = float(line[1].strip('all:')) + tel_id = tel_id+1 + if line[0] == 'STORE_PHOTOELECTRONS': + next_nsb = True + else: + next_nsb = False + log.warning('Original MC night sky background extracted from the config history in the simtel file.\n' + 'This is done for existing LST MC such as the one created using: ' + 'https://github.com/cta-observatory/lst-sim-config/tree/sim-tel_LSTProd2_MAGICST0316' + '\nExtracted values are: ' + str(np.asarray(nsb)) + 'GHz. Check that it corresponds to expectations.') return nsb diff --git a/lstchain/io/tests/test_io.py b/lstchain/io/tests/test_io.py index 8f65bd674f..6442eff66e 100644 --- a/lstchain/io/tests/test_io.py +++ b/lstchain/io/tests/test_io.py @@ -156,8 +156,7 @@ def test_extract_simulation_nsb(mc_gamma_testfile): from lstchain.io.io import extract_simulation_nsb nsb = extract_simulation_nsb(mc_gamma_testfile) - assert np.isclose(nsb[0], 0.246, rtol=0.1) - assert np.isclose(nsb[1], 0.217, rtol=0.1) + assert np.isclose(nsb[1], 0.246, rtol=0.1) def test_remove_duplicated_events(): diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py index 1dbe3be1aa..3f2db65a45 100644 --- a/lstchain/reco/r0_to_dl1.py +++ b/lstchain/reco/r0_to_dl1.py @@ -433,10 +433,8 @@ def r0_to_dl1( charge_spe_cumulative_pdf = interp1d(spe_integral, spe[0], kind='cubic', bounds_error=False, fill_value=0., assume_sorted=True) - allowed_tel = np.zeros(len(nsb_original), dtype=bool) - allowed_tel[np.array(config['source_config']['LSTEventSource']['allowed_tels'])] = True logger.info('Tuning NSB on MC waveform from ' - + str(np.asarray(nsb_original)[allowed_tel]) + + str(np.asarray(nsb_original)) + 'GHz to {0:d}%'.format(int(nsb_tuning_ratio * 100 + 100.5)) + ' for telescopes ids ' + str(config['source_config']['LSTEventSource']['allowed_tels'])) nsb_tuning_args = [nsb_tuning_ratio, nsb_original, pulse_template, charge_spe_cumulative_pdf] diff --git a/lstchain/tools/lstchain_create_irf_files.py b/lstchain/tools/lstchain_create_irf_files.py index 1186570289..d3c2e464b3 100644 --- a/lstchain/tools/lstchain_create_irf_files.py +++ b/lstchain/tools/lstchain_create_irf_files.py @@ -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 @@ -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", @@ -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 = { @@ -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( @@ -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: diff --git a/lstchain/tools/tests/test_tools.py b/lstchain/tools/tests/test_tools.py index c02a4e202b..3e5b96b900 100644 --- a/lstchain/tools/tests/test_tools.py +++ b/lstchain/tools/tests/test_tools.py @@ -381,11 +381,16 @@ def test_index_dl3_files(temp_dir_observed_files): assert 2008 in data.obs_table["OBS_ID"] for hdu_name in [ - 'EVENTS', 'GTI', 'POINTING', - 'EFFECTIVE AREA', 'ENERGY DISPERSION', - 'BACKGROUND', 'PSF' + "EVENTS", + "GTI", + "POINTING", + "EFFECTIVE AREA", + "ENERGY DISPERSION", + "BACKGROUND", + "PSF", ]: - assert hdu_name in data.hdu_table['HDU_NAME'] + assert hdu_name in data.hdu_table["HDU_NAME"] + @pytest.mark.private_data def test_index_srcdep_dl3_files(temp_dir_observed_srcdep_files): @@ -411,7 +416,116 @@ def test_index_srcdep_dl3_files(temp_dir_observed_srcdep_files): assert 2008 in data.obs_table["OBS_ID"] for hdu_name in [ - 'EVENTS', 'GTI', 'POINTING', - 'EFFECTIVE AREA', 'ENERGY DISPERSION' + "EVENTS", + "GTI", + "POINTING", + "EFFECTIVE AREA", + "ENERGY DISPERSION", ]: - assert hdu_name in data.hdu_table['HDU_NAME'] + assert hdu_name in data.hdu_table["HDU_NAME"] + + +@pytest.mark.private_data +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. + """ + + import astropy.units as u + from gammapy.irf import EffectiveAreaTable2D, EnergyDispersion2D + from lstchain.tools.lstchain_create_irf_files import IRFFITSWriter + + 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]