From e31408b685a8626f6ce534bada9abb6c0f74ab65 Mon Sep 17 00:00:00 2001 From: SeiyaNozaki Date: Mon, 3 Jul 2023 14:14:01 +0000 Subject: [PATCH 01/12] use effective focal length to compute src-dep parameters --- lstchain/reco/dl1_to_dl2.py | 52 +++++++++++++------ ...stchain_add_source_dependent_parameters.py | 17 ++++-- lstchain/scripts/lstchain_dl1_to_dl2.py | 22 +++++--- 3 files changed, 61 insertions(+), 30 deletions(-) diff --git a/lstchain/reco/dl1_to_dl2.py b/lstchain/reco/dl1_to_dl2.py index 116094898f..cadf9c4c73 100644 --- a/lstchain/reco/dl1_to_dl2.py +++ b/lstchain/reco/dl1_to_dl2.py @@ -31,6 +31,7 @@ from ctapipe.image.hillas import camera_to_shower_coordinates from ctapipe.instrument import SubarrayDescription +from ctapipe_io_lst import OPTICS __all__ = [ 'apply_models', @@ -349,10 +350,18 @@ def build_models(filegammas, fileprotons, src_dep_df_gamma = get_srcdep_params(filegammas) else: - subarray_info = SubarrayDescription.from_hdf(filegammas) - tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 - focal_length = subarray_info.tel[tel_id].optics.equivalent_focal_length - src_dep_df_gamma = get_source_dependent_parameters(df_gamma, config, focal_length=focal_length) + try: + subarray_info = SubarrayDescription.from_hdf(filegammas) + tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 + effective_focal_length = subarray_info.tel[tel_id].optics.effective_focal_length + except OSError: + print("subarray table is not readable because of the version inompatibility.") + print("Use the effective focal lentgh for the standard LST optics") + effective_focal_length = OPTICS.effective_focal_length + + src_dep_df_gamma = get_source_dependent_parameters( + df_gamma, config, effective_focal_length=effective_focal_length + ) df_gamma = pd.concat([df_gamma, src_dep_df_gamma['on']], axis=1) @@ -362,10 +371,18 @@ def build_models(filegammas, fileprotons, src_dep_df_proton = get_srcdep_params(fileprotons) else: - subarray_info = SubarrayDescription.from_hdf(fileprotons) - tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 - focal_length = subarray_info.tel[tel_id].optics.equivalent_focal_length - src_dep_df_proton = get_source_dependent_parameters(df_proton, config, focal_length=focal_length) + try: + subarray_info = SubarrayDescription.from_hdf(fileprotons) + tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 + effective_focal_length = subarray_info.tel[tel_id].optics.effective_focal_length + except OSError: + print("subarray table is not readable because of the version inompatibility.") + print("Use the effective focal lentgh for the standard LST optics") + effective_focal_length = OPTICS.effective_focal_length + + src_dep_df_proton = get_source_dependent_parameters( + df_proton, config, effective_focal_length=effective_focal_length + ) df_proton = pd.concat([df_proton, src_dep_df_proton['on']], axis=1) @@ -536,7 +553,7 @@ def apply_models(dl1, reg_disp_vector=None, reg_disp_norm=None, cls_disp_sign=None, - focal_length=28 * u.m, + effective_focal_length=29.30565 * u.m, custom_config=None, ): """ @@ -562,7 +579,7 @@ def apply_models(dl1, cls_disp_sign: string | Path | bytes | sklearn.ensemble.RandomForestClassifier Path to the random forest filename or file or pre-loaded RandomForestClassifier object for disp sign reconstruction - focal_length: `astropy.unit` + effective_focal_length: `astropy.unit` custom_config: dictionary Modified configuration to update the standard one @@ -654,7 +671,7 @@ def apply_models(dl1, dl2.y.values * u.m, dl2.reco_disp_dx.values * u.m, dl2.reco_disp_dy.values * u.m, - focal_length, + effective_focal_length, alt_tel * u.rad, az_tel * u.rad) @@ -688,7 +705,7 @@ def apply_models(dl1, return dl2 -def get_source_dependent_parameters(data, config, focal_length=28 * u.m): +def get_source_dependent_parameters(data, config, effective_focal_length=29.30565 * u.m): """Get parameters dict for source-dependent analysis. Parameters @@ -704,8 +721,9 @@ def get_source_dependent_parameters(data, config, focal_length=28 * u.m): else: data_type = 'real_data' - expected_src_pos_x_m, expected_src_pos_y_m = get_expected_source_pos(data, data_type, config, - focal_length=focal_length) + expected_src_pos_x_m, expected_src_pos_y_m = get_expected_source_pos( + data, data_type, config, effective_focal_length=effective_focal_length + ) src_dep_params = calc_source_dependent_parameters(data, expected_src_pos_x_m, expected_src_pos_y_m) src_dep_params_dict = {'on': src_dep_params} @@ -755,7 +773,7 @@ def calc_source_dependent_parameters(data, expected_src_pos_x_m, expected_src_po return src_dep_params -def get_expected_source_pos(data, data_type, config, focal_length=28 * u.m): +def get_expected_source_pos(data, data_type, config, effective_focal_length=29.30565 * u.m): """Get expected source position for source-dependent analysis . Parameters @@ -775,7 +793,7 @@ def get_expected_source_pos(data, data_type, config, focal_length=28 * u.m): expected_src_pos = utils.sky_to_camera( u.Quantity(data['mc_alt_tel'].values + config['mc_nominal_source_x_deg'], u.deg, copy=False), u.Quantity(data['mc_az_tel'].values + config['mc_nominal_source_y_deg'], u.deg, copy=False), - focal_length, + effective_focal_length, u.Quantity(data['mc_alt_tel'].values, u.deg, copy=False), u.Quantity(data['mc_az_tel'].values, u.deg, copy=False) ) @@ -806,7 +824,7 @@ def get_expected_source_pos(data, data_type, config, focal_length=28 * u.m): obstime = Time(time, scale='utc', format='unix') pointing_alt = u.Quantity(data['alt_tel'], u.rad, copy=False) pointing_az = u.Quantity(data['az_tel'], u.rad, copy=False) - source_pos = utils.radec_to_camera(source_coord, obstime, pointing_alt, pointing_az, focal_length) + source_pos = utils.radec_to_camera(source_coord, obstime, pointing_alt, pointing_az, effective_focal_length) expected_src_pos_x_m = source_pos.x.to_value(u.m) expected_src_pos_y_m = source_pos.y.to_value(u.m) diff --git a/lstchain/scripts/lstchain_add_source_dependent_parameters.py b/lstchain/scripts/lstchain_add_source_dependent_parameters.py index b5054561ca..d633e8c790 100644 --- a/lstchain/scripts/lstchain_add_source_dependent_parameters.py +++ b/lstchain/scripts/lstchain_add_source_dependent_parameters.py @@ -18,6 +18,7 @@ import pandas as pd from ctapipe.instrument import SubarrayDescription +from ctapipe_io_lst import OPTICS from lstchain.io import ( get_standard_config, @@ -61,11 +62,17 @@ def main(): pass dl1_params = pd.read_hdf(dl1_filename, key=dl1_params_lstcam_key) - subarray_info = SubarrayDescription.from_hdf(dl1_filename) - tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 - focal_length = subarray_info.tel[tel_id].optics.equivalent_focal_length - - src_dep_df = pd.concat(get_source_dependent_parameters(dl1_params, config, focal_length=focal_length), axis=1) + try: + subarray_info = SubarrayDescription.from_hdf(dl1_filename) + tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 + effective_focal_length = subarray_info.tel[tel_id].optics.effective_focal_length + except OSError: + print("subarray table is not readable because of the version inompatibility.") + print("Use the effective focal lentgh for the standard LST optics") + effective_focal_length = OPTICS.effective_focal_length + + src_dep_df = pd.concat(get_source_dependent_parameters( + dl1_params, config, effective_focal_length=effective_focal_length), axis=1) metadata = global_metadata() write_dataframe(src_dep_df, dl1_filename, dl1_params_src_dep_lstcam_key, config=config, meta=metadata) diff --git a/lstchain/scripts/lstchain_dl1_to_dl2.py b/lstchain/scripts/lstchain_dl1_to_dl2.py index 4b78c6221d..8e98b18c0b 100644 --- a/lstchain/scripts/lstchain_dl1_to_dl2.py +++ b/lstchain/scripts/lstchain_dl1_to_dl2.py @@ -13,6 +13,7 @@ import numpy as np import pandas as pd from ctapipe.instrument import SubarrayDescription +from ctapipe_io_lst import OPTICS from tables import open_file from lstchain.io import ( @@ -95,9 +96,14 @@ def apply_to_file(filename, models_dict, output_dir, config): data.alt_tel = - np.pi / 2. data.az_tel = - np.pi / 2. - subarray_info = SubarrayDescription.from_hdf(filename) - tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 - focal_length = subarray_info.tel[tel_id].optics.equivalent_focal_length + try: + subarray_info = SubarrayDescription.from_hdf(filename) + tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 + effective_focal_length = subarray_info.tel[tel_id].optics.effective_focal_length + except OSError: + print("subarray table is not readable because of the version inompatibility.") + print("Use the effective focal lentgh for the standard LST optics") + effective_focal_length = OPTICS.effective_focal_length # Apply the models to the data @@ -116,7 +122,7 @@ def apply_to_file(filename, models_dict, output_dir, config): models_dict['cls_gh'], models_dict['reg_energy'], reg_disp_vector=models_dict['disp_vector'], - focal_length=focal_length, + effective_focal_length=effective_focal_length, custom_config=config) elif config['disp_method'] == 'disp_norm_sign': dl2 = dl1_to_dl2.apply_models(data, @@ -124,7 +130,7 @@ def apply_to_file(filename, models_dict, output_dir, config): models_dict['reg_energy'], reg_disp_norm=models_dict['disp_norm'], cls_disp_sign=models_dict['disp_sign'], - focal_length=focal_length, + effective_focal_length=effective_focal_length, custom_config=config) # Source-dependent analysis @@ -136,7 +142,7 @@ def apply_to_file(filename, models_dict, output_dir, config): # if not, source-dependent parameters are added now else: data_srcdep = pd.concat(dl1_to_dl2.get_source_dependent_parameters( - data, config, focal_length=focal_length), axis=1) + data, config, effective_focal_length=effective_focal_length), axis=1) dl2_srcdep_dict = {} srcindep_keys = data.keys() @@ -157,7 +163,7 @@ def apply_to_file(filename, models_dict, output_dir, config): models_dict['cls_gh'], models_dict['reg_energy'], reg_disp_vector=models_dict['disp_vector'], - focal_length=focal_length, + effective_focal_length=effective_focal_length, custom_config=config) elif config['disp_method'] == 'disp_norm_sign': dl2_df = dl1_to_dl2.apply_models(data_with_srcdep_param, @@ -165,7 +171,7 @@ def apply_to_file(filename, models_dict, output_dir, config): models_dict['reg_energy'], reg_disp_norm=models_dict['disp_norm'], cls_disp_sign=models_dict['disp_sign'], - focal_length=focal_length, + effective_focal_length=effective_focal_length, custom_config=config) dl2_srcdep = dl2_df.drop(srcindep_keys, axis=1) From 7c146d46b749679bf68bdf389cfc7817d3c1baac Mon Sep 17 00:00:00 2001 From: Seiya Nozaki Date: Mon, 3 Jul 2023 17:42:59 +0200 Subject: [PATCH 02/12] small bug fix, and update the value for pytest --- lstchain/scripts/lstchain_mc_rfperformance.py | 16 +++++++++++----- lstchain/tests/test_lstchain.py | 6 +++--- 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/lstchain/scripts/lstchain_mc_rfperformance.py b/lstchain/scripts/lstchain_mc_rfperformance.py index 5ec26b916d..898bd3bb6e 100644 --- a/lstchain/scripts/lstchain_mc_rfperformance.py +++ b/lstchain/scripts/lstchain_mc_rfperformance.py @@ -21,6 +21,7 @@ import matplotlib.pyplot as plt import pandas as pd from ctapipe.instrument import SubarrayDescription +from ctapipe_io_lst import OPTICS from lstchain.io import ( read_configuration_file, @@ -87,10 +88,15 @@ def main(): config = replace_config(standard_config, custom_config) - subarray_info = SubarrayDescription.from_hdf(args.gammatest) - tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 - focal_length = subarray_info.tel[tel_id].optics.equivalent_focal_length - + try: + subarray_info = SubarrayDescription.from_hdf(args.gammatest) + tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 + effective_focal_length = subarray_info.tel[tel_id].optics.equivalent_focal_length + except OSError: + print("subarray table is not readable because of the version inompatibility.") + print("Use the effective focal lentgh for the standard LST optics") + effective_focal_length = OPTICS.effective_focal_length + reg_energy, reg_disp_norm, cls_disp_sign, cls_gh = dl1_to_dl2.build_models( args.gammafile, args.protonfile, @@ -110,7 +116,7 @@ def main(): data = pd.concat([gammas, proton], ignore_index=True) dl2 = dl1_to_dl2.apply_models(data, cls_gh, reg_energy, reg_disp_norm=reg_disp_norm, - cls_disp_sign=cls_disp_sign, focal_length=focal_length, + cls_disp_sign=cls_disp_sign, effective_focal_length=effective_focal_length, custom_config=config) ####PLOT SOME RESULTS##### diff --git a/lstchain/tests/test_lstchain.py b/lstchain/tests/test_lstchain.py index 08dfecf702..ffd09aac02 100644 --- a/lstchain/tests/test_lstchain.py +++ b/lstchain/tests/test_lstchain.py @@ -196,7 +196,7 @@ def test_get_source_dependent_parameters_mc(simulated_dl1_file): assert (src_dep_df_gamma['on']['expected_src_y'] == dl1_params['src_y']).all() np.testing.assert_allclose( - src_dep_df_proton['on']['expected_src_x'], 0.195, atol=1e-2 + src_dep_df_proton['on']['expected_src_x'], 0.205, atol=1e-2 ) np.testing.assert_allclose( src_dep_df_proton['on']['expected_src_y'], 0., atol=1e-2 @@ -224,13 +224,13 @@ def test_get_source_dependent_parameters_observed(observed_dl1_files): assert (src_dep_df_on['on']['expected_src_y'] == 0).all() np.testing.assert_allclose( - src_dep_df_wobble['on']['expected_src_x'], -0.195, atol=1e-2 + src_dep_df_wobble['on']['expected_src_x'], -0.205, atol=1e-2 ) np.testing.assert_allclose( src_dep_df_wobble['on']['expected_src_y'], 0., atol=1e-2 ) np.testing.assert_allclose( - src_dep_df_wobble['off_180']['expected_src_x'], 0.195, atol=1e-2 + src_dep_df_wobble['off_180']['expected_src_x'], 0.205, atol=1e-2 ) np.testing.assert_allclose( src_dep_df_wobble['off_180']['expected_src_y'], 0., atol=1e-2 From 8543db816fd528f5d7cc073b10dd3c0b2645b62b Mon Sep 17 00:00:00 2001 From: Seiya Nozaki Date: Mon, 3 Jul 2023 18:41:59 +0200 Subject: [PATCH 03/12] Use logging instead of print --- lstchain/reco/dl1_to_dl2.py | 67 +++++++++++++++++++------------------ 1 file changed, 35 insertions(+), 32 deletions(-) diff --git a/lstchain/reco/dl1_to_dl2.py b/lstchain/reco/dl1_to_dl2.py index cadf9c4c73..1e9b6f4922 100644 --- a/lstchain/reco/dl1_to_dl2.py +++ b/lstchain/reco/dl1_to_dl2.py @@ -8,6 +8,7 @@ """ import os +import logging import astropy.units as u import joblib @@ -33,6 +34,8 @@ from ctapipe.instrument import SubarrayDescription from ctapipe_io_lst import OPTICS +logger = logging.getLogger(__name__) + __all__ = [ 'apply_models', 'build_models', @@ -68,15 +71,15 @@ def train_energy(train, custom_config=None): features = config['energy_regression_features'] model = RandomForestRegressor - print("Given features: ", features) - print("Number of events for training: ", train.shape[0]) - print("Training Random Forest Regressor for Energy Reconstruction...") + logger.info("Given features: ", features) + logger.info("Number of events for training: ", train.shape[0]) + logger.info("Training Random Forest Regressor for Energy Reconstruction...") reg = model(**energy_regression_args) reg.fit(train[features], train['log_mc_energy']) - print("Model {} trained!".format(model)) + logger.info("Model {} trained!".format(model)) return reg @@ -106,16 +109,16 @@ def train_disp_vector(train, custom_config=None, predict_features=None): features = config['disp_regression_features'] model = RandomForestRegressor - print("Given features: ", features) - print("Number of events for training: ", train.shape[0]) - print("Training model {} for disp vector regression".format(model)) + logger.info("Given features: ", features) + logger.info("Number of events for training: ", train.shape[0]) + logger.info("Training model {} for disp vector regression".format(model)) reg = model(**disp_regression_args) x = train[features] y = np.transpose([train[f] for f in predict_features]) reg.fit(x, y) - print("Model {} trained!".format(model)) + logger.info("Model {} trained!".format(model)) return reg @@ -140,16 +143,16 @@ def train_disp_norm(train, custom_config=None, predict_feature='disp_norm'): features = config['disp_regression_features'] model = RandomForestRegressor - print("Given features: ", features) - print("Number of events for training: ", train.shape[0]) - print("Training model {} for disp norm regression".format(model)) + logger.info("Given features: ", features) + logger.info("Number of events for training: ", train.shape[0]) + logger.info("Training model {} for disp norm regression".format(model)) reg = model(**disp_regression_args) x = train[features] y = np.transpose(train[predict_feature]) reg.fit(x, y) - print("Model {} trained!".format(model)) + logger.info("Model {} trained!".format(model)) return reg @@ -174,16 +177,16 @@ def train_disp_sign(train, custom_config=None, predict_feature='disp_sign'): features = config["disp_classification_features"] model = RandomForestClassifier - print("Given features: ", features) - print("Number of events for training: ", train.shape[0]) - print("Training model {} for disp sign classification".format(model)) + logger.info("Given features: ", features) + logger.info("Number of events for training: ", train.shape[0]) + logger.info("Training model {} for disp sign classification".format(model)) clf = model(**classification_args) x = train[features] y = np.transpose(train[predict_feature]) clf.fit(x, y) - print("Model {} trained!".format(model)) + logger.info("Model {} trained!".format(model)) return clf @@ -212,24 +215,24 @@ def train_reco(train, custom_config=None): disp_features = config['disp_regression_features'] model = RandomForestRegressor - print("Given energy_features: ", energy_features) - print("Number of events for training: ", train.shape[0]) - print("Training Random Forest Regressor for Energy Reconstruction...") + logger.info("Given energy_features: ", energy_features) + logger.info("Number of events for training: ", train.shape[0]) + logger.info("Training Random Forest Regressor for Energy Reconstruction...") reg_energy = model(**energy_regression_args) reg_energy.fit(train[energy_features], train['log_mc_energy']) - print("Random Forest trained!") - print("Given disp_features: ", disp_features) - print("Training Random Forest Regressor for disp_norm Reconstruction...") + logger.info("Random Forest trained!") + logger.info("Given disp_features: ", disp_features) + logger.info("Training Random Forest Regressor for disp_norm Reconstruction...") reg_disp = RandomForestRegressor(**disp_regression_args) reg_disp.fit(train[disp_features], train['disp_norm']) - print("Random Forest trained!") - print("Done!") + logger.info("Random Forest trained!") + logger.info("Done!") return reg_energy, reg_disp @@ -254,16 +257,16 @@ def train_sep(train, custom_config=None): features = config["particle_classification_features"] model = RandomForestClassifier - print("Given features: ", features) - print("Number of events for training: ", train.shape[0]) - print("Training Random Forest Classifier for", + logger.info("Given features: ", features) + logger.info("Number of events for training: ", train.shape[0]) + logger.info("Training Random Forest Classifier for", "Gamma/Hadron separation...") clf = model(**classification_args) clf.fit(train[features], train['mc_type']) - print("Random Forest trained!") + logger.info("Random Forest trained!") return clf @@ -355,8 +358,8 @@ def build_models(filegammas, fileprotons, tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 effective_focal_length = subarray_info.tel[tel_id].optics.effective_focal_length except OSError: - print("subarray table is not readable because of the version inompatibility.") - print("Use the effective focal lentgh for the standard LST optics") + logger.warning("subarray table is not readable because of the version inompatibility.") + logger.warning("Use the effective focal lentgh for the standard LST optics") effective_focal_length = OPTICS.effective_focal_length src_dep_df_gamma = get_source_dependent_parameters( @@ -376,8 +379,8 @@ def build_models(filegammas, fileprotons, tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 effective_focal_length = subarray_info.tel[tel_id].optics.effective_focal_length except OSError: - print("subarray table is not readable because of the version inompatibility.") - print("Use the effective focal lentgh for the standard LST optics") + logger.warning("subarray table is not readable because of the version inompatibility.") + logger.warning("Use the effective focal lentgh for the standard LST optics") effective_focal_length = OPTICS.effective_focal_length src_dep_df_proton = get_source_dependent_parameters( From b9b5b5550248f12d154907d5e4c6740d373f4111 Mon Sep 17 00:00:00 2001 From: SeiyaNozaki Date: Thu, 13 Jul 2023 12:45:18 +0000 Subject: [PATCH 04/12] update disp parameters using effective focal length to take into accout of aberration effect --- lstchain/reco/dl1_to_dl2.py | 70 ++++++++++++++++++++++++++++++++----- 1 file changed, 61 insertions(+), 9 deletions(-) diff --git a/lstchain/reco/dl1_to_dl2.py b/lstchain/reco/dl1_to_dl2.py index 1e9b6f4922..19244a19e0 100644 --- a/lstchain/reco/dl1_to_dl2.py +++ b/lstchain/reco/dl1_to_dl2.py @@ -47,6 +47,7 @@ 'train_energy', 'train_reco', 'train_sep', + 'update_disp' ] @@ -346,6 +347,19 @@ def build_models(filegammas, fileprotons, df_gamma = pd.read_hdf(filegammas, key=dl1_params_lstcam_key) df_proton = pd.read_hdf(fileprotons, key=dl1_params_lstcam_key) + # Update parameters related to target direction on camera frame for gamma MC + # taking into account of the abrration effect using effective focal length + try: + subarray_info = SubarrayDescription.from_hdf(filegammas) + tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 + effective_focal_length = subarray_info.tel[tel_id].optics.effective_focal_length + except OSError: + logger.warning("subarray table is not readable because of the version inompatibility.") + logger.warning("Use the effective focal lentgh for the standard LST optics") + effective_focal_length = OPTICS.effective_focal_length + + df_gamma = update_disp(df_gamma, effective_focal_length = effective_focal_length) + if config['source_dependent']: # if source-dependent parameters are already in dl1 data, just read those data # if not, source-dependent parameters are added here @@ -353,15 +367,6 @@ def build_models(filegammas, fileprotons, src_dep_df_gamma = get_srcdep_params(filegammas) else: - try: - subarray_info = SubarrayDescription.from_hdf(filegammas) - tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 - effective_focal_length = subarray_info.tel[tel_id].optics.effective_focal_length - except OSError: - logger.warning("subarray table is not readable because of the version inompatibility.") - logger.warning("Use the effective focal lentgh for the standard LST optics") - effective_focal_length = OPTICS.effective_focal_length - src_dep_df_gamma = get_source_dependent_parameters( df_gamma, config, effective_focal_length=effective_focal_length ) @@ -607,6 +612,13 @@ def apply_models(dl1, + config['disp_classification_features'], ) + # Update parameters related to target direction on camera frame for MC data + # taking into account of the abrration effect using effective focal length + is_simu = (dl2['mc_type'] >= 0).all() if 'mc_type' in dl2.columns else False + if is_simu: + dl2 = update_disp(dl2, effective_focal_length = effective_focal_length) + + # Reconstruction of Energy and disp_norm distance if isinstance(reg_energy, (str, bytes, Path)): reg_energy = joblib.load(reg_energy) @@ -838,3 +850,43 @@ def get_expected_source_pos(data, data_type, config, effective_focal_length=29.3 ) return expected_src_pos_x_m, expected_src_pos_y_m + + +def update_disp(data, effective_focal_length=29.30565 * u.m): + """Update disp parameters using effective focal length + + Parameters + ---------- + data: Pandas DataFrame + config: dictionnary containing configuration + """ + + source_pos_in_camera = utils.sky_to_camera( + u.Quantity(data['mc_alt'].values, u.rad, copy=False), + u.Quantity(data['mc_az'].values, u.rad, copy=False), + effective_focal_length, + u.Quantity(data['mc_alt_tel'].values, u.rad, copy=False), + u.Quantity(data['mc_az_tel'].values, u.rad, copy=False) + ) + + expected_src_pos_x_m = source_pos_in_camera.x.to_value(u.m) + expected_src_pos_y_m = source_pos_in_camera.y.to_value(u.m) + + data['src_x'] = expected_src_pos_x_m + data['src_y'] = expected_src_pos_y_m + + disp_dx, disp_dy, disp_norm, disp_angle, disp_sign = disp.disp( + data['x'].values, + data['y'].values, + expected_src_pos_x_m, + expected_src_pos_y_m, + data['psi'].values + ) + + data['disp_dx'] = disp_dx + data['disp_dy'] = disp_dy + data['disp_norm'] = disp_norm + data['disp_angle'] = disp_angle + data['disp_sign'] = disp_sign + + return data From 709219de8b6e166c13a9b63dc76db69847b01497 Mon Sep 17 00:00:00 2001 From: SeiyaNozaki Date: Fri, 14 Jul 2023 08:01:16 +0000 Subject: [PATCH 05/12] small bug fix --- lstchain/reco/dl1_to_dl2.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/lstchain/reco/dl1_to_dl2.py b/lstchain/reco/dl1_to_dl2.py index 19244a19e0..e9bd436825 100644 --- a/lstchain/reco/dl1_to_dl2.py +++ b/lstchain/reco/dl1_to_dl2.py @@ -806,11 +806,11 @@ def get_expected_source_pos(data, data_type, config, effective_focal_length=29.3 # For proton MC, nominal source position is one written in config file if data_type == 'mc_proton': expected_src_pos = utils.sky_to_camera( - u.Quantity(data['mc_alt_tel'].values + config['mc_nominal_source_x_deg'], u.deg, copy=False), - u.Quantity(data['mc_az_tel'].values + config['mc_nominal_source_y_deg'], u.deg, copy=False), + u.Quantity(data['mc_alt_tel'].values + np.deg2rad(config['mc_nominal_source_x_deg']), u.rad, copy=False), + u.Quantity(data['mc_az_tel'].values + np.deg2rad(config['mc_nominal_source_y_deg']), u.rad, copy=False), effective_focal_length, - u.Quantity(data['mc_alt_tel'].values, u.deg, copy=False), - u.Quantity(data['mc_az_tel'].values, u.deg, copy=False) + u.Quantity(data['mc_alt_tel'].values, u.rad, copy=False), + u.Quantity(data['mc_az_tel'].values, u.rad, copy=False) ) expected_src_pos_x_m = expected_src_pos.x.to_value(u.m) From 7871b347135dce07474c8ce4be62c9fd8b20340f Mon Sep 17 00:00:00 2001 From: SeiyaNozaki Date: Fri, 14 Jul 2023 12:57:40 +0000 Subject: [PATCH 06/12] small update --- lstchain/reco/dl1_to_dl2.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/lstchain/reco/dl1_to_dl2.py b/lstchain/reco/dl1_to_dl2.py index e9bd436825..4a1713407b 100644 --- a/lstchain/reco/dl1_to_dl2.py +++ b/lstchain/reco/dl1_to_dl2.py @@ -47,7 +47,7 @@ 'train_energy', 'train_reco', 'train_sep', - 'update_disp' + 'update_disp_with_effective_focal_length' ] @@ -358,7 +358,7 @@ def build_models(filegammas, fileprotons, logger.warning("Use the effective focal lentgh for the standard LST optics") effective_focal_length = OPTICS.effective_focal_length - df_gamma = update_disp(df_gamma, effective_focal_length = effective_focal_length) + df_gamma = update_disp_with_effective_focal_length(df_gamma, effective_focal_length = effective_focal_length) if config['source_dependent']: # if source-dependent parameters are already in dl1 data, just read those data @@ -614,9 +614,9 @@ def apply_models(dl1, # Update parameters related to target direction on camera frame for MC data # taking into account of the abrration effect using effective focal length - is_simu = (dl2['mc_type'] >= 0).all() if 'mc_type' in dl2.columns else False + is_simu = 'disp_norm' in dl2.columns if is_simu: - dl2 = update_disp(dl2, effective_focal_length = effective_focal_length) + dl2 = update_disp_with_effective_focal_length(dl2, effective_focal_length = effective_focal_length) # Reconstruction of Energy and disp_norm distance @@ -852,7 +852,7 @@ def get_expected_source_pos(data, data_type, config, effective_focal_length=29.3 return expected_src_pos_x_m, expected_src_pos_y_m -def update_disp(data, effective_focal_length=29.30565 * u.m): +def update_disp_with_effective_focal_length(data, effective_focal_length=29.30565 * u.m): """Update disp parameters using effective focal length Parameters From a34fcffea7bea89d41123580fe3b539551f303ba Mon Sep 17 00:00:00 2001 From: SeiyaNozaki Date: Fri, 14 Jul 2023 17:54:39 +0000 Subject: [PATCH 07/12] update the calculation of assumed source position for proton MC --- lstchain/reco/dl1_to_dl2.py | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/lstchain/reco/dl1_to_dl2.py b/lstchain/reco/dl1_to_dl2.py index 4a1713407b..6611247c1d 100644 --- a/lstchain/reco/dl1_to_dl2.py +++ b/lstchain/reco/dl1_to_dl2.py @@ -805,16 +805,8 @@ def get_expected_source_pos(data, data_type, config, effective_focal_length=29.3 # For proton MC, nominal source position is one written in config file if data_type == 'mc_proton': - expected_src_pos = utils.sky_to_camera( - u.Quantity(data['mc_alt_tel'].values + np.deg2rad(config['mc_nominal_source_x_deg']), u.rad, copy=False), - u.Quantity(data['mc_az_tel'].values + np.deg2rad(config['mc_nominal_source_y_deg']), u.rad, copy=False), - effective_focal_length, - u.Quantity(data['mc_alt_tel'].values, u.rad, copy=False), - u.Quantity(data['mc_az_tel'].values, u.rad, copy=False) - ) - - expected_src_pos_x_m = expected_src_pos.x.to_value(u.m) - expected_src_pos_y_m = expected_src_pos.y.to_value(u.m) + expected_src_pos_x_m = np.tan(np.deg2rad(config['mc_nominal_source_x_deg'])) * effective_focal_length + expected_src_pos_y_m = np.tan(np.deg2rad(config['mc_nominal_source_y_deg'])) * effective_focal_length # For real data if data_type == 'real_data': From fd6b0ca9ae9197714ad3bd6734483864f0053168 Mon Sep 17 00:00:00 2001 From: SeiyaNozaki Date: Wed, 19 Jul 2023 14:59:42 +0000 Subject: [PATCH 08/12] correct typo --- lstchain/reco/dl1_to_dl2.py | 8 ++++---- .../scripts/lstchain_add_source_dependent_parameters.py | 4 ++-- lstchain/scripts/lstchain_dl1_to_dl2.py | 4 ++-- lstchain/scripts/lstchain_mc_rfperformance.py | 4 ++-- 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/lstchain/reco/dl1_to_dl2.py b/lstchain/reco/dl1_to_dl2.py index 33e072700e..ccef57a44d 100644 --- a/lstchain/reco/dl1_to_dl2.py +++ b/lstchain/reco/dl1_to_dl2.py @@ -354,8 +354,8 @@ def build_models(filegammas, fileprotons, tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 effective_focal_length = subarray_info.tel[tel_id].optics.effective_focal_length except OSError: - logger.warning("subarray table is not readable because of the version inompatibility.") - logger.warning("Use the effective focal lentgh for the standard LST optics") + logger.warning("subarray table is not readable because of the version incompatibility.") + logger.warning("The effective focal length for the standard LST optics will be used.") effective_focal_length = OPTICS.effective_focal_length df_gamma = update_disp_with_effective_focal_length(df_gamma, effective_focal_length = effective_focal_length) @@ -384,8 +384,8 @@ def build_models(filegammas, fileprotons, tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 effective_focal_length = subarray_info.tel[tel_id].optics.effective_focal_length except OSError: - logger.warning("subarray table is not readable because of the version inompatibility.") - logger.warning("Use the effective focal lentgh for the standard LST optics") + logger.warning("subarray table is not readable because of the version incompatibility.") + logger.warning("The effective focal length for the standard LST optics will be used.") effective_focal_length = OPTICS.effective_focal_length src_dep_df_proton = get_source_dependent_parameters( diff --git a/lstchain/scripts/lstchain_add_source_dependent_parameters.py b/lstchain/scripts/lstchain_add_source_dependent_parameters.py index d633e8c790..b63488f1d0 100644 --- a/lstchain/scripts/lstchain_add_source_dependent_parameters.py +++ b/lstchain/scripts/lstchain_add_source_dependent_parameters.py @@ -67,8 +67,8 @@ def main(): tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 effective_focal_length = subarray_info.tel[tel_id].optics.effective_focal_length except OSError: - print("subarray table is not readable because of the version inompatibility.") - print("Use the effective focal lentgh for the standard LST optics") + print("subarray table is not readable because of the version incompatibility.") + print("The effective focal length for the standard LST optics will be used.") effective_focal_length = OPTICS.effective_focal_length src_dep_df = pd.concat(get_source_dependent_parameters( diff --git a/lstchain/scripts/lstchain_dl1_to_dl2.py b/lstchain/scripts/lstchain_dl1_to_dl2.py index 94d3551290..ac94ebc356 100644 --- a/lstchain/scripts/lstchain_dl1_to_dl2.py +++ b/lstchain/scripts/lstchain_dl1_to_dl2.py @@ -103,8 +103,8 @@ def apply_to_file(filename, models_dict, output_dir, config): tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 effective_focal_length = subarray_info.tel[tel_id].optics.effective_focal_length except OSError: - print("subarray table is not readable because of the version inompatibility.") - print("Use the effective focal lentgh for the standard LST optics") + print("subarray table is not readable because of the version incompatibility.") + print("The effective focal length for the standard LST optics will be used.") effective_focal_length = OPTICS.effective_focal_length # Normalize all azimuth angles to the range [0, 360) degrees diff --git a/lstchain/scripts/lstchain_mc_rfperformance.py b/lstchain/scripts/lstchain_mc_rfperformance.py index 898bd3bb6e..179bf692b1 100644 --- a/lstchain/scripts/lstchain_mc_rfperformance.py +++ b/lstchain/scripts/lstchain_mc_rfperformance.py @@ -93,8 +93,8 @@ def main(): tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 effective_focal_length = subarray_info.tel[tel_id].optics.equivalent_focal_length except OSError: - print("subarray table is not readable because of the version inompatibility.") - print("Use the effective focal lentgh for the standard LST optics") + print("subarray table is not readable because of the version incompatibility.") + print("The effective focal length for the standard LST optics will be used.") effective_focal_length = OPTICS.effective_focal_length reg_energy, reg_disp_norm, cls_disp_sign, cls_gh = dl1_to_dl2.build_models( From 094560162f26b0065f6888a0146d8b3d2e94cad8 Mon Sep 17 00:00:00 2001 From: SeiyaNozaki Date: Wed, 19 Jul 2023 15:15:43 +0000 Subject: [PATCH 09/12] recalculate MC gamma source position inside get_expected_source_pos function --- lstchain/reco/dl1_to_dl2.py | 1 + 1 file changed, 1 insertion(+) diff --git a/lstchain/reco/dl1_to_dl2.py b/lstchain/reco/dl1_to_dl2.py index ccef57a44d..98365557d2 100644 --- a/lstchain/reco/dl1_to_dl2.py +++ b/lstchain/reco/dl1_to_dl2.py @@ -810,6 +810,7 @@ def get_expected_source_pos(data, data_type, config, effective_focal_length=29.3 # For gamma MC, expected source position is actual one for each event if data_type == 'mc_gamma': + data = update_disp_with_effective_focal_length(data, effective_focal_length = effective_focal_length) expected_src_pos_x_m = data['src_x'].values expected_src_pos_y_m = data['src_y'].values From 09dcf48754e7503956374825b31833459e36a425 Mon Sep 17 00:00:00 2001 From: SeiyaNozaki Date: Wed, 19 Jul 2023 15:19:53 +0000 Subject: [PATCH 10/12] remove lstchain_add_source_dependent_parameters.py --- ...stchain_add_source_dependent_parameters.py | 82 ------------------- .../scripts/tests/test_lstchain_scripts.py | 10 --- 2 files changed, 92 deletions(-) delete mode 100644 lstchain/scripts/lstchain_add_source_dependent_parameters.py diff --git a/lstchain/scripts/lstchain_add_source_dependent_parameters.py b/lstchain/scripts/lstchain_add_source_dependent_parameters.py deleted file mode 100644 index b63488f1d0..0000000000 --- a/lstchain/scripts/lstchain_add_source_dependent_parameters.py +++ /dev/null @@ -1,82 +0,0 @@ -#!/usr/bin/env python3 - -""" -Script to add the source dependent parameters to a DL1 file. - -Input: DL1 data file. Source dependent parameters will be added to this file. - -Usage: - -$> python lstchain_add_source_dependent_parameters.py ---input-file dl1_LST-1.Run02033.0137.h5 ---config lstchain_src_dep_config.json - -""" - -import argparse -import os - -import pandas as pd -from ctapipe.instrument import SubarrayDescription -from ctapipe_io_lst import OPTICS - -from lstchain.io import ( - get_standard_config, - read_configuration_file, -) -from lstchain.io.io import ( - dl1_params_lstcam_key, - dl1_params_src_dep_lstcam_key, - global_metadata, - write_dataframe, -) -from lstchain.reco.dl1_to_dl2 import get_source_dependent_parameters - -parser = argparse.ArgumentParser(description="Add the source dependent parameters to a DL1 file") - -# Required arguments -parser.add_argument('--input-file', '-f', type=str, - dest='input_file', - help='path to a DL1 HDF5 file', - ) - -# Optional arguments -parser.add_argument('--config', '-c', action='store', type=str, - dest='config_file', - help='Path to a configuration file for source dependent analysis', - default=None - ) - - -def main(): - - args = parser.parse_args() - - dl1_filename = os.path.abspath(args.input_file) - - config = get_standard_config() - if args.config_file is not None: - try: - config = read_configuration_file(os.path.abspath(args.config_file)) - except("Custom configuration could not be loaded !!!"): - pass - - dl1_params = pd.read_hdf(dl1_filename, key=dl1_params_lstcam_key) - try: - subarray_info = SubarrayDescription.from_hdf(dl1_filename) - tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 - effective_focal_length = subarray_info.tel[tel_id].optics.effective_focal_length - except OSError: - print("subarray table is not readable because of the version incompatibility.") - print("The effective focal length for the standard LST optics will be used.") - effective_focal_length = OPTICS.effective_focal_length - - src_dep_df = pd.concat(get_source_dependent_parameters( - dl1_params, config, effective_focal_length=effective_focal_length), axis=1) - - metadata = global_metadata() - write_dataframe(src_dep_df, dl1_filename, dl1_params_src_dep_lstcam_key, config=config, meta=metadata) - - -if __name__ == '__main__': - main() diff --git a/lstchain/scripts/tests/test_lstchain_scripts.py b/lstchain/scripts/tests/test_lstchain_scripts.py index 3e7d00d90c..1d39b6029a 100644 --- a/lstchain/scripts/tests/test_lstchain_scripts.py +++ b/lstchain/scripts/tests/test_lstchain_scripts.py @@ -65,16 +65,6 @@ def simulated_dl1ab(temp_dir_simulated_files, simulated_dl1_file): run_program("lstchain_dl1ab", "-f", simulated_dl1_file, "-o", output_file) return output_file -def test_add_source_dependent_parameters(temp_dir_simulated_srcdep_files, simulated_dl1_file): - shutil.copy(simulated_dl1_file, temp_dir_simulated_srcdep_files / "dl1_copy.h5") - dl1_file = temp_dir_simulated_srcdep_files / "dl1_copy.h5" - run_program("lstchain_add_source_dependent_parameters", "-f", dl1_file) - dl1_params_src_dep = get_srcdep_params(dl1_file) - - assert 'alpha' in dl1_params_src_dep['on'].columns - assert 'dist' in dl1_params_src_dep['on'].columns - assert 'time_gradient_from_source' in dl1_params_src_dep['on'].columns - assert 'skewness_from_source' in dl1_params_src_dep['on'].columns @pytest.fixture(scope="session") def merged_simulated_dl1_file(simulated_dl1_file, temp_dir_simulated_files): From 4b95835deda6d3c7ade1292d6cd046edfa709939 Mon Sep 17 00:00:00 2001 From: SeiyaNozaki Date: Wed, 19 Jul 2023 16:10:17 +0000 Subject: [PATCH 11/12] update for the calculation of proton MC expected source position --- lstchain/reco/dl1_to_dl2.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/lstchain/reco/dl1_to_dl2.py b/lstchain/reco/dl1_to_dl2.py index 98365557d2..fec9342004 100644 --- a/lstchain/reco/dl1_to_dl2.py +++ b/lstchain/reco/dl1_to_dl2.py @@ -32,6 +32,7 @@ from ctapipe.image.hillas import camera_to_shower_coordinates from ctapipe.instrument import SubarrayDescription +from ctapipe.coordinates import CameraFrame, TelescopeFrame from ctapipe_io_lst import OPTICS logger = logging.getLogger(__name__) @@ -816,8 +817,15 @@ def get_expected_source_pos(data, data_type, config, effective_focal_length=29.3 # For proton MC, nominal source position is one written in config file if data_type == 'mc_proton': - expected_src_pos_x_m = np.tan(np.deg2rad(config['mc_nominal_source_x_deg'])) * effective_focal_length - expected_src_pos_y_m = np.tan(np.deg2rad(config['mc_nominal_source_y_deg'])) * effective_focal_length + source_pos = SkyCoord( + fov_lon = -1 * config['mc_nominal_source_y_deg'] * u.deg, + fov_lat = config['mc_nominal_source_x_deg'] * u.deg, + frame=TelescopeFrame() + ) + camera_frame = CameraFrame(focal_length=effective_focal_length) + source_camera = source_pos.transform_to(camera_frame) + expected_src_pos_x_m = source_camera.x.to_value(u.m) + expected_src_pos_y_m = source_camera.y.to_value(u.m) # For real data if data_type == 'real_data': From 836030ba1b8164b1a839c32696adb91302b1df10 Mon Sep 17 00:00:00 2001 From: SeiyaNozaki Date: Wed, 19 Jul 2023 16:53:17 +0000 Subject: [PATCH 12/12] remove lstchain_add_source_dependent_parameters from index.rst --- docs/lstchain_api/scripts/index.rst | 17 +---------------- 1 file changed, 1 insertion(+), 16 deletions(-) diff --git a/docs/lstchain_api/scripts/index.rst b/docs/lstchain_api/scripts/index.rst index f087fa9587..d2006908c5 100644 --- a/docs/lstchain_api/scripts/index.rst +++ b/docs/lstchain_api/scripts/index.rst @@ -9,7 +9,6 @@ The scripts to be executed from the command line are described below: Currently both `scripts` and `Tools` are meant to be run from the command line. Please see also :ref:`tools` section for more information. -* `lstchain_add_source_dependent_parameters`_ * `lstchain_check_dl1`_ * `lstchain_create_run_summary`_ * `lstchain_data_create_time_calibration_file`_ @@ -30,20 +29,6 @@ The scripts to be executed from the command line are described below: * `lstchain_tune_nsb`_ -.. _lstchain_add_source_dependent_parameters: - -lstchain_add_source_dependent_parameters -++++++++++++++++++++++++++++++++++++++++ - -.. automodule:: lstchain.scripts.lstchain_add_source_dependent_parameters - -Usage ------ -.. argparse:: - :module: lstchain.scripts.lstchain_add_source_dependent_parameters - :func: parser - :prog: lstchain_add_source_dependent_parameters - .. _lstchain_check_dl1: lstchain_check_dl1 @@ -296,4 +281,4 @@ Usage .. argparse:: :module: lstchain.scripts.lstchain_dump_config :func: build_parser - :prog: lstchain_dump_config \ No newline at end of file + :prog: lstchain_dump_config