From b92a0577af25b3a479923f39313e215743314000 Mon Sep 17 00:00:00 2001 From: SeiyaNozaki Date: Mon, 3 Jul 2023 16:19:06 +0000 Subject: [PATCH 01/14] normalize all az_tel to [0, 360), and compute sin_az_tel if there is no sin_az_tel in DL1 data --- lstchain/reco/dl1_to_dl2.py | 11 +++++++++++ lstchain/scripts/lstchain_dl1_to_dl2.py | 9 +++++++++ 2 files changed, 20 insertions(+) diff --git a/lstchain/reco/dl1_to_dl2.py b/lstchain/reco/dl1_to_dl2.py index 116094898f..8b64c2432f 100644 --- a/lstchain/reco/dl1_to_dl2.py +++ b/lstchain/reco/dl1_to_dl2.py @@ -14,6 +14,7 @@ import numpy as np import pandas as pd from astropy.coordinates import SkyCoord +from astropy.coordinates import Angle from astropy.time import Time from pathlib import Path from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor @@ -392,6 +393,16 @@ def build_models(filegammas, fileprotons, + config['disp_classification_features'], ) + # Normalize all azimuth angles to the range [0, 360) degrees + df_gamma.az_tel = Angle(df_gamma.az_tel, u.rad).wrap_at(360 * u.deg).rad + df_proton.az_tel = Angle(df_proton.az_tel, u.rad).wrap_at(360 * u.deg).rad + + # Dealing with `sin_az_tel` missing data because of the former version of lstchain + if 'sin_az_tel' not in df_gamma.columns: + df_gamma.sin_az_tel = np.sin(df_gamma.az_tel) + if 'sin_az_tel' not in df_proton.columns: + df_proton.sin_az_tel = np.sin(df_proton.az_tel) + # Training MC gammas in reduced viewcone src_r_min = config['train_gamma_src_r_deg'][0] src_r_max = config['train_gamma_src_r_deg'][1] diff --git a/lstchain/scripts/lstchain_dl1_to_dl2.py b/lstchain/scripts/lstchain_dl1_to_dl2.py index 4b78c6221d..9deaf2124d 100644 --- a/lstchain/scripts/lstchain_dl1_to_dl2.py +++ b/lstchain/scripts/lstchain_dl1_to_dl2.py @@ -12,6 +12,8 @@ import numpy as np import pandas as pd +import astropy.units as u +from astropy.coordinates import Angle from ctapipe.instrument import SubarrayDescription from tables import open_file @@ -95,6 +97,13 @@ def apply_to_file(filename, models_dict, output_dir, config): data.alt_tel = - np.pi / 2. data.az_tel = - np.pi / 2. + # Normalize all azimuth angles to the range [0, 360) degrees + data.az_tel = Angle(data.az_tel, u.rad).wrap_at(360 * u.deg).rad + + # Dealing with `sin_az_tel` missing data because of the former version of lstchain + if 'sin_az_tel' not in data.columns: + data.sin_az_tel = np.sin(data.az_tel) + 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 From 35b5c01e65850889f6e9ca58a571165f61bf830e Mon Sep 17 00:00:00 2001 From: SeiyaNozaki Date: Mon, 3 Jul 2023 17:15:13 +0000 Subject: [PATCH 02/14] small bugfix --- lstchain/reco/dl1_to_dl2.py | 4 ++-- lstchain/scripts/lstchain_dl1_to_dl2.py | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/lstchain/reco/dl1_to_dl2.py b/lstchain/reco/dl1_to_dl2.py index 8b64c2432f..4ca0a18ef4 100644 --- a/lstchain/reco/dl1_to_dl2.py +++ b/lstchain/reco/dl1_to_dl2.py @@ -399,9 +399,9 @@ def build_models(filegammas, fileprotons, # Dealing with `sin_az_tel` missing data because of the former version of lstchain if 'sin_az_tel' not in df_gamma.columns: - df_gamma.sin_az_tel = np.sin(df_gamma.az_tel) + df_gamma['sin_az_tel'] = np.sin(df_gamma.az_tel) if 'sin_az_tel' not in df_proton.columns: - df_proton.sin_az_tel = np.sin(df_proton.az_tel) + df_proton['sin_az_tel'] = np.sin(df_proton.az_tel) # Training MC gammas in reduced viewcone src_r_min = config['train_gamma_src_r_deg'][0] diff --git a/lstchain/scripts/lstchain_dl1_to_dl2.py b/lstchain/scripts/lstchain_dl1_to_dl2.py index 9deaf2124d..a10fea90af 100644 --- a/lstchain/scripts/lstchain_dl1_to_dl2.py +++ b/lstchain/scripts/lstchain_dl1_to_dl2.py @@ -102,7 +102,7 @@ def apply_to_file(filename, models_dict, output_dir, config): # Dealing with `sin_az_tel` missing data because of the former version of lstchain if 'sin_az_tel' not in data.columns: - data.sin_az_tel = np.sin(data.az_tel) + data['sin_az_tel'] = np.sin(data.az_tel) subarray_info = SubarrayDescription.from_hdf(filename) tel_id = config["allowed_tels"][0] if "allowed_tels" in config else 1 From 5ed25b61a237be97b92e1320d60bcaa2e76a9623 Mon Sep 17 00:00:00 2001 From: SeiyaNozaki Date: Tue, 4 Jul 2023 09:17:52 +0000 Subject: [PATCH 03/14] combine astropy.coordinates imports --- lstchain/reco/dl1_to_dl2.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/lstchain/reco/dl1_to_dl2.py b/lstchain/reco/dl1_to_dl2.py index 4ca0a18ef4..a78d78245b 100644 --- a/lstchain/reco/dl1_to_dl2.py +++ b/lstchain/reco/dl1_to_dl2.py @@ -13,8 +13,7 @@ import joblib import numpy as np import pandas as pd -from astropy.coordinates import SkyCoord -from astropy.coordinates import Angle +from astropy.coordinates import SkyCoord, Angle from astropy.time import Time from pathlib import Path from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor From 63d61ab584ee29704bea2fd7079ad25d87789b95 Mon Sep 17 00:00:00 2001 From: Abelardo Moralejo Date: Tue, 4 Jul 2023 14:28:55 +0200 Subject: [PATCH 04/14] Change way to access monitoring info in dl1 file The former approach np.array(f.root[].cols.charge_mean) results now in an attempt to allocate a huge array (instead of just 2x2x1855), and hence a crash. No idea when or why the change of behaviour happened - could it be because of some change in pytables? --- lstchain/calib/camera/pixel_threshold_estimation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lstchain/calib/camera/pixel_threshold_estimation.py b/lstchain/calib/camera/pixel_threshold_estimation.py index 8ebfd390b0..c8cc98ad79 100644 --- a/lstchain/calib/camera/pixel_threshold_estimation.py +++ b/lstchain/calib/camera/pixel_threshold_estimation.py @@ -24,10 +24,10 @@ def get_bias_and_std(dl1_file): """ with tables.open_file(dl1_file) as f: ped = f.root[dl1_params_tel_mon_ped_key] - ped_charge_mean = np.array(ped.cols.charge_mean) - ped_charge_std = np.array(ped.cols.charge_std) + ped_charge_mean = ped.col('charge_mean') + ped_charge_std = ped.col('charge_std') calib = f.root[dl1_params_tel_mon_cal_key] - dc_to_pe = np.array(calib.cols.dc_to_pe[ORIGINAL_CALIBRATION_ID]) + dc_to_pe = calib.col('dc_to_pe')[ORIGINAL_CALIBRATION_ID] ped_charge_mean_pe = ped_charge_mean * dc_to_pe ped_charge_std_pe = ped_charge_std * dc_to_pe From 7d17ec05f9f0a3f2dbce95ef2a7ec3fbaf0b91d1 Mon Sep 17 00:00:00 2001 From: Ruben Lopez Coto Date: Thu, 6 Jul 2023 09:47:00 +0200 Subject: [PATCH 05/14] bump ctapipe_io_lst version --- environment.yml | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/environment.yml b/environment.yml index b7043052da..2209eb9ba1 100644 --- a/environment.yml +++ b/environment.yml @@ -31,7 +31,7 @@ dependencies: - pandas>=2 - pymongo - seaborn - - ctapipe_io_lst=0.21.0 + - ctapipe_io_lst=0.22.0 - pytest - pip: - ctaplot~=0.6.2 diff --git a/setup.py b/setup.py index 7c46513481..f890ce3b2e 100644 --- a/setup.py +++ b/setup.py @@ -42,7 +42,7 @@ def find_scripts(script_dir, prefix): 'astropy~=5.0', 'bokeh~=2.0', 'ctapipe~=0.19.2', - 'ctapipe_io_lst~=0.21.0', + 'ctapipe_io_lst~=0.22.0', 'ctaplot~=0.6.2', 'eventio>=1.9.1,<2.0.0a0', # at least 1.1.1, but not 2 'gammapy~=1.1', From dceb7d007fbcfd14bae7e5d6fe76508d0bd7a4a6 Mon Sep 17 00:00:00 2001 From: Ruben Lopez-Coto Date: Thu, 6 Jul 2023 13:09:17 +0200 Subject: [PATCH 06/14] Update environment.yml --- environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index 69c2785b43..5e6a594607 100644 --- a/environment.yml +++ b/environment.yml @@ -30,7 +30,7 @@ dependencies: - pandas>=2 - pymongo - seaborn - - ctapipe_io_lst=0.22.0 + - ctapipe_io_lst=0.22 - pytest - pip: - ctaplot~=0.6.2 From 48de0b427e071ff8f2ca22f0355da6b4789b4ecf Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Fri, 7 Jul 2023 10:06:27 +0200 Subject: [PATCH 07/14] Add missing v for release drafter template --- .github/release-drafter.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/release-drafter.yml b/.github/release-drafter.yml index 6cdad7a800..32bb1a8cfe 100644 --- a/.github/release-drafter.yml +++ b/.github/release-drafter.yml @@ -1,5 +1,5 @@ -name-template: '$NEXT_PATCH_VERSION' -tag-template: '$NEXT_PATCH_VERSION' +name-template: 'v$NEXT_PATCH_VERSION' +tag-template: 'v$NEXT_PATCH_VERSION' template: | ## What’s Changed From 0b965009e6a748efa2acf591591eb215ec28074c Mon Sep 17 00:00:00 2001 From: Daniel Morcuende Date: Fri, 7 Jul 2023 19:18:24 +0200 Subject: [PATCH 08/14] add option to overwite summary file --- lstchain/scripts/lstchain_create_run_summary.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/lstchain/scripts/lstchain_create_run_summary.py b/lstchain/scripts/lstchain_create_run_summary.py index f679fe44db..bc5b958fa0 100644 --- a/lstchain/scripts/lstchain_create_run_summary.py +++ b/lstchain/scripts/lstchain_create_run_summary.py @@ -54,6 +54,13 @@ default="/fefs/aswg/data/real/monitoring/RunSummary", ) +parser.add_argument( + "--overwrite", + action="store_true", + help="Overwrite existing Run Summary file", + default=False, +) + dtypes = { "ucts_timestamp": np.int64, "run_start": np.int64, @@ -306,7 +313,10 @@ def main(): run_summary.add_column(n_subruns, name="n_subruns", index=1) run_summary.add_column(run_types, name="run_type", index=2) run_summary.write( - args.output_dir / f"RunSummary_{args.date}.ecsv", format="ascii.ecsv", delimiter=",", + args.output_dir / f"RunSummary_{args.date}.ecsv", + format="ascii.ecsv", + delimiter=",", + overwrite=args.overwrite, ) From 7468fd00968af048c7532a8a295abbfe74258cde Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 11 Jul 2023 10:51:18 +0200 Subject: [PATCH 09/14] Get pyirf from conda in env file --- environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index 5e6a594607..7c8f397e74 100644 --- a/environment.yml +++ b/environment.yml @@ -32,6 +32,6 @@ dependencies: - seaborn - ctapipe_io_lst=0.22 - pytest + - pyirf=0.8 - pip: - ctaplot~=0.6.2 - - pyirf~=0.8.0 From 3592e084f123997183b174d0864364a5151d94c9 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 12 Jul 2023 15:01:54 +0000 Subject: [PATCH 10/14] Add detection QE spread for outlier estimation --- lstchain/calib/camera/calibration_calculator.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/lstchain/calib/camera/calibration_calculator.py b/lstchain/calib/camera/calibration_calculator.py index ad93c2b1de..75592cd320 100644 --- a/lstchain/calib/camera/calibration_calculator.py +++ b/lstchain/calib/camera/calibration_calculator.py @@ -210,10 +210,12 @@ def calculate_calibration_coefficients(self, event): # define unusables on number of estimated pe npe_deviation = calib_data.n_pe - npe_median[:,np.newaxis] - # cut on the base of pe statistical uncertanty (neglect the 7% spread due to detection QE) + # cut on the base of pe statistical uncertainty (adding a 7% spread due to different detection QE among PMs) + tot_std = np.sqrt(npe_median + (0.07 * npe_median)**2) + npe_outliers = ( - np.logical_or(npe_deviation < self.npe_median_cut_outliers[0] * np.sqrt(npe_median)[:,np.newaxis], - npe_deviation > self.npe_median_cut_outliers[1] * np.sqrt(npe_median)[:,np.newaxis])) + np.logical_or(npe_deviation < self.npe_median_cut_outliers[0] * tot_std[:,np.newaxis], + npe_deviation > self.npe_median_cut_outliers[1] * tot_std[:,np.newaxis])) # calibration unusable pixels are an OR of all masks calib_data.unusable_pixels = np.logical_or(unusable_pixels, npe_outliers) From dc9ca3503cb778a1d44b54d1391623e9a7180518 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Fri, 14 Jul 2023 12:08:52 +0000 Subject: [PATCH 11/14] Improve info and warning messages --- lstchain/calib/camera/flatfield.py | 1 + .../data/onsite_camera_calibration_param.json | 3 ++ .../onsite/onsite_create_calibration_file.py | 2 +- .../tools/lstchain_create_calibration_file.py | 28 +++++++++++++------ 4 files changed, 24 insertions(+), 10 deletions(-) diff --git a/lstchain/calib/camera/flatfield.py b/lstchain/calib/camera/flatfield.py index 4f7d5f6ae7..fdf40fd157 100644 --- a/lstchain/calib/camera/flatfield.py +++ b/lstchain/calib/camera/flatfield.py @@ -183,6 +183,7 @@ def calculate_relative_gain(self, event): self.collect_sample(charge, pixel_mask, arrival_time) sample_age = (self.trigger_time - self.time_start).to_value(u.s) + # check if to create a calibration event if (self.num_events_seen > 0 and (sample_age > self.sample_duration or diff --git a/lstchain/data/onsite_camera_calibration_param.json b/lstchain/data/onsite_camera_calibration_param.json index 656a848d7c..59b24d7e20 100644 --- a/lstchain/data/onsite_camera_calibration_param.json +++ b/lstchain/data/onsite_camera_calibration_param.json @@ -7,6 +7,9 @@ "log_level": "INFO" }, "LSTEventSource": { + "MultiFiles": { + "all_subruns": true + }, "allowed_tels": [ 1 ], diff --git a/lstchain/scripts/onsite/onsite_create_calibration_file.py b/lstchain/scripts/onsite/onsite_create_calibration_file.py index 2f8ccddb34..fd553152f5 100755 --- a/lstchain/scripts/onsite/onsite_create_calibration_file.py +++ b/lstchain/scripts/onsite/onsite_create_calibration_file.py @@ -236,7 +236,7 @@ def main(): f"--PedestalCalculator.sample_size={stat_events}", f"--config={config_file}", f"--log-file={log_file}", - "--log-file-level=DEBUG", + "--log-file-level=INFO", *remaining_args, ] diff --git a/lstchain/tools/lstchain_create_calibration_file.py b/lstchain/tools/lstchain_create_calibration_file.py index 05ab8c0fd9..f6f7485e60 100644 --- a/lstchain/tools/lstchain_create_calibration_file.py +++ b/lstchain/tools/lstchain_create_calibration_file.py @@ -103,10 +103,11 @@ def __init__(self, **kwargs): self.processor = None self.writer = None self.simulation = False + self.n_calib = 0 def setup(self): - self.log.debug("Opening file") + self.log.info("Opening file") self.eventsource = EventSource.from_config(parent=self) self.processor = CalibrationCalculator.from_name( @@ -130,7 +131,7 @@ def setup(self): group_name = 'tel_' + str(tel_id) - self.log.debug(f"Open output file {self.output_file}") + self.log.info(f"Open output file {self.output_file}") self.writer = HDF5TableWriter( filename=self.output_file, group_name=group_name, overwrite=True @@ -202,16 +203,16 @@ def start(self): # write flatfield results when enough statistics (also for pedestals) if (new_ff and new_ped): - self.log.debug(f"Write calibration at event n. {count+1}, event id {event.index.event_id} ") + self.log.info(f"Write calibration at event n. {count+1}, event id {event.index.event_id} ") - self.log.debug(f"Ready flatfield data at event n. {count_ff} " - f"stat = {ff_data.n_events} events") + self.log.info(f"Ready flatfield data at event n. {count_ff} " + f"stat = {ff_data.n_events} events") # write on file self.writer.write('flatfield', ff_data) - self.log.debug(f"Ready pedestal data at event n. {count_ped}" - f"stat = {ped_data.n_events} events") + self.log.info(f"Ready pedestal data at event n. {count_ped} " + f"stat = {ped_data.n_events} events") # write only pedestal data used for calibration self.writer.write('pedestal', ped_data) @@ -223,15 +224,24 @@ def start(self): self.processor.calculate_calibration_coefficients(event) # write calib and pixel status - self.log.debug("Write pixel_status data") + self.log.info("Write pixel_status data") self.writer.write('pixel_status', status_data) - self.log.debug("Write calibration data") + self.log.info("Write calibration data") self.writer.write('calibration', calib_data) + self.n_calib += 1 + if self.one_event: break def finish(self): + + self.log.info(f"Written {self.n_calib} calibration events") + if self.n_calib == 0: + self.log.warning(f"!!! No calibration events in the output file !!! : ") + self.log.warning(f"flatfiels collected statistics = {self.processor.flatfield.num_events_seen} events") + self.log.warning(f"pedestal collected statistics = {self.processor.pedestal.num_events_seen} events") + Provenance().add_output_file( self.output_file, role='mon.tel.calibration' From 2044145d43d2cfa703e0af427a8b27494aaee928 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 17 Jul 2023 08:15:56 +0000 Subject: [PATCH 12/14] React strongly in the case of missing calibration events --- lstchain/tools/lstchain_create_calibration_file.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/lstchain/tools/lstchain_create_calibration_file.py b/lstchain/tools/lstchain_create_calibration_file.py index f6f7485e60..850605ceba 100644 --- a/lstchain/tools/lstchain_create_calibration_file.py +++ b/lstchain/tools/lstchain_create_calibration_file.py @@ -238,9 +238,10 @@ def finish(self): self.log.info(f"Written {self.n_calib} calibration events") if self.n_calib == 0: - self.log.warning(f"!!! No calibration events in the output file !!! : ") - self.log.warning(f"flatfiels collected statistics = {self.processor.flatfield.num_events_seen} events") - self.log.warning(f"pedestal collected statistics = {self.processor.pedestal.num_events_seen} events") + self.log.critical(f"!!! No calibration events in the output file !!! : ") + self.log.critical(f"flatfiels collected statistics = {self.processor.flatfield.num_events_seen} events") + self.log.critical(f"pedestal collected statistics = {self.processor.pedestal.num_events_seen} events") + self.exit(1) Provenance().add_output_file( self.output_file, From 0fb2b572c5820db34581a30ebccb9458b90d291b Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 17 Jul 2023 09:11:00 +0000 Subject: [PATCH 13/14] Minor corrections --- lstchain/tools/lstchain_create_calibration_file.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lstchain/tools/lstchain_create_calibration_file.py b/lstchain/tools/lstchain_create_calibration_file.py index 850605ceba..357968d54a 100644 --- a/lstchain/tools/lstchain_create_calibration_file.py +++ b/lstchain/tools/lstchain_create_calibration_file.py @@ -238,8 +238,8 @@ def finish(self): self.log.info(f"Written {self.n_calib} calibration events") if self.n_calib == 0: - self.log.critical(f"!!! No calibration events in the output file !!! : ") - self.log.critical(f"flatfiels collected statistics = {self.processor.flatfield.num_events_seen} events") + self.log.critical("!!! No calibration events in the output file !!! : ") + self.log.critical(f"flatfield collected statistics = {self.processor.flatfield.num_events_seen} events") self.log.critical(f"pedestal collected statistics = {self.processor.pedestal.num_events_seen} events") self.exit(1) From 4092af3386a9062759951d4755fd938850d4aac8 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 17 Jul 2023 10:19:34 +0000 Subject: [PATCH 14/14] Make qe_dispersion configurable --- lstchain/calib/camera/calibration_calculator.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/lstchain/calib/camera/calibration_calculator.py b/lstchain/calib/camera/calibration_calculator.py index 75592cd320..1cb9a089ad 100644 --- a/lstchain/calib/camera/calibration_calculator.py +++ b/lstchain/calib/camera/calibration_calculator.py @@ -48,6 +48,11 @@ class CalibrationCalculator(Component): help='Excess noise factor squared: 1+ Var(gain)/Mean(Gain)**2' ).tag(config=True) + relative_qe_dispersion = traits.Float( + 0.07, + help='Relative (effective) quantum efficiency dispersion of PMs over the camera' + ).tag(config=True) + pedestal_product = traits.create_class_enum_trait( PedestalCalculator, default_value='PedestalIntegrator' @@ -211,7 +216,7 @@ def calculate_calibration_coefficients(self, event): npe_deviation = calib_data.n_pe - npe_median[:,np.newaxis] # cut on the base of pe statistical uncertainty (adding a 7% spread due to different detection QE among PMs) - tot_std = np.sqrt(npe_median + (0.07 * npe_median)**2) + tot_std = np.sqrt(npe_median + (self.relative_qe_dispersion * npe_median)**2) npe_outliers = ( np.logical_or(npe_deviation < self.npe_median_cut_outliers[0] * tot_std[:,np.newaxis],