From ff66286b84c5be6249fc580643e704ccf8b0e95d Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 15 Apr 2020 15:41:12 +0200 Subject: [PATCH 01/45] Add config for calibration calculator --- lstchain/data/lstchain_standard_config.json | 34 +++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/lstchain/data/lstchain_standard_config.json b/lstchain/data/lstchain_standard_config.json index 591ba65523..4243500e16 100644 --- a/lstchain/data/lstchain_standard_config.json +++ b/lstchain/data/lstchain_standard_config.json @@ -102,6 +102,7 @@ "gain_selector_config": { "threshold": 3500 }, + "LocalPeakWindowSum":{ "window_shift": 4, "window_width":8 @@ -115,5 +116,38 @@ "algorithm": null, "parameters": { } + }, + "calibration_product": "LSTCalibrationCalculator", + + "LSTCalibrationCalculator":{ + "minimum_hg_charge_median": 5000, + "maximum_lg_charge_std": 300, + "flatfield_product": "FlasherFlatFieldCalculator", + "pedestal_product": "PedestalIntegrator", + "PedestalCalculator":{ + "sample_size": 100, + "sample_duration":100000, + "tel_id":1, + "charge_median_cut_outliers": [-10,10], + "charge_std_cut_outliers": [-10,10], + "charge_product":"FixedWindowSum" + }, + "FlatFieldCalculator":{ + "sample_size": 100, + "sample_duration":100000, + "tel_id":1, + "charge_product":"LocalPeakWindowSum", + "charge_median_cut_outliers": [-0.5,0.5], + "charge_std_cut_outliers": [-10,10], + "time_cut_outliers": [2,38] + }, + "LocalPeakWindowSum":{ + "window_shift": 5, + "window_width":12 + }, + "FixedWindowSum":{ + "window_start": 12, + "window_width":12 + } } } From 48f353b2964542a25d484d8bd22235212f98ab19 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 15 Apr 2020 15:41:45 +0200 Subject: [PATCH 02/45] Add new Component with performs camera calibration calculations --- .../calib/camera/calibration_calculator.py | 251 ++++++++++++++++++ 1 file changed, 251 insertions(+) create mode 100644 lstchain/calib/camera/calibration_calculator.py diff --git a/lstchain/calib/camera/calibration_calculator.py b/lstchain/calib/camera/calibration_calculator.py new file mode 100644 index 0000000000..70bc6146cd --- /dev/null +++ b/lstchain/calib/camera/calibration_calculator.py @@ -0,0 +1,251 @@ +""" +Component for the estimation of the calibration coefficients events +""" + + +from abc import abstractmethod +import numpy as np +from ctapipe.core import Component +from ctapipe.core import traits +from ctapipe.core.traits import Int, Float, List +from lstchain.calib.camera.flatfield import FlatFieldCalculator +from lstchain.calib.camera.pedestals import PedestalCalculator + + +__all__ = [ + 'CalibrationCalculator', + 'LSTCalibrationCalculator' + +] + + +class CalibrationCalculator(Component): + """ + Parent class for the camera calibration calculators. + Fills the MonitoringCameraContainer on the base of calibration events + + Parameters + ---------- + + flatfield_calculator: lstchain.calib.camera.flatfield + The flatfield to use. If None, then FlatFieldCalculator + will be used by default. + + pedestal_calculator: lstchain.calib.camera.pedestal + The pedestal to use. If None, then + PedestalCalculator will be used by default. + + kwargs + + """ + + pedestal_product = traits.enum_trait( + PedestalCalculator, + default='PedestalIntegrator' + ) + + flatfield_product = traits.enum_trait( + FlatFieldCalculator, + default='FlasherFlatFieldCalculator' + ) + + classes = List([ + FlatFieldCalculator, + PedestalCalculator + ] + + traits.classes_with_traits(FlatFieldCalculator) + + traits.classes_with_traits(PedestalCalculator) + + ) + + def __init__( + self, + **kwargs + + ): + + """ + Parent class for the camera calibration calculators. + Fills the MonitoringCameraContainer on the base of calibration events + + Parameters + ---------- + + flatfield_calculator: lstchain.calib.camera.flatfield + The flatfield to use. If None, then FlatFieldCalculator + will be used by default. + + pedestal_calculator: lstchain.calib.camera.pedestal + The pedestal to use. If None, then + PedestalCalculator will be used by default. + + kwargs + + """ + + super().__init__(**kwargs) + + self.flatfield = FlatFieldCalculator.from_name( + self.flatfield_product, + **kwargs + ) + self.pedestal = PedestalCalculator.from_name( + self.pedestal_product, + **kwargs + ) + + msg = "tel_id not the same for all calibration components" + assert self.pedestal.tel_id == self.flatfield.tel_id, msg + + self.tel_id = self.flatfield.tel_id + + self.log.debug(f"{self.pedestal}") + self.log.debug(f"{self.flatfield}") + + + + + +class LSTCalibrationCalculator(CalibrationCalculator): + """ + Calibration calculator for LST camera + Fills the MonitoringCameraContainer on the base of calibration events + + Parameters: + ---------- + minimum_hg_charge_median : + Temporary cut on HG charge till the calibox TIB do not work + (default for filter 5.2) + + maximum_lg_charge_std + Temporary cut on LG std against Lidar events till the calibox TIB do not work + (default for filter 5.2) + """ + + minimum_hg_charge_median = Float( + 5000, + help='Temporary cut on HG charge till the calibox TIB do not work (default for filter 5.2)' + ).tag(config=True) + + maximum_lg_charge_std = Float( + 300, + help='Temporary cut on LG std against Lidar events till the calibox TIB do not work (default for filter 5.2) ' + ).tag(config=True) + + def __init__(self, **kwargs): + """ + Calibration calculator for LST camera + Fills the MonitoringCameraContainer on the base of calibration events + + Parameters: + ---------- + minimum_hg_charge_median : + Temporary cut on HG charge till the calibox TIB do not work + (default for filter 5.2) + + maximum_lg_charge_std + Temporary cut on LG std against Lidar events till the calibox TIB do not work + (default for filter 5.2) + + """ + super().__init__(**kwargs) + + + + + def calculate_calibration_coefficients(self, event): + """ + Calculate calibration coefficients from flatfield and pedestal statistics + associated to the present event + + Parameters + ---------- + + """ + + ped_data = event.mon.tel[self.tel_id].pedestal + ff_data = event.mon.tel[self.tel_id].flatfield + status_data = event.mon.tel[self.tel_id].pixel_status + calib_data = event.mon.tel[self.tel_id].calibration + + + # mask from pedestal and flat-field data + monitoring_unusable_pixels = np.logical_or(status_data.pedestal_failing_pixels, + status_data.flatfield_failing_pixels) + + # calibration unusable pixels are an OR of all masks + calib_data.unusable_pixels = np.logical_or(monitoring_unusable_pixels, + status_data.hardware_failing_pixels) + + # Extract calibration coefficients with F-factor method + # Assume fix F2 factor, F2=1+Var(gain)/Mean(Gain)**2 must be known from elsewhere + F2 = 1.2 + + # calculate photon-electrons + n_pe = F2 * (ff_data.charge_median - ped_data.charge_median) ** 2 / ( + ff_data.charge_std ** 2 - ped_data.charge_std ** 2) + + # fill WaveformCalibrationContainer (this is an example) + calib_data.time = ff_data.sample_time + calib_data.time_range = ff_data.sample_time_range + calib_data.n_pe = n_pe + + # find signal median of good pixels + masked_npe = np.ma.array(n_pe, mask=calib_data.unusable_pixels) + npe_signal_median = np.ma.median(masked_npe, axis=1) + + # Flat field factor + ff = np.ma.getdata(npe_signal_median[:, np.newaxis] / n_pe) + + calib_data.dc_to_pe = n_pe / (ff_data.charge_median - ped_data.charge_median) * ff + + # put the time around zero + camera_time_median = np.median(ff_data.time_median, axis=1) + calib_data.time_correction = -ff_data.relative_time_median - camera_time_median[:, np.newaxis] + + ped_extractor_name = self.config.get("PedestalCalculator").get("charge_product") + ped_width = self.config.get(ped_extractor_name).get("window_width") + calib_data.pedestal_per_sample = ped_data.charge_median / ped_width + + # put to zero unusable pixels + calib_data.dc_to_pe[calib_data.unusable_pixels] = 0 + calib_data.pedestal_per_sample[calib_data.unusable_pixels] = 0 + + # eliminate inf values id any (still necessary?) + calib_data.dc_to_pe[np.isinf(calib_data.dc_to_pe)] = 0 + + + + def process_interleaved(self, event): + """ + Process interleaved calibration events (pedestals and FF) + Parameters + ---------- + """ + new_calibration = False + new_ped = False + + # if pedestal event + if event.r1.tel[self.tel_id].trigger_type == 32: + + new_ped = self.pedestal.calculate_pedestals(event) + + + # if flat-field event: no calibration TIB for the moment, + # use a cut on the charge for ff events and on std for rejecting Magic Lidar events + elif event.r1.tel[self.tel_id].trigger_type == 4 or ( + np.median(np.sum(event.r1.tel[self.tel_id].waveform[0], axis=1)) + > self.minimum_hg_charge_median + and np.std(np.sum(event.r1.tel[self.tel_id].waveform[1], axis=1)) + < self.maximum_lg_charge_std): + + new_ff = self.flatfield.calculate_relative_gain(event) + + + # if new ff, calculate new calibration coefficients + if new_ff: + self.calculate_calibration_coefficients(event) + new_calibration=True + + return new_calibration, new_ped + From 5e4fe9767cf85f69f17457458988eee9eaa00cbd Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 15 Apr 2020 15:44:34 +0200 Subject: [PATCH 03/45] Change the name of the calibration tool: it works also for MC data --- lstchain/scripts/onsite/onsite_create_calibration_file.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lstchain/scripts/onsite/onsite_create_calibration_file.py b/lstchain/scripts/onsite/onsite_create_calibration_file.py index 907693ffe2..51c39aa3b4 100755 --- a/lstchain/scripts/onsite/onsite_create_calibration_file.py +++ b/lstchain/scripts/onsite/onsite_create_calibration_file.py @@ -105,7 +105,7 @@ def main(): if ff_calibration is 'yes': # run lstchain script - cmd = f"lstchain_data_create_calibration_file " \ + cmd = f"lstchain_create_calibration_file " \ f"--input_file={input_file} --output_file={output_file} --pedestal_file={pedestal_file} " \ f"--FlatFieldCalculator.sample_size={stat_events} --PedestalCalculator.sample_size={stat_events} " \ f"--EventSource.max_events={max_events} --config={config_file} > {log_file} 2>&1" From ac174a943ac9e32533f344300d13e125778f81d9 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 15 Apr 2020 15:50:21 +0200 Subject: [PATCH 04/45] Changed in order to use the new calibration calculator component --- .../lstchain_data_create_calibration_file.py | 139 +++++------------- 1 file changed, 36 insertions(+), 103 deletions(-) diff --git a/lstchain/tools/lstchain_data_create_calibration_file.py b/lstchain/tools/lstchain_data_create_calibration_file.py index ab49ad69b7..0f306e27b8 100644 --- a/lstchain/tools/lstchain_data_create_calibration_file.py +++ b/lstchain/tools/lstchain_data_create_calibration_file.py @@ -10,8 +10,7 @@ from ctapipe.core import Tool from ctapipe.io import EventSource from ctapipe.io.containers import PixelStatusContainer -from lstchain.calib.camera.flatfield import FlatFieldCalculator -from lstchain.calib.camera.pedestals import PedestalCalculator +from lstchain.calib.camera.calibration_calculator import CalibrationCalculator from lstchain.calib.camera import CameraR0Calibrator __all__ = [ @@ -24,7 +23,7 @@ class CalibrationHDF5Writer(Tool): Tool that generates a HDF5 file with camera calibration coefficients. This is just an example on how the monitoring containers can be filled using the calibration Components in calib/camera. - This example is based on an input file with interleaved pedestal and flat-field events + This example is based on an input file with pedestal and flat-field events For getting help run: python calc_camera_calibration.py --help @@ -34,15 +33,6 @@ class CalibrationHDF5Writer(Tool): name = "CalibrationHDF5Writer" description = "Generate a HDF5 file with camera calibration coefficients" - minimum_hg_charge_median = Float( - 5000, - help='Temporary cut on HG charge till the calibox TIB do not work (default for filter 5.2)' - ).tag(config=True) - - maximum_lg_charge_std = Float( - 300, - help='Temporary cut on LG std against Lidar events till the calibox TIB do not work (default for filter 5.2) ' - ).tag(config=True) one_event = Bool( False, @@ -59,14 +49,9 @@ class CalibrationHDF5Writer(Tool): help='Name of the log file' ).tag(config=True) - pedestal_product = traits.enum_trait( - PedestalCalculator, - default='PedestalIntegrator' - ) - - flatfield_product = traits.enum_trait( - FlatFieldCalculator, - default='FlasherFlatFieldCalculator' + calibration_product = traits.enum_trait( + CalibrationCalculator, + default='LSTCalibrationCalculator' ) r0calibrator_product =traits.enum_trait( @@ -80,18 +65,15 @@ class CalibrationHDF5Writer(Tool): log_file='CalibrationHDF5Writer.log_file', max_events='EventSource.max_events', pedestal_file= 'LSTR0Corrections.pedestal_path', - flatfield_product='CalibrationHDF5Writer.flatfield_product', - pedestal_product='CalibrationHDF5Writer.pedestal_product', + calibration_product='CalibrationHDF5Writer.calibration_product', r0calibrator_product='CalibrationHDF5Writer.r0calibrator_product', )) classes = List([EventSource, - FlatFieldCalculator, - PedestalCalculator + CalibrationCalculator ] + traits.classes_with_traits(CameraR0Calibrator) - + traits.classes_with_traits(FlatFieldCalculator) - + traits.classes_with_traits(PedestalCalculator) + + traits.classes_with_traits(CalibrationCalculator) ) @@ -106,12 +88,10 @@ def __init__(self, **kwargs): """ self.eventsource = None - self.flatfield = None - self.pedestal = None + self.processor = None self.container = None self.writer = None self.r0calibrator = None - self.tel_id = None self.tot_events = 0 self.simulation = False @@ -129,12 +109,8 @@ def setup(self): self.tot_events = self.eventsource.max_events self.simulation = True - self.flatfield = FlatFieldCalculator.from_name( - self.flatfield_product, - **kwargs - ) - self.pedestal = PedestalCalculator.from_name( - self.pedestal_product, + self.processor = CalibrationCalculator.from_name( + self.calibration_product, **kwargs ) @@ -144,12 +120,8 @@ def setup(self): **kwargs ) - msg = "tel_id not the same for all calibration components" - assert self.r0calibrator.tel_id == self.pedestal.tel_id == self.flatfield.tel_id, msg - self.tel_id = self.flatfield.tel_id - - group_name = 'tel_' + str(self.tel_id) + group_name = 'tel_' + str(self.processor.tel_id) self.log.debug(f"Open output file {self.output_file}") @@ -160,15 +132,16 @@ def setup(self): def start(self): '''Calibration coefficient calculator''' + tel_id = self.processor.tel_id new_ped = False new_ff = False end_of_file = False - try: + try: self.log.debug(f"Start loop") for count, event in enumerate(self.eventsource): - if count % 1000 == 0: + if count % 100 == 0: self.log.debug(f"Event {count}") # if last event write results @@ -180,53 +153,53 @@ def start(self): if count == 0: if self.simulation: + initialize_pixel_status(event.mon.tel[tel_id],event.r1.tel[tel_id].waveform.shape) - initialize_pixel_status(event.mon.tel[self.tel_id],event.r1.tel[self.tel_id].waveform.shape) - - ped_data = event.mon.tel[self.tel_id].pedestal + ped_data = event.mon.tel[tel_id].pedestal ped_data.meta['config'] = self.config - ff_data = event.mon.tel[self.tel_id].flatfield + ff_data = event.mon.tel[tel_id].flatfield ff_data.meta['config'] = self.config - status_data = event.mon.tel[self.tel_id].pixel_status + status_data = event.mon.tel[tel_id].pixel_status status_data.meta['config'] = self.config - calib_data = event.mon.tel[self.tel_id].calibration + calib_data = event.mon.tel[tel_id].calibration calib_data.meta['config'] = self.config # correct for low level calibration self.r0calibrator.calibrate(event) # reject event without trigger type - if event.r1.tel[self.tel_id].trigger_type == -1: + if event.r1.tel[tel_id].trigger_type == -1: continue # if pedestal event - if event.r1.tel[self.tel_id].trigger_type == 32 or ( + + if event.r1.tel[tel_id].trigger_type == 32 or ( self.simulation and - np.median(np.sum(event.r1.tel[self.tel_id].waveform[0], axis=1)) - < self.minimum_hg_charge_median): + np.median(np.sum(event.r1.tel[tel_id].waveform[0], axis=1)) + < self.processor.minimum_hg_charge_median): - new_ped = self.pedestal.calculate_pedestals(event) + new_ped = self.processor.pedestal.calculate_pedestals(event) # if flat-field event: no calibration TIB for the moment, # use a cut on the charge for ff events and on std for rejecting Magic Lidar events - elif event.r1.tel[self.tel_id].trigger_type == 4 or ( - np.median(np.sum(event.r1.tel[self.tel_id].waveform[0], axis=1)) - > self.minimum_hg_charge_median - and np.std(np.sum(event.r1.tel[self.tel_id].waveform[1], axis=1)) - < self.maximum_lg_charge_std): + elif event.r1.tel[tel_id].trigger_type == 4 or ( + np.median(np.sum(event.r1.tel[tel_id].waveform[0], axis=1)) + > self.processor.minimum_hg_charge_median + and np.std(np.sum(event.r1.tel[tel_id].waveform[1], axis=1)) + < self.processor.maximum_lg_charge_std): - new_ff = self.flatfield.calculate_relative_gain(event) + new_ff = self.processor.flatfield.calculate_relative_gain(event) # write pedestal results when enough statistics or end of file if new_ped or end_of_file: # update the monitoring container with the last statistics if end_of_file: - self.pedestal.store_results(event) + self.processor.pedestal.store_results(event) # write the event self.log.debug(f"Write pedestal data at event n. {count+1}, id {event.r0.event_id} " @@ -242,7 +215,7 @@ def start(self): # update the monitoring container with the last statistics if end_of_file: - self.flatfield.store_results(event) + self.processor.flatfield.store_results(event) self.log.debug(f"Write flatfield data at event n. {count+1}, id {event.r0.event_id} " f"stat = {ff_data.n_events} events") @@ -253,49 +226,9 @@ def start(self): new_ff = False # Then, calculate calibration coefficients - - # mask from pedestal and flat-field data - monitoring_unusable_pixels = np.logical_or(status_data.pedestal_failing_pixels, - status_data.flatfield_failing_pixels) - - # calibration unusable pixels are an OR of all masks - calib_data.unusable_pixels = np.logical_or(monitoring_unusable_pixels, - status_data.hardware_failing_pixels) - - # Extract calibration coefficients with F-factor method - # Assume fix F2 factor, F2=1+Var(gain)/Mean(Gain)**2 must be known from elsewhere - F2 = 1.2 - - # calculate photon-electrons - n_pe = F2 * (ff_data.charge_median - ped_data.charge_median) ** 2 / ( - ff_data.charge_std ** 2 - ped_data.charge_std ** 2) - - - - # fill WaveformCalibrationContainer (this is an example) - calib_data.time = ff_data.sample_time - calib_data.time_range = ff_data.sample_time_range - calib_data.n_pe = n_pe - - # find signal median of good pixels - masked_npe = np.ma.array(n_pe,mask=calib_data.unusable_pixels) - npe_signal_median = np.ma.median(masked_npe, axis=1) - - # Flat field factor - ff = npe_signal_median[:, np.newaxis]/n_pe - - calib_data.dc_to_pe = n_pe/(ff_data.charge_median-ped_data.charge_median)*ff - - # put the time around zero - camera_time_median = np.median(ff_data.time_median, axis=1) - calib_data.time_correction = -ff_data.relative_time_median - camera_time_median[:, np.newaxis] - - ped_extractor_name = self.config.get("PedestalCalculator").get("charge_product") - ped_width=self.config.get(ped_extractor_name).get("window_width") - calib_data.pedestal_per_sample = ped_data.charge_median/ped_width + self.processor.calculate_calibration_coefficients(event) # write calib and pixel status - self.log.debug(f"Write pixel_status data") self.writer.write('pixel_status',status_data) @@ -304,7 +237,7 @@ def start(self): if self.one_event: break - #self.writer.write('mon', event.mon.tel[self.tel_id]) + #self.writer.write('mon', event.mon.tel[tel_id]) except ValueError as e: self.log.error(e) From c13ee32cd3db22533e35c65e0fc878acd0ffc4ee Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 15 Apr 2020 15:53:15 +0200 Subject: [PATCH 05/45] Rename module because it works also for MC data --- .../tools/lstchain_create_calibration_file.py | 273 ++++++++++++++++++ 1 file changed, 273 insertions(+) create mode 100644 lstchain/tools/lstchain_create_calibration_file.py diff --git a/lstchain/tools/lstchain_create_calibration_file.py b/lstchain/tools/lstchain_create_calibration_file.py new file mode 100644 index 0000000000..0f306e27b8 --- /dev/null +++ b/lstchain/tools/lstchain_create_calibration_file.py @@ -0,0 +1,273 @@ +""" +Extract flat field coefficients from flasher data files. +""" +import numpy as np +from traitlets import Dict, List, Unicode, Float, Bool + + +from ctapipe.core import Provenance, traits +from ctapipe.io import HDF5TableWriter +from ctapipe.core import Tool +from ctapipe.io import EventSource +from ctapipe.io.containers import PixelStatusContainer +from lstchain.calib.camera.calibration_calculator import CalibrationCalculator +from lstchain.calib.camera import CameraR0Calibrator + +__all__ = [ + 'CalibrationHDF5Writer' +] + + +class CalibrationHDF5Writer(Tool): + """ + Tool that generates a HDF5 file with camera calibration coefficients. + This is just an example on how the monitoring containers can be filled using + the calibration Components in calib/camera. + This example is based on an input file with pedestal and flat-field events + + For getting help run: + python calc_camera_calibration.py --help + + """ + + name = "CalibrationHDF5Writer" + description = "Generate a HDF5 file with camera calibration coefficients" + + + one_event = Bool( + False, + help='Stop after first calibration event' + ).tag(config=True) + + output_file = Unicode( + 'calibration.hdf5', + help='Name of the output file' + ).tag(config=True) + + log_file = Unicode( + 'None', + help='Name of the log file' + ).tag(config=True) + + calibration_product = traits.enum_trait( + CalibrationCalculator, + default='LSTCalibrationCalculator' + ) + + r0calibrator_product =traits.enum_trait( + CameraR0Calibrator, + default='NullR0Calibrator' + ) + + aliases = Dict(dict( + input_file='EventSource.input_url', + output_file='CalibrationHDF5Writer.output_file', + log_file='CalibrationHDF5Writer.log_file', + max_events='EventSource.max_events', + pedestal_file= 'LSTR0Corrections.pedestal_path', + calibration_product='CalibrationHDF5Writer.calibration_product', + r0calibrator_product='CalibrationHDF5Writer.r0calibrator_product', + )) + + classes = List([EventSource, + CalibrationCalculator + ] + + traits.classes_with_traits(CameraR0Calibrator) + + traits.classes_with_traits(CalibrationCalculator) + + ) + + def __init__(self, **kwargs): + super().__init__(**kwargs) + """ + Tool that generates a HDF5 file with camera calibration coefficients. + Input file must contain interleaved pedestal and flat-field events + + For getting help run: + python calc_camera_calibration.py --help + + """ + self.eventsource = None + self.processor = None + self.container = None + self.writer = None + self.r0calibrator = None + self.tot_events = 0 + self.simulation = False + + def setup(self): + kwargs = dict(parent=self) + self.log.debug(f"Open file") + self.eventsource = EventSource.from_config(**kwargs) + + # if data remember how many event in the files + if "LSTEventSource" in str(type(self.eventsource)): + + self.tot_events = len(self.eventsource.multi_file) + self.log.debug(f"Input file has file {self.tot_events} events") + else: + self.tot_events = self.eventsource.max_events + self.simulation = True + + self.processor = CalibrationCalculator.from_name( + self.calibration_product, + **kwargs + ) + + if self.r0calibrator_product: + self.r0calibrator = CameraR0Calibrator.from_name( + self.r0calibrator_product, + **kwargs + ) + + + group_name = 'tel_' + str(self.processor.tel_id) + + self.log.debug(f"Open output file {self.output_file}") + + self.writer = HDF5TableWriter( + filename=self.output_file, group_name=group_name, overwrite=True + ) + + def start(self): + '''Calibration coefficient calculator''' + + tel_id = self.processor.tel_id + new_ped = False + new_ff = False + end_of_file = False + + try: + self.log.debug(f"Start loop") + for count, event in enumerate(self.eventsource): + + if count % 100 == 0: + self.log.debug(f"Event {count}") + + # if last event write results + if count == self.tot_events-1 or count == self.eventsource.max_events-1: + self.log.debug(f"Last event, count = {count}") + end_of_file = True + + # save the config, to be retrieved as data.meta['config'] + if count == 0: + + if self.simulation: + initialize_pixel_status(event.mon.tel[tel_id],event.r1.tel[tel_id].waveform.shape) + + ped_data = event.mon.tel[tel_id].pedestal + ped_data.meta['config'] = self.config + + ff_data = event.mon.tel[tel_id].flatfield + ff_data.meta['config'] = self.config + + status_data = event.mon.tel[tel_id].pixel_status + status_data.meta['config'] = self.config + + calib_data = event.mon.tel[tel_id].calibration + calib_data.meta['config'] = self.config + + # correct for low level calibration + self.r0calibrator.calibrate(event) + + # reject event without trigger type + if event.r1.tel[tel_id].trigger_type == -1: + continue + + # if pedestal event + + if event.r1.tel[tel_id].trigger_type == 32 or ( + self.simulation and + np.median(np.sum(event.r1.tel[tel_id].waveform[0], axis=1)) + < self.processor.minimum_hg_charge_median): + + new_ped = self.processor.pedestal.calculate_pedestals(event) + + + # if flat-field event: no calibration TIB for the moment, + # use a cut on the charge for ff events and on std for rejecting Magic Lidar events + elif event.r1.tel[tel_id].trigger_type == 4 or ( + np.median(np.sum(event.r1.tel[tel_id].waveform[0], axis=1)) + > self.processor.minimum_hg_charge_median + and np.std(np.sum(event.r1.tel[tel_id].waveform[1], axis=1)) + < self.processor.maximum_lg_charge_std): + + new_ff = self.processor.flatfield.calculate_relative_gain(event) + + # write pedestal results when enough statistics or end of file + if new_ped or end_of_file: + + # update the monitoring container with the last statistics + if end_of_file: + self.processor.pedestal.store_results(event) + + # write the event + self.log.debug(f"Write pedestal data at event n. {count+1}, id {event.r0.event_id} " + f"stat = {ped_data.n_events} events") + + # write on file + self.writer.write('pedestal', ped_data) + + new_ped = False + + # write flatfield results when enough statistics (also for pedestals) or end of file + if (new_ff and ped_data.n_events > 0) or end_of_file: + + # update the monitoring container with the last statistics + if end_of_file: + self.processor.flatfield.store_results(event) + + self.log.debug(f"Write flatfield data at event n. {count+1}, id {event.r0.event_id} " + f"stat = {ff_data.n_events} events") + + # write on file + self.writer.write('flatfield', ff_data) + + new_ff = False + + # Then, calculate calibration coefficients + self.processor.calculate_calibration_coefficients(event) + + # write calib and pixel status + self.log.debug(f"Write pixel_status data") + self.writer.write('pixel_status',status_data) + + self.log.debug(f"Write calibration data") + self.writer.write('calibration', calib_data) + if self.one_event: + break + + #self.writer.write('mon', event.mon.tel[tel_id]) + except ValueError as e: + self.log.error(e) + + def finish(self): + Provenance().add_output_file( + self.output_file, + role='mon.tel.calibration' + ) + self.writer.close() + +def initialize_pixel_status(mon_camera_container,shape): + """ + Initialize the pixel status container in the case of + simulation events (this should be done in the event source, but + added here for the moment) + """ + # initialize the container + status_container = PixelStatusContainer() + status_container.hardware_failing_pixels = np.zeros((shape[0],shape[1]), dtype=bool) + status_container.pedestal_failing_pixels = np.zeros((shape[0],shape[1]), dtype=bool) + status_container.flatfield_failing_pixels = np.zeros((shape[0],shape[1]), dtype=bool) + + mon_camera_container.pixel_status = status_container + + +def main(): + exe = CalibrationHDF5Writer() + + exe.run() + + +if __name__ == '__main__': + main() From 99cf8841c06c9a6d4b2bcd11ef7e53ee45eddbfd Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 15 Apr 2020 15:54:32 +0200 Subject: [PATCH 06/45] Minor changes --- .../data/onsite_camera_calibration_param.json | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/lstchain/data/onsite_camera_calibration_param.json b/lstchain/data/onsite_camera_calibration_param.json index 7406681950..b71f884c20 100644 --- a/lstchain/data/onsite_camera_calibration_param.json +++ b/lstchain/data/onsite_camera_calibration_param.json @@ -10,20 +10,20 @@ "r0calibrator_product": "LSTR0Corrections", "log_level":"DEBUG" }, - "PedestalCalculator":{ + "PedestalCalculator":{ "sample_duration":100000, "tel_id":1, - "charge_median_cut_outliers": [-10,10], - "charge_std_cut_outliers": [-10,10], - "charge_product":"FixedWindowSum" + "charge_median_cut_outliers": [-10,10], + "charge_std_cut_outliers": [-10,10], + "charge_product":"FixedWindowSum" }, - "FlatFieldCalculator":{ + "FlatFieldCalculator":{ "sample_duration":100000, - "tel_id":1, + "tel_id":1, "charge_product":"LocalPeakWindowSum", "charge_median_cut_outliers": [-0.5,0.5], - "charge_std_cut_outliers": [-10,10], - "time_cut_outliers": [2,38] + "charge_std_cut_outliers": [-10,10], + "time_cut_outliers": [2,38] }, "LSTR0Corrections": { "r1_sample_start": 2, From e1a6ac9faa8ca7beb1c53ea1f22158997c34fc07 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 15 Apr 2020 15:59:43 +0200 Subject: [PATCH 07/45] Initialize all montoring containers (changes for using interleaved calibration events) --- lstchain/calib/camera/calibrator.py | 32 +++++++++++++++++++---------- 1 file changed, 21 insertions(+), 11 deletions(-) diff --git a/lstchain/calib/camera/calibrator.py b/lstchain/calib/camera/calibrator.py index 26e0a69b94..0f804e4f78 100644 --- a/lstchain/calib/camera/calibrator.py +++ b/lstchain/calib/camera/calibrator.py @@ -116,19 +116,23 @@ def _initialize_correction(self): with HDF5TableReader(self.calibration_path) as h5_table: assert h5_table._h5file.isopen == True for telid in self.allowed_tels: - # read the calibration data for the moment only one event + # read the calibration data table = '/tel_' + str(telid) + '/calibration' next(h5_table.read(table, self.mon_data.tel[telid].calibration)) - dc_to_pe=self.mon_data.tel[telid].calibration.dc_to_pe - # put to zero unusable pixels - dc_to_pe[self.mon_data.tel[telid].calibration.unusable_pixels] = 0 - # eliminate inf values id any (should be done probably before) - dc_to_pe[np.isinf(dc_to_pe)] = 0 + # read pedestal data + table = '/tel_' + str(telid) + '/pedestal' + next(h5_table.read(table, self.mon_data.tel[telid].pedestal)) + + # read flat-field data + table = '/tel_' + str(telid) + '/flatfield' + next(h5_table.read(table, self.mon_data.tel[telid].flatfield)) # read the pixel_status container table = '/tel_' + str(telid) + '/pixel_status' next(h5_table.read(table, self.mon_data.tel[telid].pixel_status)) + + except: self.log.error(f"Problem in reading calibration file {self.calibration_path}") @@ -141,17 +145,22 @@ def _calibrate_dl0(self, event, telid): return event.dl0.event_id = event.r1.event_id - event.mon.tel[telid].calibration = self.mon_data.tel[telid].calibration - event.mon.tel[telid].pixel_status = self.mon_data.tel[telid].pixel_status + + # if not already done, initialize the event monitoring containers + if event.mon.tel[telid].calibration.dc_to_pe is None: + event.mon.tel[telid].calibration = self.mon_data.tel[telid].calibration + event.mon.tel[telid].flatfield = self.mon_data.tel[telid].flatfield + event.mon.tel[telid].pedestal = self.mon_data.tel[telid].pedestal + event.mon.tel[telid].pixel_status = self.mon_data.tel[telid].pixel_status # subtract the pedestal per sample (should we do it?) and multiply for the calibration coefficients # + event.dl0.tel[telid].waveform = ( + (event.r1.tel[telid].waveform - event.mon.tel[telid].calibration.pedestal_per_sample[:, :, np.newaxis]) + * event.mon.tel[telid].calibration.dc_to_pe[:, :, np.newaxis]) - event.dl0.tel[telid].waveform = ( - (event.r1.tel[telid].waveform - self.mon_data.tel[telid].calibration.pedestal_per_sample[:, :, np.newaxis]) - * self.mon_data.tel[telid].calibration.dc_to_pe[:, :, np.newaxis]) def _calibrate_dl1(self, event, telid): """ @@ -165,6 +174,7 @@ def _calibrate_dl1(self, event, telid): if self.image_extractor.requires_neighbors(): camera = event.inst.subarray.tel[telid].camera self.image_extractor.neighbors = camera.neighbor_matrix_where + charge, pulse_time = self.image_extractor(waveforms) # correct time with drs4 correction if available From d3b050ea2591e296a15dde41dfaf07b9d5e5d4e7 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 15 Apr 2020 16:02:31 +0200 Subject: [PATCH 08/45] Add function to write calibration events containers in DL1 --- lstchain/io/io.py | 33 ++++++++++++++++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/lstchain/io/io.py b/lstchain/io/io.py index 3bacfe6f44..be62baaa6c 100644 --- a/lstchain/io/io.py +++ b/lstchain/io/io.py @@ -33,7 +33,8 @@ 'write_subarray_tables', 'write_metadata', 'write_dataframe', - 'write_dl2_dataframe' + 'write_dl2_dataframe', + 'write_calibration_data' ] @@ -667,3 +668,33 @@ def recursive_copy_node(src_file, dir_file, path): recursive_path = os.path.join(recursive_path, p) +def write_calibration_data(writer,mon_index, mon_event, new_ped=False, new_calib=False): + mon_event.pedestal.prefix = '' + mon_event.flatfield.prefix = '' + mon_event.calibration.prefix = '' + mon_index.prefix = '' + + if new_ped: + # write ped container + writer.write( + table_name=f'telescope/monitoring/pedestal', + containers=[mon_index, mon_event.pedestal] + ) + # update index + mon_index.pedestal_id += 1 + + if new_calib: + # write calibration container + writer.write( + table_name="telescope/monitoring/flatfield", + containers=[mon_index, mon_event.flatfield] + ) + + # write ff container + writer.write( + table_name="telescope/monitoring/calibration", + containers=[mon_index, mon_event.calibration] + ) + # update index + mon_index.flatfield_id += 1 + mon_index.calibration_id += 1 \ No newline at end of file From fe1920629eda672edfad1d9d5eaba91e39e1de9b Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 15 Apr 2020 16:05:04 +0200 Subject: [PATCH 09/45] Add the processing and writing of the interleaved calibration data --- lstchain/io/lstcontainers.py | 28 +++++++++++---- lstchain/reco/r0_to_dl1.py | 68 ++++++++++++++++++++++++++++++------ 2 files changed, 79 insertions(+), 17 deletions(-) diff --git a/lstchain/io/lstcontainers.py b/lstchain/io/lstcontainers.py index bf473e196b..84f89f06f9 100644 --- a/lstchain/io/lstcontainers.py +++ b/lstchain/io/lstcontainers.py @@ -9,6 +9,7 @@ from ctapipe.image import timing_parameters as time from ctapipe.image import leakage from ctapipe.image.cleaning import number_of_islands + from ..reco import utils from numpy import nan @@ -16,9 +17,11 @@ 'DL1ParametersContainer', 'DispContainer', 'MetaData', - 'ThrownEventsHistogram' + 'ThrownEventsHistogram', + 'DL1MonitoringEventIndexContainer' ] + class DL1ParametersContainer(Container): """ TODO: maybe fields could be inherited from ctapipe containers definition @@ -53,6 +56,7 @@ class DL1ParametersContainer(Container): obs_id = Field(None, 'Observation ID') event_id = Field(None, 'Event ID') + calibration_id = Field(None, 'ID of the employed calibration event') gps_time = Field(None, 'GPS time event trigger') dragon_time = Field(None, 'Dragon time event trigger') ucts_time = Field(None, 'UCTS time event trigger') @@ -82,6 +86,14 @@ class DL1ParametersContainer(Container): tel_pos_y = Field(None, "Telescope y position in the ground") tel_pos_z = Field(None, "Telescope z position in the ground") + trigger_type = Field(None, "trigger type") + ucts_trigger_type = Field(None, "UCTS trigger type") + trigger_time = Field(None, "trigger time") + + # info not available in data + #num_trig_pix = Field(None, "Number of trigger groups (sectors) listed") + #trig_pix_id = Field(None, "pixels involved in the camera trigger") + def fill_hillas(self, hillas): """ fill Hillas parameters @@ -194,11 +206,7 @@ class ExtraMCInfo(Container): class ExtraImageInfo(Container): """ attach the tel_id """ tel_id = Field(0, "Telescope ID") - trigger_type = Field(None, "trigger type") - ucts_trigger_type = Field(None, "UCTS trigger type") - trigger_time = Field(None, "trigger time") - num_trig_pix = Field(None, "Number of trigger groups (sectors) listed") - trig_pix_id = Field(None, "pixels involved in the camera trigger") + low_gain_mask = Field(None, "pixels with charge from low gain channel") class ThrownEventsHistogram(Container): @@ -234,3 +242,11 @@ class MetaData(Container): CONTACT = Field(None, "Person or institution responsible for this data product") +class DL1MonitoringEventIndexContainer(Container): + """ + Container with the calibration coefficients + """ + tel_id = Field(1, 'Index of telescope') + calibration_id = Field(0, 'Index of calibration event for DL1 file') + pedestal_id = Field(0, 'Index of pedestal event for DL1 file') + flatfield_id = Field(0, 'Index of flat-field event for DL1 file') \ No newline at end of file diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py index ece4d38dd6..644f024a1c 100644 --- a/lstchain/reco/r0_to_dl1.py +++ b/lstchain/reco/r0_to_dl1.py @@ -24,8 +24,9 @@ import math from . import utils from .volume_reducer import check_and_apply_volume_reduction -from ..io.lstcontainers import ExtraImageInfo +from ..io.lstcontainers import ExtraImageInfo, DL1MonitoringEventIndexContainer from ..calib.camera import lst_calibration, load_calibrator_from_config +from ..calib.camera.calibration_calculator import CalibrationCalculator from ..io import DL1ParametersContainer, standard_config, replace_config from ..image.muon import analyze_muon_event, tag_pix_thr from ..image.muon import create_muon_table, fill_muon_event @@ -36,7 +37,7 @@ import tables from functools import partial from ..io import write_simtel_energy_histogram, write_mcheader, write_array_info, global_metadata -from ..io import add_global_metadata, write_metadata, write_subarray_tables +from ..io import add_global_metadata, write_metadata, write_subarray_tables, write_calibration_data from ..io.io import add_column_table import pandas as pd @@ -70,7 +71,7 @@ ) -def get_dl1(calibrated_event, telescope_id, dl1_container = None, +def get_dl1(calibrated_event, telescope_id, dl1_container = None, custom_config = {}, use_main_island = True): """ Return a DL1ParametersContainer of extracted features from a calibrated event. @@ -136,6 +137,7 @@ def get_dl1(calibrated_event, telescope_id, dl1_container = None, dl1_container.n_islands = num_islands dl1_container.set_telescope_info(calibrated_event, telescope_id) + return dl1_container else: @@ -241,7 +243,14 @@ def r0_to_dl1(input_filename = get_dataset_path('gamma_test_large.simtel.gz'), config = Config(config), allowed_tels = [1],) - + # Component to process interleaved pedestal and flat-fields + calibration_calculator = CalibrationCalculator.from_name( + config['calibration_product'], + config=Config(config[config['calibration_product']]) + ) + + + calibration_index = DL1MonitoringEventIndexContainer() dl1_container = DL1ParametersContainer() if pointing_file_path: @@ -290,6 +299,7 @@ def r0_to_dl1(input_filename = get_dataset_path('gamma_test_large.simtel.gz'), transform = tel_list_transform ) + ### EVENT LOOP ### for i, event in enumerate(source): if i % 100 == 0: @@ -299,6 +309,7 @@ def r0_to_dl1(input_filename = get_dataset_path('gamma_test_large.simtel.gz'), event.mc.prefix = 'mc' event.trig.prefix = '' + # write sub tables if is_simu: write_subarray_tables(writer, event, metadata) @@ -306,9 +317,40 @@ def r0_to_dl1(input_filename = get_dataset_path('gamma_test_large.simtel.gz'), cal_mc(event) else: + if i==0: + # iniitalize the telescope + # FIXME? LST calibrator is only for one telescope + # it should be inside the telescope loop + + tel_id = calibration_calculator.tel_id + + # write the first calibration event (initialized from h5 file) + write_calibration_data(writer, + calibration_index, + r1_dl1_calibrator.mon_data.tel[tel_id], + new_ped=True, new_calib=True) + + # drs4 calibrations r0_r1_calibrator.calibrate(event) + + # process interleaved events (pedestals, ff, calibration) + new_calib_event, new_ped_event = calibration_calculator.process_interleaved(event) + + # write monitoring containers if updated + if new_ped_event or new_calib_event: + write_calibration_data(writer, + calibration_index, + event.mon.tel[tel_id], + new_ped=new_ped_event, new_calib=new_calib_event) + + # calibrate and extract image from event r1_dl1_calibrator(event) + # update the calibration index in the dl1 event container + dl1_container.calibration_id = calibration_index.calibration_id + + + # Temporal volume reducer for lstchain - dl1 level must be filled and dl0 will be overwritten. # When the last version of the method is implemented, vol. reduction will be done at dl0 check_and_apply_volume_reduction(event, config) @@ -324,6 +366,7 @@ def r0_to_dl1(input_filename = get_dataset_path('gamma_test_large.simtel.gz'), tel_name = str(event.inst.subarray.tel[telescope_id])[4:] tel_name = tel_name.replace('-003', '') + if custom_calibration: lst_calibration(event, telescope_id) @@ -418,8 +461,14 @@ def r0_to_dl1(input_filename = get_dataset_path('gamma_test_large.simtel.gz'), # Until the TIB trigger_type is fully reliable, we also add # the ucts_trigger_type to the data - extra_im.ucts_trigger_type = event.lst.tel[telescope_id].evt.ucts_trigger_type + dl1_container.ucts_trigger_type = event.lst.tel[telescope_id].evt.ucts_trigger_type + dl1_container.trigger_time = event.r0.tel[telescope_id].trigger_time + dl1_container.trigger_type = event.r0.tel[telescope_id].trigger_type + + # info not available in data + # dl1_container.num_trig_pix = calibrated_event.r0.tel[telescope_id].num_trig_pix + # dl1_container.trig_pix_id = calibrated_event.r0.tel[telescope_id].trig_pix_id # FIXME: no need to read telescope characteristics like foclen for every event! foclen = event.inst.subarray.tel[telescope_id].optics.equivalent_focal_length @@ -428,14 +477,11 @@ def r0_to_dl1(input_filename = get_dataset_path('gamma_test_large.simtel.gz'), length = np.rad2deg(np.arctan2(dl1_container.length, foclen)) dl1_container.width = width.value dl1_container.length = length.value - dl1_container.prefix = tel.prefix + # extra info for the image table extra_im.tel_id = telescope_id - extra_im.num_trig_pix = event.r0.tel[telescope_id].num_trig_pix - extra_im.trigger_time = event.r0.tel[telescope_id].trigger_time - extra_im.trigger_type = event.r0.tel[telescope_id].trigger_type - extra_im.trig_pix_id = event.r0.tel[telescope_id].trig_pix_id + extra_im.low_gain_mask = event.lst.tel[telescope_id].evt.pixel_status >> 2 & 1 for container in [extra_im, dl1_container, event.r0, tel]: add_global_metadata(container, metadata) @@ -445,7 +491,7 @@ def r0_to_dl1(input_filename = get_dataset_path('gamma_test_large.simtel.gz'), writer.write(table_name = f'telescope/image/{tel_name}', containers = [event.r0, tel, extra_im]) writer.write(table_name = f'telescope/parameters/{tel_name}', - containers = [dl1_container, extra_im]) + containers = [dl1_container]) # Muon ring analysis, for real data only (MC is done starting from DL1 files) From dd4306cae5f33117aa65d035603cca540cd526fa Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 20 Apr 2020 15:52:17 +0200 Subject: [PATCH 10/45] Add ff time corrction at calibrator level --- lstchain/calib/camera/calibrator.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/lstchain/calib/camera/calibrator.py b/lstchain/calib/camera/calibrator.py index 0f804e4f78..a0dff88a1c 100644 --- a/lstchain/calib/camera/calibrator.py +++ b/lstchain/calib/camera/calibrator.py @@ -131,8 +131,6 @@ def _initialize_correction(self): # read the pixel_status container table = '/tel_' + str(telid) + '/pixel_status' next(h5_table.read(table, self.mon_data.tel[telid].pixel_status)) - - except: self.log.error(f"Problem in reading calibration file {self.calibration_path}") @@ -153,8 +151,8 @@ def _calibrate_dl0(self, event, telid): event.mon.tel[telid].pedestal = self.mon_data.tel[telid].pedestal event.mon.tel[telid].pixel_status = self.mon_data.tel[telid].pixel_status - - # subtract the pedestal per sample (should we do it?) and multiply for the calibration coefficients + # + # subtract the pedestal per sample and multiply for the calibration coefficients # event.dl0.tel[telid].waveform = ( (event.r1.tel[telid].waveform - event.mon.tel[telid].calibration.pedestal_per_sample[:, :, np.newaxis]) @@ -179,19 +177,17 @@ def _calibrate_dl1(self, event, telid): # correct time with drs4 correction if available if self.time_corrector: - pulse_corr_array = self.time_corrector.get_corr_pulse(event, pulse_time) + pulse_time = self.time_corrector.get_corr_pulse(event, pulse_time) - # otherwise use the ff time correction (not drs4 corrected) - else: - pulse_corr_array = pulse_time + self.mon_data.tel[telid].calibration.time_correction + # add flat-fielding time correction + pulse_time_ff_corrected = pulse_time + self.mon_data.tel[telid].calibration.time_correction # perform the gain selection if the threshold is defined - if self.gain_threshold: waveforms, gain_mask = self.gain_selector(event.r1.tel[telid].waveform) event.dl1.tel[telid].image = charge[gain_mask, np.arange(charge.shape[1])] - event.dl1.tel[telid].pulse_time = pulse_corr_array[gain_mask, np.arange(pulse_corr_array.shape[1])] + event.dl1.tel[telid].pulse_time = pulse_time_ff_corrected[gain_mask, np.arange(pulse_time_ff_corrected.shape[1])] # remember the mask in the lst pixel_status array (this info is missing for the moment in the # r1 container). I follow the prescription given in the document @@ -212,6 +208,6 @@ def _calibrate_dl1(self, event, telid): # if threshold == None else: event.dl1.tel[telid].image = charge - event.dl1.tel[telid].pulse_time = pulse_corr_array + event.dl1.tel[telid].pulse_time = pulse_time_ff_corrected From ccdc50f393b2c76608ffb873a939c2ab61396e55 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 20 Apr 2020 15:54:04 +0200 Subject: [PATCH 11/45] Add drs4 charge time correction --- lstchain/calib/camera/flatfield.py | 27 ++++++++++++++++++++++++--- 1 file changed, 24 insertions(+), 3 deletions(-) diff --git a/lstchain/calib/camera/flatfield.py b/lstchain/calib/camera/flatfield.py index c8a12b0dfa..02815a2129 100644 --- a/lstchain/calib/camera/flatfield.py +++ b/lstchain/calib/camera/flatfield.py @@ -2,12 +2,12 @@ Factory for the estimation of the flat field coefficients """ - +import os import numpy as np from astropy import units as u from ctapipe.calib.camera.flatfield import FlatFieldCalculator -from ctapipe.core.traits import List - +from ctapipe.core.traits import List, Unicode +from lstchain.calib.camera.pulse_time_correction import PulseTimeCorrection __all__ = [ 'FlasherFlatFieldCalculator' @@ -43,6 +43,11 @@ class FlasherFlatFieldCalculator(FlatFieldCalculator): [-3, 3], help='Interval (number of std) of accepted charge standard deviation around camera median value' ).tag(config=True) + time_calibration_path = Unicode( + '', + allow_none = True, + help = 'Path to drs4 time calibration file' + ).tag(config = True) def __init__(self, **kwargs): """Calculates flat-field parameters from flasher data @@ -74,6 +79,15 @@ def __init__(self, **kwargs): self.arrival_times = None # arrival time per event in sample self.sample_masked_pixels = None # masked pixels per event in sample + # declare time calibrator if correction file exist + if os.path.exists(self.time_calibration_path): + self.time_corrector = PulseTimeCorrection( + calib_file_path = self.time_calibration_path) + else: + self.time_corrector = None + self.log.info(f"File {self.time_calibration_path} not found. No drs4 time corrections") + + def _extract_charge(self, event): """ Extract the charge and the time from a calibration event @@ -96,6 +110,10 @@ def _extract_charge(self, event): charge, peak_pos = self.extractor(waveforms) + # correct time with drs4 correction if available + if self.time_corrector: + peak_pos = self.time_corrector.get_corr_pulse(event, peak_pos) + return charge, peak_pos def calculate_relative_gain(self, event): @@ -141,6 +159,9 @@ def calculate_relative_gain(self, event): # the peak position (assumed as time for the moment) charge, arrival_time = self._extract_charge(event) + # correct pulse time with drs4 corrections + + self.collect_sample(charge, pixel_mask, arrival_time) sample_age = self.trigger_time - self.time_start From e9cc86bff21da3fe59ee1a03ec3d3bd2def14316 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 20 Apr 2020 15:55:09 +0200 Subject: [PATCH 12/45] Add drs4 charge time corrections --- lstchain/calib/camera/pedestals.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lstchain/calib/camera/pedestals.py b/lstchain/calib/camera/pedestals.py index c530323b18..00aa65dd37 100644 --- a/lstchain/calib/camera/pedestals.py +++ b/lstchain/calib/camera/pedestals.py @@ -6,7 +6,7 @@ import numpy as np from astropy import units as u from ctapipe.calib.camera.pedestals import PedestalCalculator -from ctapipe.core.traits import Int, Unicode, List +from ctapipe.core.traits import List __all__ = [ 'PedestalIntegrator' From dbcba0cdb8986197caa44931b4c123f88c6b057a Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 20 Apr 2020 15:56:17 +0200 Subject: [PATCH 13/45] Remove time flat-field correction at drs4 level --- lstchain/calib/camera/pulse_time_correction.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/lstchain/calib/camera/pulse_time_correction.py b/lstchain/calib/camera/pulse_time_correction.py index a168cff948..ef07e32730 100644 --- a/lstchain/calib/camera/pulse_time_correction.py +++ b/lstchain/calib/camera/pulse_time_correction.py @@ -147,7 +147,9 @@ def get_first_capacitor(self, event, nr): @njit() def get_corr_time_jit(first_cap, fan, fbn, n_harmonics, n_cap): - time = fan[0] / 2. + + #time = fan[0] / 2. #commented because time flat-fielding is performed in charge calibrator + time = 0 for n in prange(1, n_harmonics): time += fan[n] * np.cos((first_cap * n * 2 * np.pi) / n_cap) time += fbn[n] * np.sin((first_cap * n * 2 * np.pi) / n_cap) From fff004f0f926c2cee1d40677ceb7194798a81da9 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 20 Apr 2020 15:58:10 +0200 Subject: [PATCH 14/45] Add the possiblity to force calibration output --- .../calib/camera/calibration_calculator.py | 35 +++++++++++-------- 1 file changed, 21 insertions(+), 14 deletions(-) diff --git a/lstchain/calib/camera/calibration_calculator.py b/lstchain/calib/camera/calibration_calculator.py index 70bc6146cd..d03c7efbca 100644 --- a/lstchain/calib/camera/calibration_calculator.py +++ b/lstchain/calib/camera/calibration_calculator.py @@ -5,9 +5,10 @@ from abc import abstractmethod import numpy as np +import os from ctapipe.core import Component from ctapipe.core import traits -from ctapipe.core.traits import Int, Float, List +from ctapipe.core.traits import Unicode, Float, List from lstchain.calib.camera.flatfield import FlatFieldCalculator from lstchain.calib.camera.pedestals import PedestalCalculator @@ -15,7 +16,6 @@ __all__ = [ 'CalibrationCalculator', 'LSTCalibrationCalculator' - ] @@ -103,9 +103,6 @@ def __init__( self.log.debug(f"{self.flatfield}") - - - class LSTCalibrationCalculator(CalibrationCalculator): """ Calibration calculator for LST camera @@ -120,6 +117,9 @@ class LSTCalibrationCalculator(CalibrationCalculator): maximum_lg_charge_std Temporary cut on LG std against Lidar events till the calibox TIB do not work (default for filter 5.2) + + time_calibration_path: + Path with the drs4 time calibration corrections """ minimum_hg_charge_median = Float( @@ -150,9 +150,6 @@ def __init__(self, **kwargs): """ super().__init__(**kwargs) - - - def calculate_calibration_coefficients(self, event): """ Calculate calibration coefficients from flatfield and pedestal statistics @@ -160,6 +157,7 @@ def calculate_calibration_coefficients(self, event): Parameters ---------- + event: EventAndMonDataContainer """ @@ -168,15 +166,12 @@ def calculate_calibration_coefficients(self, event): status_data = event.mon.tel[self.tel_id].pixel_status calib_data = event.mon.tel[self.tel_id].calibration - # mask from pedestal and flat-field data monitoring_unusable_pixels = np.logical_or(status_data.pedestal_failing_pixels, status_data.flatfield_failing_pixels) - # calibration unusable pixels are an OR of all masks calib_data.unusable_pixels = np.logical_or(monitoring_unusable_pixels, status_data.hardware_failing_pixels) - # Extract calibration coefficients with F-factor method # Assume fix F2 factor, F2=1+Var(gain)/Mean(Gain)**2 must be known from elsewhere F2 = 1.2 @@ -222,8 +217,8 @@ def process_interleaved(self, event): Parameters ---------- """ - new_calibration = False new_ped = False + new_ff = False # if pedestal event if event.r1.tel[self.tel_id].trigger_type == 32: @@ -245,7 +240,19 @@ def process_interleaved(self, event): # if new ff, calculate new calibration coefficients if new_ff: self.calculate_calibration_coefficients(event) - new_calibration=True - return new_calibration, new_ped + + return new_ped, new_ff + + def force_interleaved_results(self, event): + """ + Force output + + """ + #store results + self.pedestal.store_results(event) + self.flatfield.store_results(event) + + # calculates calibration values + self.calculate_calibration_coefficients(event) From b34bc2259468e06f39e7091c27bc565d85959c47 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 20 Apr 2020 16:01:44 +0200 Subject: [PATCH 15/45] Force interleaved calibraiton output at the end of file --- lstchain/reco/r0_to_dl1.py | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py index 644f024a1c..9fc57479d4 100644 --- a/lstchain/reco/r0_to_dl1.py +++ b/lstchain/reco/r0_to_dl1.py @@ -243,10 +243,14 @@ def r0_to_dl1(input_filename = get_dataset_path('gamma_test_large.simtel.gz'), config = Config(config), allowed_tels = [1],) + # Component to process interleaved pedestal and flat-fields + calib_config = Config(config[config['calibration_product']]) + calib_config.merge({"time_calibration_path": time_calibration_path}) + calibration_calculator = CalibrationCalculator.from_name( config['calibration_product'], - config=Config(config[config['calibration_product']]) + config=calib_config ) @@ -318,30 +322,39 @@ def r0_to_dl1(input_filename = get_dataset_path('gamma_test_large.simtel.gz'), else: if i==0: - # iniitalize the telescope + # initialize the telescope # FIXME? LST calibrator is only for one telescope - # it should be inside the telescope loop + # it should be inside the telescope loop (?) tel_id = calibration_calculator.tel_id - # write the first calibration event (initialized from h5 file) + # write the first calibration event (initialized from calibration h5 file) write_calibration_data(writer, calibration_index, r1_dl1_calibrator.mon_data.tel[tel_id], - new_ped=True, new_calib=True) + new_ped=True, new_ff=True) # drs4 calibrations r0_r1_calibrator.calibrate(event) - # process interleaved events (pedestals, ff, calibration) - new_calib_event, new_ped_event = calibration_calculator.process_interleaved(event) + + if i < source.max_events-1: + # process interleaved events (pedestals, ff, calibration) + new_ped_event, new_ff_event = calibration_calculator.process_interleaved(event) + + else: + # if end of file calculate results anyway + calibration_calculator.force_interleaved_results(event) + new_ped_event = True + new_ff_event = True + # write monitoring containers if updated - if new_ped_event or new_calib_event: + if new_ped_event or new_ff_event: write_calibration_data(writer, calibration_index, event.mon.tel[tel_id], - new_ped=new_ped_event, new_calib=new_calib_event) + new_ped=new_ped_event, new_ff=new_ff_event) # calibrate and extract image from event r1_dl1_calibrator(event) From 60403d072affc47314ebf0c247b8b0aa0ec81109 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 20 Apr 2020 16:03:21 +0200 Subject: [PATCH 16/45] First produce the time calibration file, which is now necessary for the charge calibration file --- .../onsite/onsite_create_calibration_file.py | 53 ++++++++++--------- 1 file changed, 27 insertions(+), 26 deletions(-) diff --git a/lstchain/scripts/onsite/onsite_create_calibration_file.py b/lstchain/scripts/onsite/onsite_create_calibration_file.py index 51c39aa3b4..d515a9c510 100755 --- a/lstchain/scripts/onsite/onsite_create_calibration_file.py +++ b/lstchain/scripts/onsite/onsite_create_calibration_file.py @@ -77,12 +77,6 @@ def main(): print(f">>> Error: The pedestal file {pedestal_file} do not exist.\n Exit") exit(0) - - # - # produce ff calibration file - # - - # define charge file names output_file = f"{output_dir}/calibration.Run{run}.0000.hdf5" log_file = f"{output_dir}/log/calibration.Run{run}.0000.log" @@ -103,22 +97,6 @@ def main(): exit(1) print(f"\n--> Config file {config_file}") - if ff_calibration is 'yes': - # run lstchain script - cmd = f"lstchain_create_calibration_file " \ - f"--input_file={input_file} --output_file={output_file} --pedestal_file={pedestal_file} " \ - f"--FlatFieldCalculator.sample_size={stat_events} --PedestalCalculator.sample_size={stat_events} " \ - f"--EventSource.max_events={max_events} --config={config_file} > {log_file} 2>&1" - - print("\n--> RUNNING...") - os.system(cmd) - - # plot and save some results - plot_file=f"{output_dir}/log/calibration.Run{run}.0000.pedestal.Run{ped_run}.0000.pdf" - print(f"\n--> PRODUCING PLOTS in {plot_file} ...") - calib.read_file(output_file,tel_id) - calib.plot_all(calib.ped_data, calib.ff_data, calib.calib_data, run, plot_file) - # # produce drs4 time calibration file # @@ -133,20 +111,43 @@ def main(): else: # otherwise perform a link to the default time calibration file print(f"\n--> PRODUCING LINK TO DEFAULT TIME CALIBRATION (run {default_time_run})") - file_list = sorted(Path(f"{base_dir}/calibration/").rglob(f'*/{prod_id}/time_calibration.Run{default_time_run}*')) - + file_list = sorted( + Path(f"{base_dir}/calibration/").rglob(f'*/{prod_id}/time_calibration.Run{default_time_run}*')) if len(file_list) == 0: print(f">>> Error: time calibration file for run {default_time_run} not found\n") raise NameError() else: time_calibration_file = file_list[0] - input_dir, name = os.path.split(os.path.abspath(time_calibration_file )) - cmd=f"ln -sf {time_calibration_file} {time_file}" + input_dir, name = os.path.split(os.path.abspath(time_calibration_file)) + cmd = f"ln -sf {time_calibration_file} {time_file}" os.system(cmd) print(f"\n--> Time calibration file: {time_file}") + + # + # produce ff calibration file + # + if ff_calibration is 'yes': + # run lstchain script + cmd = f"lstchain_create_calibration_file " \ + f"--input_file={input_file} --output_file={output_file} --pedestal_file={pedestal_file} " \ + f"--FlatFieldCalculator.time_calibration_path={time_file} --FlatFieldCalculator.sample_size={stat_events} "\ + f"--PedestalCalculator.sample_size={stat_events} " \ + f"--EventSource.max_events={max_events} --config={config_file} > {log_file} 2>&1" + + print("\n--> RUNNING...") + os.system(cmd) + + # plot and save some results + plot_file=f"{output_dir}/log/calibration.Run{run}.0000.pedestal.Run{ped_run}.0000.pdf" + print(f"\n--> PRODUCING PLOTS in {plot_file} ...") + calib.read_file(output_file,tel_id) + calib.plot_all(calib.ped_data, calib.ff_data, calib.calib_data, run, plot_file) + + + print("\n--> END") except Exception as e: From 874887996b68126e2557854baa7326094208abb9 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 20 Apr 2020 16:05:17 +0200 Subject: [PATCH 17/45] Container initalization changed --- lstchain/io/lstcontainers.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lstchain/io/lstcontainers.py b/lstchain/io/lstcontainers.py index 84f89f06f9..3e141c8815 100644 --- a/lstchain/io/lstcontainers.py +++ b/lstchain/io/lstcontainers.py @@ -247,6 +247,6 @@ class DL1MonitoringEventIndexContainer(Container): Container with the calibration coefficients """ tel_id = Field(1, 'Index of telescope') - calibration_id = Field(0, 'Index of calibration event for DL1 file') - pedestal_id = Field(0, 'Index of pedestal event for DL1 file') - flatfield_id = Field(0, 'Index of flat-field event for DL1 file') \ No newline at end of file + calibration_id = Field(-1, 'Index of calibration event for DL1 file') + pedestal_id = Field(-1, 'Index of pedestal event for DL1 file') + flatfield_id = Field(-1, 'Index of flat-field event for DL1 file') \ No newline at end of file From 6a2388a9bf0b062c871da67bc19569f50665206f Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 20 Apr 2020 16:06:15 +0200 Subject: [PATCH 18/45] Minor changes --- lstchain/io/io.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/lstchain/io/io.py b/lstchain/io/io.py index be62baaa6c..a365de863c 100644 --- a/lstchain/io/io.py +++ b/lstchain/io/io.py @@ -668,7 +668,7 @@ def recursive_copy_node(src_file, dir_file, path): recursive_path = os.path.join(recursive_path, p) -def write_calibration_data(writer,mon_index, mon_event, new_ped=False, new_calib=False): +def write_calibration_data(writer,mon_index, mon_event, new_ped=False, new_ff=False): mon_event.pedestal.prefix = '' mon_event.flatfield.prefix = '' mon_event.calibration.prefix = '' @@ -676,14 +676,16 @@ def write_calibration_data(writer,mon_index, mon_event, new_ped=False, new_calib if new_ped: # write ped container + mon_index.pedestal_id += 1 writer.write( table_name=f'telescope/monitoring/pedestal', containers=[mon_index, mon_event.pedestal] ) - # update index - mon_index.pedestal_id += 1 - if new_calib: + if new_ff: + mon_index.flatfield_id += 1 + mon_index.calibration_id += 1 + # write calibration container writer.write( table_name="telescope/monitoring/flatfield", @@ -695,6 +697,3 @@ def write_calibration_data(writer,mon_index, mon_event, new_ped=False, new_calib table_name="telescope/monitoring/calibration", containers=[mon_index, mon_event.calibration] ) - # update index - mon_index.flatfield_id += 1 - mon_index.calibration_id += 1 \ No newline at end of file From 47716281f711dd5bf888f32405c04331167e2573 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 20 Apr 2020 16:07:00 +0200 Subject: [PATCH 19/45] Ask high statistic for pedestals and flatfields, so only one calibration event will be created per sub-run and all events wll be calibrated with the beginning of night run --- lstchain/data/lstchain_standard_config.json | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lstchain/data/lstchain_standard_config.json b/lstchain/data/lstchain_standard_config.json index 4243500e16..2a3bfbf0c1 100644 --- a/lstchain/data/lstchain_standard_config.json +++ b/lstchain/data/lstchain_standard_config.json @@ -125,7 +125,7 @@ "flatfield_product": "FlasherFlatFieldCalculator", "pedestal_product": "PedestalIntegrator", "PedestalCalculator":{ - "sample_size": 100, + "sample_size": 10000, "sample_duration":100000, "tel_id":1, "charge_median_cut_outliers": [-10,10], @@ -133,7 +133,7 @@ "charge_product":"FixedWindowSum" }, "FlatFieldCalculator":{ - "sample_size": 100, + "sample_size": 10000, "sample_duration":100000, "tel_id":1, "charge_product":"LocalPeakWindowSum", From eb449b5803fe0ee3075965e7da8507eeffa934db Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 20 Apr 2020 16:10:18 +0200 Subject: [PATCH 20/45] Eliminate the "fits" postfix in the h5 file output name --- lstchain/scripts/lstchain_data_r0_to_dl1.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/lstchain/scripts/lstchain_data_r0_to_dl1.py b/lstchain/scripts/lstchain_data_r0_to_dl1.py index b2f4778d24..1f3b491a2e 100644 --- a/lstchain/scripts/lstchain_data_r0_to_dl1.py +++ b/lstchain/scripts/lstchain_data_r0_to_dl1.py @@ -94,7 +94,12 @@ def main(): os.makedirs(args.outdir, exist_ok=True) r0_to_dl1.allowed_tels = {1, 2, 3, 4} - output_filename = args.outdir + '/dl1_' + os.path.basename(args.infile).rsplit('.', 1)[0] + '.h5' + + # if simu + output_filename = args.outdir + '/dl1_' + os.path.basename(args.infile).replace(".gz",".h5") + + # if data + output_filename = args.outdir + '/dl1_' + os.path.basename(args.infile).replace(".fits.fz",".h5") config = {} if args.config_file is not None: From 814ee9dc2534603b4c552dca3aaec8e8b4bc77a3 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 20 Apr 2020 16:19:12 +0200 Subject: [PATCH 21/45] Improve prints --- .../onsite/onsite_create_calibration_file.py | 43 ++++++++++--------- 1 file changed, 23 insertions(+), 20 deletions(-) diff --git a/lstchain/scripts/onsite/onsite_create_calibration_file.py b/lstchain/scripts/onsite/onsite_create_calibration_file.py index d515a9c510..807394ad11 100755 --- a/lstchain/scripts/onsite/onsite_create_calibration_file.py +++ b/lstchain/scripts/onsite/onsite_create_calibration_file.py @@ -77,31 +77,12 @@ def main(): print(f">>> Error: The pedestal file {pedestal_file} do not exist.\n Exit") exit(0) - # define charge file names - output_file = f"{output_dir}/calibration.Run{run}.0000.hdf5" - log_file = f"{output_dir}/log/calibration.Run{run}.0000.log" - print(f"\n--> Output file {output_file}") - if os.path.exists(output_file) and ff_calibration is 'yes': - if query_yes_no(">>> Output file exists already. Do you want to remove it?"): - os.remove(output_file) - else: - print(f"\n--> Stop") - exit(1) - - print(f"\n--> Log file {log_file}") - - # define config file - config_file = os.path.join(os.path.dirname(__file__), "../../data/onsite_camera_calibration_param.json") - if not os.path.exists(config_file): - print(f">>> Config file {config_file} do not exists. \n Exit ") - exit(1) - print(f"\n--> Config file {config_file}") # # produce drs4 time calibration file # time_file = f"{output_dir}/time_calibration.Run{run}.0000.hdf5" - + print(f"\n***** PRODUCE TIME CALIBRATION FILE ***** ") if default_time_run is 0: print(f"\n--> PRODUCING TIME CALIBRATION in {time_file} ...") cmd = f"lstchain_data_create_time_calibration_file --input_file {input_file} " \ @@ -125,6 +106,28 @@ def main(): print(f"\n--> Time calibration file: {time_file}") + # define charge file names + print(f"\n***** PRODUCE CHARGE CALIBRATION FILE ***** ") + output_file = f"{output_dir}/calibration.Run{run}.0000.hdf5" + log_file = f"{output_dir}/log/calibration.Run{run}.0000.log" + print(f"\n--> Output file {output_file}") + if os.path.exists(output_file) and ff_calibration is 'yes': + if query_yes_no(">>> Output file exists already. Do you want to remove it?"): + os.remove(output_file) + else: + print(f"\n--> Stop") + exit(1) + + print(f"\n--> Log file {log_file}") + + # define config file + config_file = os.path.join(os.path.dirname(__file__), "../../data/onsite_camera_calibration_param.json") + if not os.path.exists(config_file): + print(f">>> Config file {config_file} do not exists. \n Exit ") + exit(1) + print(f"\n--> Config file {config_file}") + + # # produce ff calibration file From 4bf3fdb5b494b1f51339dde322c40f362c7cc09f Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Mon, 20 Apr 2020 16:23:36 +0200 Subject: [PATCH 22/45] Changed name to script --- .../lstchain_data_create_calibration_file.py | 273 ------------------ 1 file changed, 273 deletions(-) delete mode 100644 lstchain/tools/lstchain_data_create_calibration_file.py diff --git a/lstchain/tools/lstchain_data_create_calibration_file.py b/lstchain/tools/lstchain_data_create_calibration_file.py deleted file mode 100644 index 0f306e27b8..0000000000 --- a/lstchain/tools/lstchain_data_create_calibration_file.py +++ /dev/null @@ -1,273 +0,0 @@ -""" -Extract flat field coefficients from flasher data files. -""" -import numpy as np -from traitlets import Dict, List, Unicode, Float, Bool - - -from ctapipe.core import Provenance, traits -from ctapipe.io import HDF5TableWriter -from ctapipe.core import Tool -from ctapipe.io import EventSource -from ctapipe.io.containers import PixelStatusContainer -from lstchain.calib.camera.calibration_calculator import CalibrationCalculator -from lstchain.calib.camera import CameraR0Calibrator - -__all__ = [ - 'CalibrationHDF5Writer' -] - - -class CalibrationHDF5Writer(Tool): - """ - Tool that generates a HDF5 file with camera calibration coefficients. - This is just an example on how the monitoring containers can be filled using - the calibration Components in calib/camera. - This example is based on an input file with pedestal and flat-field events - - For getting help run: - python calc_camera_calibration.py --help - - """ - - name = "CalibrationHDF5Writer" - description = "Generate a HDF5 file with camera calibration coefficients" - - - one_event = Bool( - False, - help='Stop after first calibration event' - ).tag(config=True) - - output_file = Unicode( - 'calibration.hdf5', - help='Name of the output file' - ).tag(config=True) - - log_file = Unicode( - 'None', - help='Name of the log file' - ).tag(config=True) - - calibration_product = traits.enum_trait( - CalibrationCalculator, - default='LSTCalibrationCalculator' - ) - - r0calibrator_product =traits.enum_trait( - CameraR0Calibrator, - default='NullR0Calibrator' - ) - - aliases = Dict(dict( - input_file='EventSource.input_url', - output_file='CalibrationHDF5Writer.output_file', - log_file='CalibrationHDF5Writer.log_file', - max_events='EventSource.max_events', - pedestal_file= 'LSTR0Corrections.pedestal_path', - calibration_product='CalibrationHDF5Writer.calibration_product', - r0calibrator_product='CalibrationHDF5Writer.r0calibrator_product', - )) - - classes = List([EventSource, - CalibrationCalculator - ] - + traits.classes_with_traits(CameraR0Calibrator) - + traits.classes_with_traits(CalibrationCalculator) - - ) - - def __init__(self, **kwargs): - super().__init__(**kwargs) - """ - Tool that generates a HDF5 file with camera calibration coefficients. - Input file must contain interleaved pedestal and flat-field events - - For getting help run: - python calc_camera_calibration.py --help - - """ - self.eventsource = None - self.processor = None - self.container = None - self.writer = None - self.r0calibrator = None - self.tot_events = 0 - self.simulation = False - - def setup(self): - kwargs = dict(parent=self) - self.log.debug(f"Open file") - self.eventsource = EventSource.from_config(**kwargs) - - # if data remember how many event in the files - if "LSTEventSource" in str(type(self.eventsource)): - - self.tot_events = len(self.eventsource.multi_file) - self.log.debug(f"Input file has file {self.tot_events} events") - else: - self.tot_events = self.eventsource.max_events - self.simulation = True - - self.processor = CalibrationCalculator.from_name( - self.calibration_product, - **kwargs - ) - - if self.r0calibrator_product: - self.r0calibrator = CameraR0Calibrator.from_name( - self.r0calibrator_product, - **kwargs - ) - - - group_name = 'tel_' + str(self.processor.tel_id) - - self.log.debug(f"Open output file {self.output_file}") - - self.writer = HDF5TableWriter( - filename=self.output_file, group_name=group_name, overwrite=True - ) - - def start(self): - '''Calibration coefficient calculator''' - - tel_id = self.processor.tel_id - new_ped = False - new_ff = False - end_of_file = False - - try: - self.log.debug(f"Start loop") - for count, event in enumerate(self.eventsource): - - if count % 100 == 0: - self.log.debug(f"Event {count}") - - # if last event write results - if count == self.tot_events-1 or count == self.eventsource.max_events-1: - self.log.debug(f"Last event, count = {count}") - end_of_file = True - - # save the config, to be retrieved as data.meta['config'] - if count == 0: - - if self.simulation: - initialize_pixel_status(event.mon.tel[tel_id],event.r1.tel[tel_id].waveform.shape) - - ped_data = event.mon.tel[tel_id].pedestal - ped_data.meta['config'] = self.config - - ff_data = event.mon.tel[tel_id].flatfield - ff_data.meta['config'] = self.config - - status_data = event.mon.tel[tel_id].pixel_status - status_data.meta['config'] = self.config - - calib_data = event.mon.tel[tel_id].calibration - calib_data.meta['config'] = self.config - - # correct for low level calibration - self.r0calibrator.calibrate(event) - - # reject event without trigger type - if event.r1.tel[tel_id].trigger_type == -1: - continue - - # if pedestal event - - if event.r1.tel[tel_id].trigger_type == 32 or ( - self.simulation and - np.median(np.sum(event.r1.tel[tel_id].waveform[0], axis=1)) - < self.processor.minimum_hg_charge_median): - - new_ped = self.processor.pedestal.calculate_pedestals(event) - - - # if flat-field event: no calibration TIB for the moment, - # use a cut on the charge for ff events and on std for rejecting Magic Lidar events - elif event.r1.tel[tel_id].trigger_type == 4 or ( - np.median(np.sum(event.r1.tel[tel_id].waveform[0], axis=1)) - > self.processor.minimum_hg_charge_median - and np.std(np.sum(event.r1.tel[tel_id].waveform[1], axis=1)) - < self.processor.maximum_lg_charge_std): - - new_ff = self.processor.flatfield.calculate_relative_gain(event) - - # write pedestal results when enough statistics or end of file - if new_ped or end_of_file: - - # update the monitoring container with the last statistics - if end_of_file: - self.processor.pedestal.store_results(event) - - # write the event - self.log.debug(f"Write pedestal data at event n. {count+1}, id {event.r0.event_id} " - f"stat = {ped_data.n_events} events") - - # write on file - self.writer.write('pedestal', ped_data) - - new_ped = False - - # write flatfield results when enough statistics (also for pedestals) or end of file - if (new_ff and ped_data.n_events > 0) or end_of_file: - - # update the monitoring container with the last statistics - if end_of_file: - self.processor.flatfield.store_results(event) - - self.log.debug(f"Write flatfield data at event n. {count+1}, id {event.r0.event_id} " - f"stat = {ff_data.n_events} events") - - # write on file - self.writer.write('flatfield', ff_data) - - new_ff = False - - # Then, calculate calibration coefficients - self.processor.calculate_calibration_coefficients(event) - - # write calib and pixel status - self.log.debug(f"Write pixel_status data") - self.writer.write('pixel_status',status_data) - - self.log.debug(f"Write calibration data") - self.writer.write('calibration', calib_data) - if self.one_event: - break - - #self.writer.write('mon', event.mon.tel[tel_id]) - except ValueError as e: - self.log.error(e) - - def finish(self): - Provenance().add_output_file( - self.output_file, - role='mon.tel.calibration' - ) - self.writer.close() - -def initialize_pixel_status(mon_camera_container,shape): - """ - Initialize the pixel status container in the case of - simulation events (this should be done in the event source, but - added here for the moment) - """ - # initialize the container - status_container = PixelStatusContainer() - status_container.hardware_failing_pixels = np.zeros((shape[0],shape[1]), dtype=bool) - status_container.pedestal_failing_pixels = np.zeros((shape[0],shape[1]), dtype=bool) - status_container.flatfield_failing_pixels = np.zeros((shape[0],shape[1]), dtype=bool) - - mon_camera_container.pixel_status = status_container - - -def main(): - exe = CalibrationHDF5Writer() - - exe.run() - - -if __name__ == '__main__': - main() From 7461bff5a1b7b86c5c9a15931c7d4e0a11355273 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Tue, 21 Apr 2020 09:24:49 +0200 Subject: [PATCH 23/45] Changed logic to assure that last interleaved statistics is written on file --- lstchain/reco/r0_to_dl1.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py index f53f8a8ddd..31c3b886a1 100644 --- a/lstchain/reco/r0_to_dl1.py +++ b/lstchain/reco/r0_to_dl1.py @@ -356,16 +356,9 @@ def r0_to_dl1( r0_r1_calibrator.calibrate(event) - if i < source.max_events-1: - # process interleaved events (pedestals, ff, calibration) - new_ped_event, new_ff_event = calibration_calculator.process_interleaved(event) - - else: - # if end of file calculate results anyway - calibration_calculator.force_interleaved_results(event) - new_ped_event = True - new_ff_event = True + # process interleaved events (pedestals, ff, calibration) + new_ped_event, new_ff_event = calibration_calculator.process_interleaved(event) # write monitoring containers if updated if new_ped_event or new_ff_event: @@ -590,6 +583,17 @@ def r0_to_dl1( containers = [event.mc.tel[telescope_id], extra_im] ) + + + # at the end of event loop force calculation if interleaved statistics + calibration_calculator.force_interleaved_results(event) + # write monitoring events + write_calibration_data(writer, + calibration_index, + event.mon.tel[tel_id], + new_ped=True, new_ff=True) + + if is_simu: ### Reconstruct source position from disp for all events and write the result in the output file for tel_name in ['LST_LSTCam']: From 4890cab04e5a5a40926eab60a0eeb06c714f7d1e Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Tue, 21 Apr 2020 09:42:09 +0200 Subject: [PATCH 24/45] Minor changes --- lstchain/reco/r0_to_dl1.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py index af525e3092..4b3dfbf497 100644 --- a/lstchain/reco/r0_to_dl1.py +++ b/lstchain/reco/r0_to_dl1.py @@ -358,8 +358,6 @@ def r0_to_dl1( # drs4 calibrations r0_r1_calibrator.calibrate(event) - - # process interleaved events (pedestals, ff, calibration) new_ped_event, new_ff_event = calibration_calculator.process_interleaved(event) From cc92303cfa7b0c6952102399fb3708433ce8d1d2 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Tue, 21 Apr 2020 10:01:45 +0200 Subject: [PATCH 25/45] Correct index update and merge problem --- lstchain/io/io.py | 14 ++++++++++---- lstchain/reco/r0_to_dl1.py | 4 +--- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/lstchain/io/io.py b/lstchain/io/io.py index 5b4f877ce3..219348efe5 100644 --- a/lstchain/io/io.py +++ b/lstchain/io/io.py @@ -677,18 +677,24 @@ def write_calibration_data(writer,mon_index, mon_event, new_ped=False, new_ff=Fa mon_event.calibration.prefix = '' mon_index.prefix = '' + + # update index if new_ped: - # write ped container mon_index.pedestal_id += 1 + + if new_ff: + mon_index.flatfield_id += 1 + mon_index.calibration_id += 1 + + + if new_ped: + # write ped container writer.write( table_name=f'telescope/monitoring/pedestal', containers=[mon_index, mon_event.pedestal] ) if new_ff: - mon_index.flatfield_id += 1 - mon_index.calibration_id += 1 - # write calibration container writer.write( table_name="telescope/monitoring/flatfield", diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py index 1981250c15..69a5bb6d22 100644 --- a/lstchain/reco/r0_to_dl1.py +++ b/lstchain/reco/r0_to_dl1.py @@ -330,10 +330,8 @@ def r0_to_dl1( writer._h5file.filters = filters print("USING FILTERS: ", writer._h5file.filters) - for i, event in enumerate(chain([first_event], event_iter)): - ### EVENT LOOP ### - for i, event in enumerate(source): + for i, event in enumerate(chain([first_event], event_iter)): if i % 100 == 0: print(i) From c633f05994b9e1578bf5077660b1a6f69550ad4d Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Tue, 21 Apr 2020 17:17:17 +0200 Subject: [PATCH 26/45] Add trailet for excess noise factor. Improve function name --- .../calib/camera/calibration_calculator.py | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/lstchain/calib/camera/calibration_calculator.py b/lstchain/calib/camera/calibration_calculator.py index d03c7efbca..42eab658f9 100644 --- a/lstchain/calib/camera/calibration_calculator.py +++ b/lstchain/calib/camera/calibration_calculator.py @@ -2,13 +2,10 @@ Component for the estimation of the calibration coefficients events """ - -from abc import abstractmethod import numpy as np -import os from ctapipe.core import Component from ctapipe.core import traits -from ctapipe.core.traits import Unicode, Float, List +from ctapipe.core.traits import Float, List from lstchain.calib.camera.flatfield import FlatFieldCalculator from lstchain.calib.camera.pedestals import PedestalCalculator @@ -38,6 +35,10 @@ class CalibrationCalculator(Component): kwargs """ + squared_excess_noise_factor = Float( + 1.2, + help='Excess noise factor squared: 1+ Var(gain)/Mean(Gain)**2' + ).tag(config=True) pedestal_product = traits.enum_trait( PedestalCalculator, @@ -173,14 +174,13 @@ def calculate_calibration_coefficients(self, event): calib_data.unusable_pixels = np.logical_or(monitoring_unusable_pixels, status_data.hardware_failing_pixels) # Extract calibration coefficients with F-factor method - # Assume fix F2 factor, F2=1+Var(gain)/Mean(Gain)**2 must be known from elsewhere - F2 = 1.2 + # Assume fixed excess noise factor must be known from elsewhere # calculate photon-electrons - n_pe = F2 * (ff_data.charge_median - ped_data.charge_median) ** 2 / ( + n_pe = self.squared_excess_noise_factor * (ff_data.charge_median - ped_data.charge_median) ** 2 / ( ff_data.charge_std ** 2 - ped_data.charge_std ** 2) - # fill WaveformCalibrationContainer (this is an example) + # fill WaveformCalibrationContainer calib_data.time = ff_data.sample_time calib_data.time_range = ff_data.sample_time_range calib_data.n_pe = n_pe @@ -244,9 +244,9 @@ def process_interleaved(self, event): return new_ped, new_ff - def force_interleaved_results(self, event): + def output_interleaved_results(self, event): """ - Force output + Output interleaved results on request """ #store results From e001a230679f890cfd2262c27663af65eafd4f69 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Tue, 21 Apr 2020 17:18:36 +0200 Subject: [PATCH 27/45] Require time calibation file as input --- lstchain/calib/camera/calibrator.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/lstchain/calib/camera/calibrator.py b/lstchain/calib/camera/calibrator.py index a0dff88a1c..9da4472ee1 100644 --- a/lstchain/calib/camera/calibrator.py +++ b/lstchain/calib/camera/calibrator.py @@ -95,8 +95,7 @@ def __init__(self, **kwargs): self.time_corrector = PulseTimeCorrection( calib_file_path = self.time_calibration_path) else: - self.time_corrector = None - self.log.info(f"File {self.time_calibration_path} not found. No drs4 time corrections") + raise IOError(f"Time calibration file {self.time_calibration_path} not found!") # calibration data container self.mon_data = MonitoringContainer() From 3544032ca1cb220a725a4e68038aa4620996e323 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Tue, 21 Apr 2020 17:19:45 +0200 Subject: [PATCH 28/45] Require time camibration file as input --- lstchain/calib/camera/flatfield.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/lstchain/calib/camera/flatfield.py b/lstchain/calib/camera/flatfield.py index 02815a2129..9dd8ac4d62 100644 --- a/lstchain/calib/camera/flatfield.py +++ b/lstchain/calib/camera/flatfield.py @@ -6,7 +6,7 @@ import numpy as np from astropy import units as u from ctapipe.calib.camera.flatfield import FlatFieldCalculator -from ctapipe.core.traits import List, Unicode +from ctapipe.core.traits import List, Unicode, Path from lstchain.calib.camera.pulse_time_correction import PulseTimeCorrection __all__ = [ @@ -45,7 +45,6 @@ class FlasherFlatFieldCalculator(FlatFieldCalculator): ).tag(config=True) time_calibration_path = Unicode( '', - allow_none = True, help = 'Path to drs4 time calibration file' ).tag(config = True) @@ -84,8 +83,7 @@ def __init__(self, **kwargs): self.time_corrector = PulseTimeCorrection( calib_file_path = self.time_calibration_path) else: - self.time_corrector = None - self.log.info(f"File {self.time_calibration_path} not found. No drs4 time corrections") + raise IOError(f"Time calibration file {self.time_calibration_path} not found!") def _extract_charge(self, event): From 19fa10304d14b17c6f1ec516a75b1fd6f2715ced Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Tue, 21 Apr 2020 17:21:01 +0200 Subject: [PATCH 29/45] Change function name. Correct trailet initialization --- lstchain/reco/r0_to_dl1.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py index 69a5bb6d22..7bab4caa99 100644 --- a/lstchain/reco/r0_to_dl1.py +++ b/lstchain/reco/r0_to_dl1.py @@ -261,8 +261,8 @@ def r0_to_dl1( # Component to process interleaved pedestal and flat-fields calib_config = Config(config[config['calibration_product']]) - calib_config.merge({"time_calibration_path": time_calibration_path}) - + # time calibration path to flatfield trailet + calib_config.FlatFieldCalculator.time_calibration_path = time_calibration_path calibration_calculator = CalibrationCalculator.from_name( config['calibration_product'], config=calib_config @@ -592,8 +592,8 @@ def r0_to_dl1( - # at the end of event loop force calculation if interleaved statistics - calibration_calculator.force_interleaved_results(event) + # at the end of event loop ask calculation of remaining interleaved statistics + calibration_calculator.output_interleaved_results(event) # write monitoring events write_calibration_data(writer, calibration_index, From c2e6267ba1e3bd852e4527ad4aee9d3142d6013f Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 22 Apr 2020 10:06:16 +0200 Subject: [PATCH 30/45] Add LSTEventType class to recognize the event type --- .../calib/camera/calibration_calculator.py | 6 +-- lstchain/calib/camera/calibrator.py | 1 + lstchain/io/lstcontainers.py | 43 ++++++++++++++++++- .../tools/lstchain_create_calibration_file.py | 7 +-- 4 files changed, 50 insertions(+), 7 deletions(-) diff --git a/lstchain/calib/camera/calibration_calculator.py b/lstchain/calib/camera/calibration_calculator.py index 42eab658f9..c801252063 100644 --- a/lstchain/calib/camera/calibration_calculator.py +++ b/lstchain/calib/camera/calibration_calculator.py @@ -8,7 +8,7 @@ from ctapipe.core.traits import Float, List from lstchain.calib.camera.flatfield import FlatFieldCalculator from lstchain.calib.camera.pedestals import PedestalCalculator - +from lstchain.io.lstcontainers import LSTEventType __all__ = [ 'CalibrationCalculator', @@ -221,14 +221,14 @@ def process_interleaved(self, event): new_ff = False # if pedestal event - if event.r1.tel[self.tel_id].trigger_type == 32: + if LSTEventType.is_pedestal(event.r1.tel[self.tel_id].trigger_type): new_ped = self.pedestal.calculate_pedestals(event) # if flat-field event: no calibration TIB for the moment, # use a cut on the charge for ff events and on std for rejecting Magic Lidar events - elif event.r1.tel[self.tel_id].trigger_type == 4 or ( + elif LSTEventType.is_calibration(event.r1.tel[self.tel_id].trigger_type) or ( np.median(np.sum(event.r1.tel[self.tel_id].waveform[0], axis=1)) > self.minimum_hg_charge_median and np.std(np.sum(event.r1.tel[self.tel_id].waveform[1], axis=1)) diff --git a/lstchain/calib/camera/calibrator.py b/lstchain/calib/camera/calibrator.py index 9da4472ee1..a8e3de892b 100644 --- a/lstchain/calib/camera/calibrator.py +++ b/lstchain/calib/camera/calibrator.py @@ -9,6 +9,7 @@ from ctapipe.calib.camera import gainselection from lstchain.calib.camera.pulse_time_correction import PulseTimeCorrection + __all__ = ['LSTCameraCalibrator'] diff --git a/lstchain/io/lstcontainers.py b/lstchain/io/lstcontainers.py index badba13f56..4bcb7b7261 100644 --- a/lstchain/io/lstcontainers.py +++ b/lstchain/io/lstcontainers.py @@ -18,7 +18,8 @@ 'DispContainer', 'MetaData', 'ThrownEventsHistogram', - 'DL1MonitoringEventIndexContainer' + 'DL1MonitoringEventIndexContainer', + 'LSTEventType' ] @@ -266,3 +267,43 @@ class DL1MonitoringEventIndexContainer(Container): calibration_id = Field(-1, 'Index of calibration event for DL1 file') pedestal_id = Field(-1, 'Index of pedestal event for DL1 file') flatfield_id = Field(-1, 'Index of flat-field event for DL1 file') + +class LSTEventType: + """ + Class to recognize event type from trigger bits + bit 0: Mono + bit 1: stereo + bit 2: Calibration + bit 3: Single Phe + bit 4: Softrig(from the UCTS) + bit 5: Pedestal + bit 6: slow control + bit 7: busy + """ + + def is_mono(trigger_type): + return trigger_type >> 0 & 1 + + def is_stereo(trigger_type): + return trigger_type >> 1 & 1 + + def is_calibration(trigger_type): + return trigger_type >> 2 & 1 + + def is_single_pe(trigger_type): + return trigger_type >> 3 & 1 + + def is_soft_trig(trigger_type): + return trigger_type >> 4 & 1 + + def is_pedestal(trigger_type): + return trigger_type >> 5 & 1 + + def is_slow_control(trigger_type): + return trigger_type >> 6 & 1 + + def is_busy(trigger_type): + return trigger_type >> 7 & 1 + + def is_unknown(trigger_type): + return trigger_type == -1 diff --git a/lstchain/tools/lstchain_create_calibration_file.py b/lstchain/tools/lstchain_create_calibration_file.py index 0f306e27b8..b6d9b00134 100644 --- a/lstchain/tools/lstchain_create_calibration_file.py +++ b/lstchain/tools/lstchain_create_calibration_file.py @@ -12,6 +12,7 @@ from ctapipe.io.containers import PixelStatusContainer from lstchain.calib.camera.calibration_calculator import CalibrationCalculator from lstchain.calib.camera import CameraR0Calibrator +from lstchain.io.lstcontainers import LSTEventType __all__ = [ 'CalibrationHDF5Writer' @@ -171,12 +172,12 @@ def start(self): self.r0calibrator.calibrate(event) # reject event without trigger type - if event.r1.tel[tel_id].trigger_type == -1: + if LSTEventType.is_unknown(event.r1.tel[tel_id].trigger_type): continue # if pedestal event - if event.r1.tel[tel_id].trigger_type == 32 or ( + if LSTEventType.is_pedestal(event.r1.tel[tel_id].trigger_type) or ( self.simulation and np.median(np.sum(event.r1.tel[tel_id].waveform[0], axis=1)) < self.processor.minimum_hg_charge_median): @@ -186,7 +187,7 @@ def start(self): # if flat-field event: no calibration TIB for the moment, # use a cut on the charge for ff events and on std for rejecting Magic Lidar events - elif event.r1.tel[tel_id].trigger_type == 4 or ( + elif LSTEventType.is_calibration(event.r1.tel[tel_id].trigger_type) or ( np.median(np.sum(event.r1.tel[tel_id].waveform[0], axis=1)) > self.processor.minimum_hg_charge_median and np.std(np.sum(event.r1.tel[tel_id].waveform[1], axis=1)) From 0abd6131f357f56c0200e587b29b47bb5fb116b5 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 22 Apr 2020 10:42:48 +0200 Subject: [PATCH 31/45] Add "static" property to static class --- lstchain/io/lstcontainers.py | 1 + 1 file changed, 1 insertion(+) diff --git a/lstchain/io/lstcontainers.py b/lstchain/io/lstcontainers.py index 4bcb7b7261..92c5fddb42 100644 --- a/lstchain/io/lstcontainers.py +++ b/lstchain/io/lstcontainers.py @@ -268,6 +268,7 @@ class DL1MonitoringEventIndexContainer(Container): pedestal_id = Field(-1, 'Index of pedestal event for DL1 file') flatfield_id = Field(-1, 'Index of flat-field event for DL1 file') +@staticmethod class LSTEventType: """ Class to recognize event type from trigger bits From c5c5e3c2c34cf1d6580874691486eaa71ea862ad Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 22 Apr 2020 10:49:26 +0200 Subject: [PATCH 32/45] Add static property to static function --- lstchain/io/lstcontainers.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/lstchain/io/lstcontainers.py b/lstchain/io/lstcontainers.py index 92c5fddb42..1c2e33cefd 100644 --- a/lstchain/io/lstcontainers.py +++ b/lstchain/io/lstcontainers.py @@ -268,7 +268,7 @@ class DL1MonitoringEventIndexContainer(Container): pedestal_id = Field(-1, 'Index of pedestal event for DL1 file') flatfield_id = Field(-1, 'Index of flat-field event for DL1 file') -@staticmethod + class LSTEventType: """ Class to recognize event type from trigger bits @@ -282,29 +282,38 @@ class LSTEventType: bit 7: busy """ + @staticmethod def is_mono(trigger_type): return trigger_type >> 0 & 1 + @staticmethod def is_stereo(trigger_type): return trigger_type >> 1 & 1 + @staticmethod def is_calibration(trigger_type): return trigger_type >> 2 & 1 + @staticmethod def is_single_pe(trigger_type): return trigger_type >> 3 & 1 + @staticmethod def is_soft_trig(trigger_type): return trigger_type >> 4 & 1 + @staticmethod def is_pedestal(trigger_type): return trigger_type >> 5 & 1 + @staticmethod def is_slow_control(trigger_type): return trigger_type >> 6 & 1 + @staticmethod def is_busy(trigger_type): return trigger_type >> 7 & 1 + @staticmethod def is_unknown(trigger_type): return trigger_type == -1 From 39b9269d71d96dc2f42dd7c4c2ee07547fe85670 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 22 Apr 2020 13:18:58 +0200 Subject: [PATCH 33/45] Replace kargs with config and parent in the CalibrationCalculator Component Replace kargs with config and parent in the CalibrationCalculator Component --- .../calib/camera/calibration_calculator.py | 34 ++++--------------- lstchain/data/lstchain_standard_config.json | 4 +-- .../data/onsite_camera_calibration_param.json | 4 +-- lstchain/reco/r0_to_dl1.py | 8 +++-- .../onsite/onsite_create_calibration_file.py | 4 +-- .../tools/lstchain_create_calibration_file.py | 9 +++-- 6 files changed, 21 insertions(+), 42 deletions(-) diff --git a/lstchain/calib/camera/calibration_calculator.py b/lstchain/calib/camera/calibration_calculator.py index c801252063..37e6a8fc38 100644 --- a/lstchain/calib/camera/calibration_calculator.py +++ b/lstchain/calib/camera/calibration_calculator.py @@ -61,8 +61,9 @@ class CalibrationCalculator(Component): def __init__( self, + parent=None, + config=None, **kwargs - ): """ @@ -80,19 +81,17 @@ def __init__( The pedestal to use. If None, then PedestalCalculator will be used by default. - kwargs - """ - super().__init__(**kwargs) + super().__init__(parent=parent, config=config,**kwargs) self.flatfield = FlatFieldCalculator.from_name( self.flatfield_product, - **kwargs + parent=self ) self.pedestal = PedestalCalculator.from_name( self.pedestal_product, - **kwargs + parent=self ) msg = "tel_id not the same for all calibration components" @@ -133,23 +132,6 @@ class LSTCalibrationCalculator(CalibrationCalculator): help='Temporary cut on LG std against Lidar events till the calibox TIB do not work (default for filter 5.2) ' ).tag(config=True) - def __init__(self, **kwargs): - """ - Calibration calculator for LST camera - Fills the MonitoringCameraContainer on the base of calibration events - - Parameters: - ---------- - minimum_hg_charge_median : - Temporary cut on HG charge till the calibox TIB do not work - (default for filter 5.2) - - maximum_lg_charge_std - Temporary cut on LG std against Lidar events till the calibox TIB do not work - (default for filter 5.2) - - """ - super().__init__(**kwargs) def calculate_calibration_coefficients(self, event): """ @@ -198,9 +180,7 @@ def calculate_calibration_coefficients(self, event): camera_time_median = np.median(ff_data.time_median, axis=1) calib_data.time_correction = -ff_data.relative_time_median - camera_time_median[:, np.newaxis] - ped_extractor_name = self.config.get("PedestalCalculator").get("charge_product") - ped_width = self.config.get(ped_extractor_name).get("window_width") - calib_data.pedestal_per_sample = ped_data.charge_median / ped_width + calib_data.pedestal_per_sample = ped_data.charge_median / self.pedestal.extractor.window_width # put to zero unusable pixels calib_data.dc_to_pe[calib_data.unusable_pixels] = 0 @@ -209,8 +189,6 @@ def calculate_calibration_coefficients(self, event): # eliminate inf values id any (still necessary?) calib_data.dc_to_pe[np.isinf(calib_data.dc_to_pe)] = 0 - - def process_interleaved(self, event): """ Process interleaved calibration events (pedestals and FF) diff --git a/lstchain/data/lstchain_standard_config.json b/lstchain/data/lstchain_standard_config.json index 9c2881e6a2..bfe8c1ff79 100644 --- a/lstchain/data/lstchain_standard_config.json +++ b/lstchain/data/lstchain_standard_config.json @@ -130,7 +130,7 @@ "maximum_lg_charge_std": 300, "flatfield_product": "FlasherFlatFieldCalculator", "pedestal_product": "PedestalIntegrator", - "PedestalCalculator":{ + "PedestalIntegrator":{ "sample_size": 10000, "sample_duration":100000, "tel_id":1, @@ -138,7 +138,7 @@ "charge_std_cut_outliers": [-10,10], "charge_product":"FixedWindowSum" }, - "FlatFieldCalculator":{ + "FlasherFlatFieldCalculator":{ "sample_size": 10000, "sample_duration":100000, "tel_id":1, diff --git a/lstchain/data/onsite_camera_calibration_param.json b/lstchain/data/onsite_camera_calibration_param.json index b71f884c20..ac47d58c18 100644 --- a/lstchain/data/onsite_camera_calibration_param.json +++ b/lstchain/data/onsite_camera_calibration_param.json @@ -10,14 +10,14 @@ "r0calibrator_product": "LSTR0Corrections", "log_level":"DEBUG" }, - "PedestalCalculator":{ + "PedestalIntegrator":{ "sample_duration":100000, "tel_id":1, "charge_median_cut_outliers": [-10,10], "charge_std_cut_outliers": [-10,10], "charge_product":"FixedWindowSum" }, - "FlatFieldCalculator":{ + "FlasherFlatFieldCalculator":{ "sample_duration":100000, "tel_id":1, "charge_product":"LocalPeakWindowSum", diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py index 7bab4caa99..c46b3541e0 100644 --- a/lstchain/reco/r0_to_dl1.py +++ b/lstchain/reco/r0_to_dl1.py @@ -260,9 +260,11 @@ def r0_to_dl1( # Component to process interleaved pedestal and flat-fields - calib_config = Config(config[config['calibration_product']]) - # time calibration path to flatfield trailet - calib_config.FlatFieldCalculator.time_calibration_path = time_calibration_path + calib_config = Config(config) + + # set time calibration path for flatfield trailet () + calib_config.FlasherFlatFieldCalculator.time_calibration_path = time_calibration_path + calibration_calculator = CalibrationCalculator.from_name( config['calibration_product'], config=calib_config diff --git a/lstchain/scripts/onsite/onsite_create_calibration_file.py b/lstchain/scripts/onsite/onsite_create_calibration_file.py index 807394ad11..769519e3b5 100755 --- a/lstchain/scripts/onsite/onsite_create_calibration_file.py +++ b/lstchain/scripts/onsite/onsite_create_calibration_file.py @@ -136,8 +136,8 @@ def main(): # run lstchain script cmd = f"lstchain_create_calibration_file " \ f"--input_file={input_file} --output_file={output_file} --pedestal_file={pedestal_file} " \ - f"--FlatFieldCalculator.time_calibration_path={time_file} --FlatFieldCalculator.sample_size={stat_events} "\ - f"--PedestalCalculator.sample_size={stat_events} " \ + f"--FlasherFlatFieldCalculator.time_calibration_path={time_file} --FlasherFlatFieldCalculator.sample_size={stat_events} "\ + f"--PedestalIntegrator.sample_size={stat_events} " \ f"--EventSource.max_events={max_events} --config={config_file} > {log_file} 2>&1" print("\n--> RUNNING...") diff --git a/lstchain/tools/lstchain_create_calibration_file.py b/lstchain/tools/lstchain_create_calibration_file.py index b6d9b00134..abf6554715 100644 --- a/lstchain/tools/lstchain_create_calibration_file.py +++ b/lstchain/tools/lstchain_create_calibration_file.py @@ -75,7 +75,6 @@ class CalibrationHDF5Writer(Tool): ] + traits.classes_with_traits(CameraR0Calibrator) + traits.classes_with_traits(CalibrationCalculator) - ) def __init__(self, **kwargs): @@ -97,9 +96,9 @@ def __init__(self, **kwargs): self.simulation = False def setup(self): - kwargs = dict(parent=self) + self.log.debug(f"Open file") - self.eventsource = EventSource.from_config(**kwargs) + self.eventsource = EventSource.from_config(parent=self) # if data remember how many event in the files if "LSTEventSource" in str(type(self.eventsource)): @@ -112,13 +111,13 @@ def setup(self): self.processor = CalibrationCalculator.from_name( self.calibration_product, - **kwargs + parent=self ) if self.r0calibrator_product: self.r0calibrator = CameraR0Calibrator.from_name( self.r0calibrator_product, - **kwargs + parent=self ) From 4de59e4a57ae62261ff45bcde67f116481ada4ba Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 22 Apr 2020 18:57:04 +0200 Subject: [PATCH 34/45] Use the R1CameraContainer.selected_gain_chanel field --- lstchain/calib/camera/calibrator.py | 17 ++--------------- lstchain/io/lstcontainers.py | 2 +- lstchain/reco/r0_to_dl1.py | 2 +- 3 files changed, 4 insertions(+), 17 deletions(-) diff --git a/lstchain/calib/camera/calibrator.py b/lstchain/calib/camera/calibrator.py index a8e3de892b..654b7d5c0c 100644 --- a/lstchain/calib/camera/calibrator.py +++ b/lstchain/calib/camera/calibrator.py @@ -189,21 +189,8 @@ def _calibrate_dl1(self, event, telid): event.dl1.tel[telid].image = charge[gain_mask, np.arange(charge.shape[1])] event.dl1.tel[telid].pulse_time = pulse_time_ff_corrected[gain_mask, np.arange(pulse_time_ff_corrected.shape[1])] - # remember the mask in the lst pixel_status array (this info is missing for the moment in the - # r1 container). I follow the prescription given in the document - # "R1 & DL0 Telescope Event Interfaces and Prototype Evaluation" of K. Kosack - - # bit 2 = LG - gain_mask *= 4 - - # bit 3 = HG - gain_mask[np.where(gain_mask == 0)] = 8 - - # bit 1 = pixel broken pixel (coming from the EvB) - gain_mask += event.lst.tel[telid].evt.pixel_status >> 1 & 1 - - # update pixel status - event.lst.tel[telid].evt.pixel_status = gain_mask + # remember which channel has been selected + event.r1.tel[telid].selected_gain_channel = gain_mask # if threshold == None else: diff --git a/lstchain/io/lstcontainers.py b/lstchain/io/lstcontainers.py index 1c2e33cefd..fffbd2a24a 100644 --- a/lstchain/io/lstcontainers.py +++ b/lstchain/io/lstcontainers.py @@ -223,7 +223,7 @@ class ExtraMCInfo(Container): class ExtraImageInfo(Container): """ attach the tel_id """ tel_id = Field(0, "Telescope ID") - low_gain_mask = Field(None, "pixels with charge from low gain channel") + selected_gain_channel = Field(None, "Selected gain channel") class ThrownEventsHistogram(Container): diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py index c46b3541e0..f154f5fb57 100644 --- a/lstchain/reco/r0_to_dl1.py +++ b/lstchain/reco/r0_to_dl1.py @@ -513,7 +513,7 @@ def r0_to_dl1( # extra info for the image table extra_im.tel_id = telescope_id - extra_im.low_gain_mask = event.lst.tel[telescope_id].evt.pixel_status >> 2 & 1 + extra_im.selected_gain_channel = event.r1.tel[telescope_id].selected_gain_channel for container in [extra_im, dl1_container, event.r0, tel]: add_global_metadata(container, metadata) From 412b3e299b28c72367aefea934d72c8318fc6f9e Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Fri, 24 Apr 2020 08:50:12 +0200 Subject: [PATCH 35/45] Correct program flow in case of MC data --- lstchain/reco/r0_to_dl1.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py index f154f5fb57..9dbbd99076 100644 --- a/lstchain/reco/r0_to_dl1.py +++ b/lstchain/reco/r0_to_dl1.py @@ -593,14 +593,14 @@ def r0_to_dl1( ) - - # at the end of event loop ask calculation of remaining interleaved statistics - calibration_calculator.output_interleaved_results(event) - # write monitoring events - write_calibration_data(writer, - calibration_index, - event.mon.tel[tel_id], - new_ped=True, new_ff=True) + if not is_simu: + # at the end of event loop ask calculation of remaining interleaved statistics + calibration_calculator.output_interleaved_results(event) + # write monitoring events + write_calibration_data(writer, + calibration_index, + event.mon.tel[tel_id], + new_ped=True, new_ff=True) if is_simu: From 1a05e04d1b8fac3bb11362f05c4ecb4344366a2d Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Fri, 24 Apr 2020 13:34:43 +0200 Subject: [PATCH 36/45] Updated calibration json parameter files --- .../data/onsite_camera_calibration_param.json | 59 +++---- lstchain/tools/camera_calibration_param.json | 19 ++- .../lstchain_data_create_pedestal_file.py | 146 ------------------ lstchain/tools/pedestal_param.json | 32 ---- 4 files changed, 43 insertions(+), 213 deletions(-) delete mode 100644 lstchain/tools/lstchain_data_create_pedestal_file.py delete mode 100644 lstchain/tools/pedestal_param.json diff --git a/lstchain/data/onsite_camera_calibration_param.json b/lstchain/data/onsite_camera_calibration_param.json index ac47d58c18..b76b489db5 100644 --- a/lstchain/data/onsite_camera_calibration_param.json +++ b/lstchain/data/onsite_camera_calibration_param.json @@ -2,44 +2,49 @@ "version": 1, "CalibrationHDF5Writer": { - "minimum_hg_charge_median": 5000, - "maximum_lg_charge_std": 300, "one_event": true, - "flatfield_product": "FlasherFlatFieldCalculator", - "pedestal_product": "PedestalIntegrator", + "calibration_product": "LSTCalibrationCalculator", "r0calibrator_product": "LSTR0Corrections", "log_level":"DEBUG" - }, + }, + "LSTCalibrationCalculator": { + "minimum_hg_charge_median": 5000, + "maximum_lg_charge_std": 300, + "flatfield_product": "FlasherFlatFieldCalculator", + "pedestal_product": "PedestalIntegrator" + }, "PedestalIntegrator":{ - "sample_duration":100000, - "tel_id":1, - "charge_median_cut_outliers": [-10,10], - "charge_std_cut_outliers": [-10,10], - "charge_product":"FixedWindowSum" + "sample_size": 10000, + "sample_duration":100000, + "tel_id":1, + "charge_median_cut_outliers": [-10,10], + "charge_std_cut_outliers": [-10,10], + "charge_product":"FixedWindowSum" }, "FlasherFlatFieldCalculator":{ - "sample_duration":100000, - "tel_id":1, - "charge_product":"LocalPeakWindowSum", - "charge_median_cut_outliers": [-0.5,0.5], - "charge_std_cut_outliers": [-10,10], - "time_cut_outliers": [2,38] + "sample_size": 10000, + "sample_duration":100000, + "tel_id":1, + "charge_product":"LocalPeakWindowSum", + "charge_median_cut_outliers": [-0.5,0.5], + "charge_std_cut_outliers": [-10,10], + "time_cut_outliers": [2,38] }, - "LSTR0Corrections": { + "LSTR0Corrections": { "r1_sample_start": 2, - "tel_id":1, + "tel_id":1, "r1_sample_end": 38 }, - "TimeCorrectionCalculate": { + "TimeCorrectionCalculate": { "tel_id":1, "charge_product":"LocalPeakWindowSum" }, - "LocalPeakWindowSum":{ - "window_shift": 5, - "window_width":12 - }, - "FixedWindowSum":{ - "window_start": 12, - "window_width":12 - } + "LocalPeakWindowSum":{ + "window_shift": 5, + "window_width":12 + }, + "FixedWindowSum":{ + "window_start": 12, + "window_width":12 + } } diff --git a/lstchain/tools/camera_calibration_param.json b/lstchain/tools/camera_calibration_param.json index ed9e4d020b..ae99da7a62 100644 --- a/lstchain/tools/camera_calibration_param.json +++ b/lstchain/tools/camera_calibration_param.json @@ -2,14 +2,17 @@ "version": 1, "CalibrationHDF5Writer": { - "minimum_charge": 800, - "flatfield_product": "FlasherFlatFieldCalculator", - "pedestal_product": "PedestalIntegrator", + "one_event": true, + "calibration_product": "LSTCalibrationCalculator", "r0calibrator_product": "LSTR0Corrections", - "output_file":"calibration.hdf5", - "log_file": "log.txt", "log_level":"DEBUG" - }, + }, + "LSTCalibrationCalculator": { + "minimum_hg_charge_median": 5000, + "maximum_lg_charge_std": 300, + "flatfield_product": "FlasherFlatFieldCalculator", + "pedestal_product": "PedestalIntegrator" + }, "EventSource": { "max_events": 1000, "input_url": "/ctadata/franca/LST/LST-1.1.Run00472.0000.fits.fz" @@ -17,7 +20,7 @@ "PedestalCalculator":{ "sample_size": 100, "sample_duration":1000, - "tel_id":0, + "tel_id":1, "charge_median_cut_outliers": [-5,5], "charge_std_cut_outliers": [-5,5], "charge_product":"FixedWindowSum" @@ -33,7 +36,7 @@ }, "LSTR0Corrections": { "pedestal_path": "/ctadata/franca/LST/pedestal_file_run446_0000.fits", - "tel_id": 0, + "tel_id": 1, "r1_sample_start": 2, "r1_sample_end": 38 }, diff --git a/lstchain/tools/lstchain_data_create_pedestal_file.py b/lstchain/tools/lstchain_data_create_pedestal_file.py deleted file mode 100644 index 7c986dc856..0000000000 --- a/lstchain/tools/lstchain_data_create_pedestal_file.py +++ /dev/null @@ -1,146 +0,0 @@ -""" -Extract pedestals from pedestal file -""" -from traitlets import Dict, List, Unicode - -from ctapipe.core import Provenance -from ctapipe.io import HDF5TableWriter -from ctapipe.core import Tool, traits -from ctapipe.io import EventSource - - -from ctapipe.calib.camera.pedestals import PedestalCalculator -from ctapipe.io.containers import PedestalContainer -from lstchain.calib.camera.r0 import CameraR0Calibrator - -__all__ = ['PedestalCalculator', - 'PedestalHDF5Writer', - 'PedestalContainer', - ] - - -class PedestalHDF5Writer(Tool): - ''' - Example of tool that extract the pedestal value per pixel and write the pedestal - container to disk in a hdf5 file - ''' - - name = "PedestalHDF5Writer" - description = "Generate a HDF5 file with pedestal values" - - output_file = Unicode( - 'pedestal.hdf5', - help='Name of the output file' - ).tag(config=True) - - calculator_product = traits.enum_trait( - PedestalCalculator, - default='PedestalIntegrator' - ) - r0calibrator_product = traits.enum_trait( - CameraR0Calibrator, - default='NullR0Calibrator' - ) - - aliases = Dict(dict( - input_file='EventSource.input_url', - output_file='CalibrationHDF5Writer.output_file', - max_events='EventSource.max_events', - tel_id='PedestalCalculator.tel_id', - sample_duration='PedestalCalculator.sample_duration', - sample_size='PedestalCalculator.sample_size', - n_channels='PedestalCalculator.n_channels', - charge_product = 'PedestalCalculator.charge_product' - )) - - classes = List([EventSource, - PedestalCalculator, - PedestalContainer, - CameraR0Calibrator, - HDF5TableWriter - ] + traits.classes_with_traits(PedestalCalculator) - + traits.classes_with_traits(CameraR0Calibrator)) - - def __init__(self, **kwargs): - ''' - Example of tool that extract the pedestal value per pixel and write the pedestal - container to disk - ''' - - super().__init__(**kwargs) - self.eventsource = None - self.pedestal = None - self.container = None - self.writer = None - self.group_name = None - self.r0calibrator = None - - def setup(self): - - kwargs = dict(parent=self) - self.eventsource = EventSource.from_config(**kwargs) - self.pedestal = PedestalCalculator.from_name( - self.calculator_product, - **kwargs - ) - self.r0calibrator = CameraR0Calibrator.from_name( - self.r0calibrator_product, - **kwargs - ) - self.group_name = 'tel_' + str(self.pedestal.tel_id) - - self.writer = HDF5TableWriter( - filename=self.output_file, group_name=self.group_name, overwrite=True - ) - - def start(self, tel_id=1): - ''' - Example of tool that extract the pedestal value per pixel and write the pedestal - container to disk - ''' - - - write_config = True - - # loop on events - for count, event in enumerate(self.eventsource): - # select only pedestal events - - if event.r0.tel[self.pedestal.tel_id].trigger_type != 32: - continue - - # perform R0->R1 - self.r0calibrator.calibrate(event) - - # fill pedestal monitoring container - if self.pedestal.calculate_pedestals(event): - - ped_data = event.mon.tel[self.pedestal.tel_id].pedestal - - self.log.debug(f" r0 {event.r0.tel[tel_id].waveform.shape}") - self.log.debug(f" r1 {event.r1.tel[tel_id].waveform.shape}") - - if write_config: - ped_data.meta['config']= self.config - write_config = False - - self.log.debug(f"write event in table: /{self.group_name}/pedestal") - - # write data to file - self.writer.write('pedestal', ped_data) - - def finish(self): - Provenance().add_output_file( - self.output_file, - role='mon.tel.pedestal' - ) - self.writer.close() - - -def main(): - exe = PedestalHDF5Writer() - exe.run() - - -if __name__ == '__main__': - main() diff --git a/lstchain/tools/pedestal_param.json b/lstchain/tools/pedestal_param.json deleted file mode 100644 index 2a6ca6020c..0000000000 --- a/lstchain/tools/pedestal_param.json +++ /dev/null @@ -1,32 +0,0 @@ -{ - "version": 1, - - "PedestalHDF5Writer": { - "calculator_product": "PedestalIntegrator", - "r0calibrator_product": "LSTR0Corrections", - "output_file":"/astro/users/cassol/soft/python/lstchain-test/pedestal.hdf5", - "log_level":"DEBUG" - }, - "EventSource": { - "input_url": "/ctadata/franca/fefs/aswg/real/R0/20191107/LST-1.1.Run01579.0000.fits.fz", - "max_events": 1000 - }, - "LSTR0Corrections": { - "pedestal_path": "/ctadata/franca/fefs/aswg/real/calibration/20191107/v00/drs4_pedestal.Run1579.0000.fits", - "r1_sample_start": 2, - "r1_sample_end": 38, - "tel_id": 1 - }, - "PedestalCalculator":{ - "sample_size": 100, - "sample_duration":1000, - "tel_id":1, - "charge_median_cut_outliers": [-3,3], - "charge_std_cut_outliers": [-3,3], - "charge_product":"FixedWindowSum" - }, - "FixedWindowSum":{ - "window_start": 12, - "window_width":12 - } -} From 287f9a4b1f97a6599d54a622bba661d4e5b6ee26 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Fri, 24 Apr 2020 14:26:03 +0200 Subject: [PATCH 37/45] Allow time_calibration_path == "None" for simulation events --- lstchain/calib/camera/flatfield.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/lstchain/calib/camera/flatfield.py b/lstchain/calib/camera/flatfield.py index 9dd8ac4d62..8737123442 100644 --- a/lstchain/calib/camera/flatfield.py +++ b/lstchain/calib/camera/flatfield.py @@ -78,12 +78,19 @@ def __init__(self, **kwargs): self.arrival_times = None # arrival time per event in sample self.sample_masked_pixels = None # masked pixels per event in sample - # declare time calibrator if correction file exist - if os.path.exists(self.time_calibration_path): - self.time_corrector = PulseTimeCorrection( - calib_file_path = self.time_calibration_path) + # declare time calibrator None of MC events + if self.time_calibration_path == "None": + self.time_corrector = None + else: - raise IOError(f"Time calibration file {self.time_calibration_path} not found!") + # look for calibration path otherwise + if os.path.exists(self.time_calibration_path): + self.time_corrector = PulseTimeCorrection( + calib_file_path = self.time_calibration_path) + else: + msg=f"Time calibration file {self.time_calibration_path} not found!" + raise IOError(msg) + def _extract_charge(self, event): From ffc249c8a7723295835e9a2f10b9014f63c9a676 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Fri, 24 Apr 2020 14:46:04 +0200 Subject: [PATCH 38/45] Use null in jason file instead of "None" --- lstchain/calib/camera/flatfield.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lstchain/calib/camera/flatfield.py b/lstchain/calib/camera/flatfield.py index 8737123442..9f238b7205 100644 --- a/lstchain/calib/camera/flatfield.py +++ b/lstchain/calib/camera/flatfield.py @@ -45,6 +45,7 @@ class FlasherFlatFieldCalculator(FlatFieldCalculator): ).tag(config=True) time_calibration_path = Unicode( '', + allow_none = True, help = 'Path to drs4 time calibration file' ).tag(config = True) @@ -78,10 +79,9 @@ def __init__(self, **kwargs): self.arrival_times = None # arrival time per event in sample self.sample_masked_pixels = None # masked pixels per event in sample - # declare time calibrator None of MC events - if self.time_calibration_path == "None": + # declare time calibrator None if MC events + if self.time_calibration_path is None: self.time_corrector = None - else: # look for calibration path otherwise if os.path.exists(self.time_calibration_path): From 1e7d05797ae784557dd264fa98ade29b61c50cbb Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Tue, 28 Apr 2020 17:49:16 +0200 Subject: [PATCH 39/45] Correct config definition position --- .../onsite/onsite_create_calibration_file.py | 16 +++++++--------- .../tools/lstchain_create_calibration_file.py | 3 --- 2 files changed, 7 insertions(+), 12 deletions(-) diff --git a/lstchain/scripts/onsite/onsite_create_calibration_file.py b/lstchain/scripts/onsite/onsite_create_calibration_file.py index 769519e3b5..676ee467ec 100755 --- a/lstchain/scripts/onsite/onsite_create_calibration_file.py +++ b/lstchain/scripts/onsite/onsite_create_calibration_file.py @@ -77,6 +77,13 @@ def main(): print(f">>> Error: The pedestal file {pedestal_file} do not exist.\n Exit") exit(0) + # define config file + config_file = os.path.join(os.path.dirname(__file__), "../../data/onsite_camera_calibration_param.json") + if not os.path.exists(config_file): + print(f">>> Config file {config_file} do not exists. \n Exit ") + exit(1) + print(f"\n--> Config file {config_file}") + # # produce drs4 time calibration file @@ -120,15 +127,6 @@ def main(): print(f"\n--> Log file {log_file}") - # define config file - config_file = os.path.join(os.path.dirname(__file__), "../../data/onsite_camera_calibration_param.json") - if not os.path.exists(config_file): - print(f">>> Config file {config_file} do not exists. \n Exit ") - exit(1) - print(f"\n--> Config file {config_file}") - - - # # produce ff calibration file # diff --git a/lstchain/tools/lstchain_create_calibration_file.py b/lstchain/tools/lstchain_create_calibration_file.py index abf6554715..533f4a19b7 100644 --- a/lstchain/tools/lstchain_create_calibration_file.py +++ b/lstchain/tools/lstchain_create_calibration_file.py @@ -102,7 +102,6 @@ def setup(self): # if data remember how many event in the files if "LSTEventSource" in str(type(self.eventsource)): - self.tot_events = len(self.eventsource.multi_file) self.log.debug(f"Input file has file {self.tot_events} events") else: @@ -180,7 +179,6 @@ def start(self): self.simulation and np.median(np.sum(event.r1.tel[tel_id].waveform[0], axis=1)) < self.processor.minimum_hg_charge_median): - new_ped = self.processor.pedestal.calculate_pedestals(event) @@ -193,7 +191,6 @@ def start(self): < self.processor.maximum_lg_charge_std): new_ff = self.processor.flatfield.calculate_relative_gain(event) - # write pedestal results when enough statistics or end of file if new_ped or end_of_file: From 87bcec091b0c8ddb04c4d108d0b160b37ea06e07 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Tue, 28 Apr 2020 17:57:59 +0200 Subject: [PATCH 40/45] Delete unused script --- .../lstchain_data_create_pedestal_file.py | 147 ------------------ 1 file changed, 147 deletions(-) delete mode 100644 lstchain/tools/lstchain_data_create_pedestal_file.py diff --git a/lstchain/tools/lstchain_data_create_pedestal_file.py b/lstchain/tools/lstchain_data_create_pedestal_file.py deleted file mode 100644 index d411fbd784..0000000000 --- a/lstchain/tools/lstchain_data_create_pedestal_file.py +++ /dev/null @@ -1,147 +0,0 @@ -""" -Extract pedestals from pedestal file -""" -from traitlets import Dict, List, Unicode - -from ctapipe.core import Provenance -from ctapipe.io import HDF5TableWriter -from ctapipe.core import Tool, traits -from ctapipe.io import EventSource - - -from ctapipe.calib.camera.pedestals import PedestalCalculator -from ctapipe.io.containers import PedestalContainer -from lstchain.calib.camera.r0 import CameraR0Calibrator - -__all__ = [ - 'PedestalCalculator', - 'PedestalHDF5Writer', - 'PedestalContainer', - ] - - -class PedestalHDF5Writer(Tool): - ''' - Example of tool that extract the pedestal value per pixel and write the pedestal - container to disk in a hdf5 file - ''' - - name = "PedestalHDF5Writer" - description = "Generate a HDF5 file with pedestal values" - - output_file = Unicode( - 'pedestal.hdf5', - help='Name of the output file' - ).tag(config=True) - - calculator_product = traits.enum_trait( - PedestalCalculator, - default='PedestalIntegrator' - ) - r0calibrator_product = traits.enum_trait( - CameraR0Calibrator, - default='NullR0Calibrator' - ) - - aliases = Dict(dict( - input_file='EventSource.input_url', - output_file='CalibrationHDF5Writer.output_file', - max_events='EventSource.max_events', - tel_id='PedestalCalculator.tel_id', - sample_duration='PedestalCalculator.sample_duration', - sample_size='PedestalCalculator.sample_size', - n_channels='PedestalCalculator.n_channels', - charge_product = 'PedestalCalculator.charge_product' - )) - - classes = List([EventSource, - PedestalCalculator, - PedestalContainer, - CameraR0Calibrator, - HDF5TableWriter - ] + traits.classes_with_traits(PedestalCalculator) - + traits.classes_with_traits(CameraR0Calibrator)) - - def __init__(self, **kwargs): - ''' - Example of tool that extract the pedestal value per pixel and write the pedestal - container to disk - ''' - - super().__init__(**kwargs) - self.eventsource = None - self.pedestal = None - self.container = None - self.writer = None - self.group_name = None - self.r0calibrator = None - - def setup(self): - - kwargs = dict(parent=self) - self.eventsource = EventSource.from_config(**kwargs) - self.pedestal = PedestalCalculator.from_name( - self.calculator_product, - **kwargs - ) - self.r0calibrator = CameraR0Calibrator.from_name( - self.r0calibrator_product, - **kwargs - ) - self.group_name = 'tel_' + str(self.pedestal.tel_id) - - self.writer = HDF5TableWriter( - filename=self.output_file, group_name=self.group_name, overwrite=True - ) - - def start(self, tel_id=1): - ''' - Example of tool that extract the pedestal value per pixel and write the pedestal - container to disk - ''' - - - write_config = True - - # loop on events - for count, event in enumerate(self.eventsource): - # select only pedestal events - - if event.r0.tel[self.pedestal.tel_id].trigger_type != 32: - continue - - # perform R0->R1 - self.r0calibrator.calibrate(event) - - # fill pedestal monitoring container - if self.pedestal.calculate_pedestals(event): - - ped_data = event.mon.tel[self.pedestal.tel_id].pedestal - - self.log.debug(f" r0 {event.r0.tel[tel_id].waveform.shape}") - self.log.debug(f" r1 {event.r1.tel[tel_id].waveform.shape}") - - if write_config: - ped_data.meta['config']= self.config - write_config = False - - self.log.debug(f"write event in table: /{self.group_name}/pedestal") - - # write data to file - self.writer.write('pedestal', ped_data) - - def finish(self): - Provenance().add_output_file( - self.output_file, - role='mon.tel.pedestal' - ) - self.writer.close() - - -def main(): - exe = PedestalHDF5Writer() - exe.run() - - -if __name__ == '__main__': - main() From 018163e52e9be83848c19ccfd8d617483f788d73 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Tue, 28 Apr 2020 18:16:34 +0200 Subject: [PATCH 41/45] Correct import --- lstchain/io/__init__.py | 2 ++ lstchain/io/io.py | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/lstchain/io/__init__.py b/lstchain/io/__init__.py index 06a7f8cd57..ee106cc75a 100644 --- a/lstchain/io/__init__.py +++ b/lstchain/io/__init__.py @@ -16,6 +16,7 @@ write_metadata, write_subarray_tables, read_simu_info_merged_hdf5, + write_calibration_data ) standard_config = get_standard_config() @@ -38,4 +39,5 @@ 'write_metadata', 'write_subarray_tables', 'read_simu_info_merged_hdf5', + 'write_calibration_data' ] diff --git a/lstchain/io/io.py b/lstchain/io/io.py index 219348efe5..c3624b2ca1 100644 --- a/lstchain/io/io.py +++ b/lstchain/io/io.py @@ -671,7 +671,7 @@ def recursive_copy_node(src_file, dir_file, path): recursive_path = os.path.join(recursive_path, p) -def write_calibration_data(writer,mon_index, mon_event, new_ped=False, new_ff=False): +def write_calibration_data(writer, mon_index, mon_event, new_ped=False, new_ff=False): mon_event.pedestal.prefix = '' mon_event.flatfield.prefix = '' mon_event.calibration.prefix = '' From 9887df34ef760e1666e0d6eeb978f99855350eaf Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 29 Apr 2020 09:43:43 +0200 Subject: [PATCH 42/45] Monir corrections --- lstchain/calib/camera/calibration_calculator.py | 4 +++- lstchain/calib/camera/flatfield.py | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/lstchain/calib/camera/calibration_calculator.py b/lstchain/calib/camera/calibration_calculator.py index 37e6a8fc38..9560d3a425 100644 --- a/lstchain/calib/camera/calibration_calculator.py +++ b/lstchain/calib/camera/calibration_calculator.py @@ -95,7 +95,9 @@ def __init__( ) msg = "tel_id not the same for all calibration components" - assert self.pedestal.tel_id == self.flatfield.tel_id, msg + if self.pedestal.tel_id != self.flatfield.tel_id: + raise ValueError(msg) + self.tel_id = self.flatfield.tel_id diff --git a/lstchain/calib/camera/flatfield.py b/lstchain/calib/camera/flatfield.py index 9f238b7205..f5b58eefa9 100644 --- a/lstchain/calib/camera/flatfield.py +++ b/lstchain/calib/camera/flatfield.py @@ -6,7 +6,7 @@ import numpy as np from astropy import units as u from ctapipe.calib.camera.flatfield import FlatFieldCalculator -from ctapipe.core.traits import List, Unicode, Path +from ctapipe.core.traits import List, Unicode from lstchain.calib.camera.pulse_time_correction import PulseTimeCorrection __all__ = [ From f21f6c6a5773f52376ebd356c9e90cdf1a756633 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Wed, 29 Apr 2020 10:21:48 +0200 Subject: [PATCH 43/45] Update imports --- lstchain/scripts/onsite/onsite_create_calibration_file.py | 1 + lstchain/tools/lstchain_create_calibration_file.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/lstchain/scripts/onsite/onsite_create_calibration_file.py b/lstchain/scripts/onsite/onsite_create_calibration_file.py index 676ee467ec..5a27710184 100755 --- a/lstchain/scripts/onsite/onsite_create_calibration_file.py +++ b/lstchain/scripts/onsite/onsite_create_calibration_file.py @@ -9,6 +9,7 @@ """ import argparse +import os from pathlib import Path from lstchain.io.data_management import * import lstchain.visualization.plot_calib as calib diff --git a/lstchain/tools/lstchain_create_calibration_file.py b/lstchain/tools/lstchain_create_calibration_file.py index 533f4a19b7..b20df0f416 100644 --- a/lstchain/tools/lstchain_create_calibration_file.py +++ b/lstchain/tools/lstchain_create_calibration_file.py @@ -11,7 +11,7 @@ from ctapipe.io import EventSource from ctapipe.io.containers import PixelStatusContainer from lstchain.calib.camera.calibration_calculator import CalibrationCalculator -from lstchain.calib.camera import CameraR0Calibrator +from lstchain.calib.camera.r0 import CameraR0Calibrator from lstchain.io.lstcontainers import LSTEventType __all__ = [ From 8d016d0776be12a1173a7b624d9926e4573a1810 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Thu, 30 Apr 2020 12:42:20 +0200 Subject: [PATCH 44/45] Skip the first 1000 events that badly drs4 corrected --- lstchain/tools/lstchain_create_calibration_file.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/lstchain/tools/lstchain_create_calibration_file.py b/lstchain/tools/lstchain_create_calibration_file.py index b20df0f416..2719bd3766 100644 --- a/lstchain/tools/lstchain_create_calibration_file.py +++ b/lstchain/tools/lstchain_create_calibration_file.py @@ -136,6 +136,9 @@ def start(self): new_ff = False end_of_file = False + # skip the first events which are badly drs4 corrected + events_to_skip = 1000 + try: self.log.debug(f"Start loop") for count, event in enumerate(self.eventsource): @@ -169,12 +172,15 @@ def start(self): # correct for low level calibration self.r0calibrator.calibrate(event) + # skip first events which are badly drs4 corrected + if not self.simulation and count < events_to_skip: + continue + # reject event without trigger type if LSTEventType.is_unknown(event.r1.tel[tel_id].trigger_type): continue # if pedestal event - if LSTEventType.is_pedestal(event.r1.tel[tel_id].trigger_type) or ( self.simulation and np.median(np.sum(event.r1.tel[tel_id].waveform[0], axis=1)) From 491fdb76a16bf93732722e825c1cba71db6e4630 Mon Sep 17 00:00:00 2001 From: Franca Cassol Date: Sat, 2 May 2020 09:38:29 +0200 Subject: [PATCH 45/45] Set None as default for time calibration file --- lstchain/calib/camera/flatfield.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lstchain/calib/camera/flatfield.py b/lstchain/calib/camera/flatfield.py index f5b58eefa9..dfa4cbe6bd 100644 --- a/lstchain/calib/camera/flatfield.py +++ b/lstchain/calib/camera/flatfield.py @@ -44,7 +44,7 @@ class FlasherFlatFieldCalculator(FlatFieldCalculator): help='Interval (number of std) of accepted charge standard deviation around camera median value' ).tag(config=True) time_calibration_path = Unicode( - '', + None, allow_none = True, help = 'Path to drs4 time calibration file' ).tag(config = True) @@ -79,7 +79,7 @@ def __init__(self, **kwargs): self.arrival_times = None # arrival time per event in sample self.sample_masked_pixels = None # masked pixels per event in sample - # declare time calibrator None if MC events + if self.time_calibration_path is None: self.time_corrector = None else: