diff --git a/ctapipe/calib/camera/__init__.py b/ctapipe/calib/camera/__init__.py index 385056a67a3..9d8467b184e 100644 --- a/ctapipe/calib/camera/__init__.py +++ b/ctapipe/calib/camera/__init__.py @@ -7,3 +7,5 @@ from .dl0 import * from .dl1 import * from .calibrator import * + + diff --git a/ctapipe/calib/camera/flatfield.py b/ctapipe/calib/camera/flatfield.py new file mode 100644 index 00000000000..f4de088bca2 --- /dev/null +++ b/ctapipe/calib/camera/flatfield.py @@ -0,0 +1,315 @@ +""" +Factory for the estimation of the flat field coefficients +""" + +from abc import abstractmethod +import numpy as np +from astropy import units as u +from ctapipe.core import Component, Factory + +from ctapipe.image import ChargeExtractorFactory, WaveformCleanerFactory +from ctapipe.core.traits import Int +from ctapipe.io.containers import FlatFieldCameraContainer + +__all__ = [ + 'FlatFieldCalculator', + 'FlasherFlatFieldCalculator', + 'FlatFieldFactory' +] + + +class FlatFieldCalculator(Component): + """ + Parent class for the flat field calculators. + Fills the MON.flatfield container. + """ + + tel_id = Int( + 0, + help='id of the telescope to calculate the flat-field coefficients' + ).tag(config=True) + sample_duration = Int( + 60, + help='sample duration in seconds' + ).tag(config=True) + sample_size = Int( + 10000, + help='sample size' + ).tag(config=True) + n_channels = Int( + 2, + help='number of channels to be treated' + ).tag(config=True) + + def __init__( + self, + config=None, + tool=None, + extractor_product=None, + cleaner_product=None, + **kwargs + ): + """ + Parent class for the flat field calculators. + Fills the MON.flatfield container. + + Parameters + ---------- + config : traitlets.loader.Config + Configuration specified by config file or cmdline arguments. + Used to set traitlet values. + Set to None if no configuration to pass. + tool : ctapipe.core.Tool + Tool executable that is calling this component. + Passes the correct logger to the component. + Set to None if no Tool to pass. + extractor_product : str + The ChargeExtractor to use. + cleaner_product : str + The WaveformCleaner to use. + kwargs + + """ + super().__init__(config=config, parent=tool, **kwargs) + + # initialize the output + self.container = FlatFieldCameraContainer() + + # load the waveform charge extractor and cleaner + kwargs_ = dict() + if extractor_product: + kwargs_['product'] = extractor_product + self.extractor = ChargeExtractorFactory.produce( + config=config, + tool=tool, + **kwargs_ + ) + self.log.info(f"extractor {self.extractor}") + + kwargs_ = dict() + if cleaner_product: + kwargs_['product'] = cleaner_product + self.cleaner = WaveformCleanerFactory.produce( + config=config, + tool=tool, + **kwargs_ + ) + self.log.info(f"cleaner {self.cleaner}") + + @abstractmethod + def calculate_relative_gain(self, event): + """calculate relative gain from event + Parameters + ---------- + event: DataContainer + + Returns: FlatFieldCameraContainer or None + + None is returned if no new flat field coefficients were calculated + e.g. due to insufficient statistics. + """ + + +class FlasherFlatFieldCalculator(FlatFieldCalculator): + + def __init__(self, config=None, tool=None, **kwargs): + """Calculates flat field coefficients from flasher data + + based on the best algorithm described by S. Fegan in MST-CAM-TN-0060 + + Parameters: see base class FlatFieldCalculator + """ + super().__init__(config=config, tool=tool, **kwargs) + + self.log.info("Used events statistics : %d", self.sample_size) + + # members to keep state in calculate_relative_gain() + self.num_events_seen = 0 + self.time_start = None # trigger time of first event in sample + self.charge_medians = None # med. charge in camera per event in sample + self.charges = None # charge per event in sample + self.arrival_times = None # arrival time per event in sample + self.sample_bad_pixels = None # bad pixels per event in sample + + def _extract_charge(self, event): + """ + Extract the charge and the time from a calibration event + + Parameters + ---------- + event : general event container + + """ + + waveforms = event.r0.tel[self.tel_id].waveform + + # Clean waveforms + if self.cleaner: + cleaned = self.cleaner.apply(waveforms) + # do nothing + else: + cleaned = waveforms + + # Extract charge and time + if self.extractor: + if self.extractor.requires_neighbours(): + g = event.inst.subarray.tel[self.tel_id].camera + self.extractor.neighbours = g.neighbor_matrix_where + + charge, peak_pos, window = self.extractor.extract_charge(cleaned) + + # sum all the samples + else: + charge = cleaned.sum(axis=2) + peak_pos = np.argmax(cleaned, axis=2) + + return charge, peak_pos + + def calculate_relative_gain(self, event): + """ + calculate the relative flat field coefficients + + Parameters + ---------- + event : general event container + + """ + + # initialize the np array at each cycle + waveform = event.r0.tel[self.tel_id].waveform + + + # patches for MC data + if not event.mcheader.simtel_version: + trigger_time = event.r0.tel[self.tel_id].trigger_time + pixel_status = event.r0.tel[self.tel_id].pixel_status + else: + trigger_time = event.trig.gps_time.unix + pixel_status = np.ones(waveform.shape[1]) + + + if self.num_events_seen == 0: + self.time_start = trigger_time + self.setup_sample_buffers(waveform, self.sample_size) + + # extract the charge of the event and + # the peak position (assumed as time for the moment) + charge, arrival_time = self._extract_charge(event) + self.collect_sample(charge, pixel_status, arrival_time) + + sample_age = trigger_time - self.time_start + + + # check if to create a calibration event + if ( + sample_age > self.sample_duration + or self.num_events_seen == self.sample_size + ): + relative_gain_results = calculate_relative_gain_results( + self.charge_medians, + self.charges, + self.sample_bad_pixels, + ) + time_results = calculate_time_results( + self.arrival_times, + self.sample_bad_pixels, + self.time_start, + trigger_time, + ) + + result = { + 'n_events': self.num_events_seen, + **relative_gain_results, + **time_results, + } + for key, value in result.items(): + setattr(self.container, key, value) + + self.num_events_seen = 0 + return self.container + + else: + + return None + + def setup_sample_buffers(self, waveform, sample_size): + n_channels = waveform.shape[0] + n_pix = waveform.shape[1] + shape = (sample_size, n_channels, n_pix) + + self.charge_medians = np.zeros((sample_size, n_channels)) + self.charges = np.zeros(shape) + self.arrival_times = np.zeros(shape) + self.sample_bad_pixels = np.zeros(shape) + + def collect_sample(self, charge, pixel_status, arrival_time): + + # extract the charge of the event and + # the peak position (assumed as time for the moment) + bad_pixels = np.zeros(charge.shape, dtype=np.bool) + bad_pixels[:] = pixel_status == 0 + + good_charge = np.ma.array(charge, mask=bad_pixels) + charge_median = np.ma.median(good_charge, axis=1) + + self.charges[self.num_events_seen] = charge + self.arrival_times[self.num_events_seen] = arrival_time + self.sample_bad_pixels[self.num_events_seen] = bad_pixels + self.charge_medians[self.num_events_seen] = charge_median + self.num_events_seen += 1 + + +def calculate_time_results( + trace_time, + bad_pixels_of_sample, + time_start, + trigger_time, +): + masked_trace_time = np.ma.array( + trace_time, + mask=bad_pixels_of_sample + ) + + # extract the average time over the camera and the events + camera_time_median = np.ma.median(masked_trace_time) + camera_time_mean = np.ma.mean(masked_trace_time) + pixel_time_median = np.ma.median(masked_trace_time, axis=0) + pixel_time_mean = np.ma.mean(masked_trace_time, axis=0) + + return { + 'time_mean': (trigger_time - time_start) / 2 * u.s, + 'time_range': [time_start, trigger_time] * u.s, + 'relative_time_median': np.ma.getdata( + pixel_time_median - camera_time_median), + 'relative_time_mean': np.ma.getdata( + pixel_time_mean - camera_time_mean), + } + + +def calculate_relative_gain_results( + event_median, + trace_integral, + bad_pixels_of_sample, +): + masked_trace_integral = np.ma.array( + trace_integral, + mask=bad_pixels_of_sample + ) + relative_gain_event = np.ma.getdata( + masked_trace_integral / event_median[:, :, np.newaxis] + ) + + return { + 'relative_gain_median': np.median(relative_gain_event, axis=0), + 'relative_gain_mean': np.mean(relative_gain_event, axis=0), + 'relative_gain_rms': np.std(relative_gain_event, axis=0), + } + + +class FlatFieldFactory(Factory): + """ + Factory to obtain flat-field coefficients + """ + base = FlatFieldCalculator + default = 'FlasherFlatFieldCalculator' + custom_product_help = ('Flat-flield method to use') diff --git a/ctapipe/calib/pedestals.py b/ctapipe/calib/camera/pedestals.py similarity index 100% rename from ctapipe/calib/pedestals.py rename to ctapipe/calib/camera/pedestals.py diff --git a/ctapipe/calib/camera/tests/test_flatfield.py b/ctapipe/calib/camera/tests/test_flatfield.py new file mode 100644 index 00000000000..1308471c1de --- /dev/null +++ b/ctapipe/calib/camera/tests/test_flatfield.py @@ -0,0 +1,29 @@ +import pytest +from ctapipe.utils import get_dataset_path +from ctapipe.io.nectarcameventsource import NectarCAMEventSource +from ctapipe.calib.camera.flatfield import FlasherFlatFieldCalculator + +# this test cannot run if protozfits is not installed +# since NectarCAMEventSource needs lsteventsource.MultiFiles which +# in turn needs protozfits, which is an optional dependency, so +# in case this is not installed the tests in this file cannot succeed +pytest.importorskip("protozfits", minversion="0.44.3") + + +def test_flasherflatfieldcalculator(): + + example_file_path = get_dataset_path("NectarCAM.Run0890.10events.fits.fz") + + inputfile_reader = NectarCAMEventSource( + input_url=example_file_path, + max_events=10 + ) + + ff_calculator = FlasherFlatFieldCalculator(sample_size=3, tel_id=0) + + for event in inputfile_reader: + + ff_data = ff_calculator.calculate_relative_gain(event) + + if ff_calculator.num_events_seen == ff_calculator.sample_size: + assert ff_data diff --git a/ctapipe/calib/tests/test_pedestals.py b/ctapipe/calib/camera/tests/test_pedestals.py similarity index 100% rename from ctapipe/calib/tests/test_pedestals.py rename to ctapipe/calib/camera/tests/test_pedestals.py diff --git a/ctapipe/image/waveform_cleaning.py b/ctapipe/image/waveform_cleaning.py index 41491ecf196..5ad6ebf6ba5 100644 --- a/ctapipe/image/waveform_cleaning.py +++ b/ctapipe/image/waveform_cleaning.py @@ -13,7 +13,7 @@ LocalPeakIntegrator) __all__ = ['WaveformCleanerFactory', 'CHECMWaveformCleanerAverage', - 'CHECMWaveformCleanerLocal', + 'CHECMWaveformCleanerLocal', 'BaselineWaveformCleaner', 'NullWaveformCleaner'] @@ -67,6 +67,25 @@ def apply(self, waveforms): return waveforms +class BaselineWaveformCleaner(WaveformCleaner): + """ + Basic waveform cleaner that simply returns the waveform subtracted + from the baseline + """ + baseline_width = Int(10, help='Define then number of samples for estimating the ' + 'baseline').tag(config=True) + + window_shift = Int(0, help='Define the shift of the window on which caluclate ' + 'the baseline.').tag(config=True) + def apply(self, waveforms): + # self.log.debug(f"calculate baseline on first {self.baseline_width} samples") + # Subtract initial baseline + baseline_sub = waveforms - np.mean(waveforms[:,:, self.window_shift:self.baseline_width], + axis=2)[:,:, None] + + return baseline_sub + + class CHECMWaveformCleaner(WaveformCleaner): """ Waveform cleaner used by CHEC-M. diff --git a/ctapipe/io/containers.py b/ctapipe/io/containers.py index 0f04d975e2b..e9bc96eac43 100644 --- a/ctapipe/io/containers.py +++ b/ctapipe/io/containers.py @@ -44,6 +44,11 @@ 'LeakageContainer', 'ConcentrationContainer', 'TimingParametersContainer', + 'FlatFieldCameraContainer', + 'PedestalCameraContainer', + 'FlatFieldContainer', + 'PedestalContainer', + 'MonitorDataContainer' ] @@ -145,7 +150,7 @@ class R0CameraContainer(Container): "(n_channels x n_pixels, n_samples)" )) num_samples = Field(None, "number of time samples for telescope") - + pixel_status = Field(0o0, "status of the pixels") class R0Container(Container): """ @@ -473,8 +478,6 @@ class NectarCAMCameraContainer(Container): "Information") - - class NectarCAMContainer(Container): """ Storage for the NectarCAMCameraContainer for each telescope @@ -487,7 +490,6 @@ class NectarCAMContainer(Container): "map of tel_id to NectarCAMCameraContainer") - class NectarCAMDataContainer(DataContainer): """ Data container including NectarCAM information @@ -495,7 +497,6 @@ class NectarCAMDataContainer(DataContainer): nectarcam = Field(NectarCAMContainer(), "NectarCAM specific Information") - class LSTServiceContainer(Container): """ Container for Fields that are specific to each LST camera configuration @@ -777,3 +778,77 @@ class TimingParametersContainer(Container): """ slope = Field(nan, 'Slope of arrival times along main shower axis') intercept = Field(nan, 'intercept of arrival times along main shower axis') + + +class FlatFieldCameraContainer(Container): + """ + Container for relative camera flat-field coefficients + + """ + + time_mean = Field(0, 'Mean time, seconds since reference', unit=u.s) + time_range = Field( + [], + 'Range of time of the calibration data [t_min, t_max]', + unit=u.s + ) + n_events = Field(0, 'Number of events used for statistics') + relative_gain_mean = Field( + None, + "np array of the relative flat-field coefficient mean (n_chan X N_pix)" + ) + relative_gain_median = Field( + None, + "np array of the relative flat-field coefficient median (n_chan X N_pix)" + ) + relative_gain_rms = Field( + None, + "np array of the relative flat-field coefficient rms (n_chan X N_pix)" + ) + relative_time_mean = Field( + None, + "np array of the relative time mean (n_chan X N_pix)", + unit=u.ns + ) + relative_time_median = Field( + None, + "np array of the relative time median (n_chan X N_pix)", + unit=u.ns) + + +class FlatFieldContainer(Container): + """ + Container for relative flat field coefficients + """ + tels_with_data = Field([], "list of telescopes with data") + tel = Field( + Map(FlatFieldCameraContainer), + "map of tel_id to FlatFiledCameraContainer") + + +class PedestalCameraContainer(Container): + """ + Container for pedestals per camera + """ + # still to be filled + + +class PedestalContainer(Container): + """ + Container for pedestal data + """ + tels_with_data = Field([], "list of telescopes with data") + tel = Field( + Map(PedestalCameraContainer), + "map of tel_id to PedestalCameraContainer") + + +class MonitorDataContainer(Container): + """ + Root container for MON data + """ + flatfield = Field(FlatFieldContainer(), "Relative flat field data") + pedestal = Field(PedestalContainer(), "Pedestal data") + + + diff --git a/ctapipe/io/lsteventsource.py b/ctapipe/io/lsteventsource.py index 0246b5fa921..482c480b97a 100644 --- a/ctapipe/io/lsteventsource.py +++ b/ctapipe/io/lsteventsource.py @@ -4,11 +4,9 @@ Needs protozfits v1.4.2 from github.com/cta-sst-1m/protozfitsreader """ +import glob import numpy as np - from astropy import units as u -import glob -from os import getcwd from ctapipe.core import Provenance from ctapipe.instrument import TelescopeDescription, SubarrayDescription, \ CameraGeometry, OpticsDescription @@ -20,15 +18,11 @@ class LSTEventSource(EventSource): - """ EventSource for LST r0 data. """ - - def __init__(self, config=None, tool=None, **kwargs): - """ Constructor Parameters @@ -47,7 +41,6 @@ def __init__(self, config=None, tool=None, **kwargs): the 'input_url' parameter. """ - # EventSource can not handle file wild cards as input_url # To overcome this we substitute the input_url with first file matching # the specified file mask (copied from MAGICEventSourceROOT). @@ -61,13 +54,17 @@ def __init__(self, config=None, tool=None, **kwargs): super().__init__(config=config, tool=tool, **kwargs) self.file_list = [self.input_url] - self.multi_file = MultiFiles(self.file_list) self.camera_config = self.multi_file.camera_config - self.log.info("Read {} input files".format(self.multi_file.num_inputs())) - + self.log.info( + "Read {} input files".format( + self.multi_file.num_inputs() + ) + ) + def rewind(self): + self.multi_file.rewind() def _generator(self): @@ -76,14 +73,13 @@ def _generator(self): self.data.meta['input_url'] = self.input_url self.data.meta['max_events'] = self.max_events - # fill LST data from the CameraConfig table self.fill_lst_service_container_from_zfile() # Instrument information for tel_id in self.data.lst.tels_with_data: - assert (tel_id == 0) # only LST1 for the moment (id = 0) + assert tel_id == 0 # only LST1 for the moment (id = 0) # optics info from standard optics.fits.gz file optics = OpticsDescription.from_name("LST") @@ -101,7 +97,6 @@ def _generator(self): # LSTs telescope position taken from MC from the moment tel_pos = {tel_id: [50., 50., 16] * u.m} - subarray = SubarrayDescription("LST1 subarray") subarray.tels = tels subarray.positions = tel_pos @@ -120,7 +115,6 @@ def _generator(self): self.fill_r0_container_from_zfile(event) yield self.data - @staticmethod def is_compatible(file_path): from .sst1meventsource import is_fits_in_header @@ -169,6 +163,7 @@ def fill_lst_service_container_from_zfile(self): svc_container.pixel_ids = self.camera_config.expected_pixels_id svc_container.data_model_version = self.camera_config.data_model_version + svc_container.num_modules = self.camera_config.lstcam.num_modules svc_container.module_ids = self.camera_config.lstcam.expected_modules_id svc_container.idaq_version = self.camera_config.lstcam.idaq_version @@ -179,7 +174,6 @@ def fill_lst_service_container_from_zfile(self): def fill_lst_event_container_from_zfile(self, event): - event_container = self.data.lst.tel[self.camera_config.telescope_id].evt event_container.configuration_id = event.configuration_id @@ -216,24 +210,28 @@ def fill_r0_camera_container_from_zfile(self, container, event): "N_chan x N_pix x N_samples= '{}'" .format(event.waveform.shape[0])) + container.pixel_status = np.zeros([self.n_camera_pixels]) + container.pixel_status[self.camera_config.expected_pixels_id] = \ + event.pixel_status reshaped_waveform = np.array( - event.waveform - ).reshape(n_gains, - self.camera_config.num_pixels, - container.num_samples) + event.waveform + ).reshape( + n_gains, + self.camera_config.num_pixels, + container.num_samples + ) # initialize the waveform container to zero container.waveform = np.zeros([n_gains, self.n_camera_pixels, container.num_samples]) - # re-order the waveform following the expected_pixels_id values (rank = pixel id) + # re-order the waveform following the expected_pixels_id values + # (rank = pixel id) container.waveform[:, self.camera_config.expected_pixels_id, :] =\ reshaped_waveform def fill_r0_container_from_zfile(self, event): - - container = self.data.r0 container.obs_id = -1 @@ -262,7 +260,6 @@ def __init__(self, file_list): self._camera_config = {} self.camera_config = None - paths = [] for file_name in file_list: paths.append(file_name) @@ -286,12 +283,11 @@ def __init__(self, file_list): if(self.camera_config is None): self.camera_config = self._camera_config[path] - except StopIteration: pass # verify that somewhere the CameraConfing is present - assert (self.camera_config) + assert self.camera_config def __iter__(self): return self @@ -325,5 +321,9 @@ def __len__(self): ) return total_length + def rewind(self): + for name, file in self._file.items(): + file.Events.protobuf_i_fits.rewind() + def num_inputs(self): return len(self._file) diff --git a/ctapipe/io/nectarcameventsource.py b/ctapipe/io/nectarcameventsource.py index b57078f782b..c6cfd9b9452 100644 --- a/ctapipe/io/nectarcameventsource.py +++ b/ctapipe/io/nectarcameventsource.py @@ -5,11 +5,15 @@ Needs protozfits v1.4.2 from github.com/cta-sst-1m/protozfitsreader """ -import numpy as np import glob +import numpy as np from astropy import units as u -from ctapipe.instrument import TelescopeDescription, SubarrayDescription, \ - CameraGeometry, OpticsDescription +from ctapipe.instrument import ( + TelescopeDescription, + SubarrayDescription, + CameraGeometry, + OpticsDescription, +) from .eventsource import EventSource from .lsteventsource import MultiFiles from .containers import NectarCAMDataContainer @@ -23,8 +27,6 @@ class NectarCAMEventSource(EventSource): """ def __init__(self, config=None, tool=None, **kwargs): - - """ Constructor Parameters @@ -46,7 +48,6 @@ def __init__(self, config=None, tool=None, **kwargs): # To overcome this we substitute the input_url with first file matching # the specified file mask (copied from MAGICEventSourceROOT). - if 'input_url' in kwargs.keys(): self.file_list = glob.glob(kwargs['input_url']) self.file_list.sort() @@ -59,8 +60,9 @@ def __init__(self, config=None, tool=None, **kwargs): self.multi_file = MultiFiles(self.file_list) self.camera_config = self.multi_file.camera_config - self.log.info("Read {} input files".format(self.multi_file.num_inputs())) - + self.log.info("Read {} input files".format( + self.multi_file.num_inputs() + )) def _generator(self): @@ -71,7 +73,6 @@ def _generator(self): # fill data from the CameraConfig table self.fill_nectarcam_service_container_from_zfile() - # Instrument information for tel_id in self.data.nectarcam.tels_with_data: assert (tel_id == 0) # only one telescope for the moment (id = 0) @@ -94,15 +95,12 @@ def _generator(self): # LSTs telescope position tel_pos = {tel_id: [0., 0., 0] * u.m} - self.subarray = SubarrayDescription("MST prototype subarray") self.subarray.tels = tels self.subarray.positions = tel_pos self.data.inst.subarray = self.subarray - - # loop on events for count, event in enumerate(self.multi_file): @@ -115,7 +113,6 @@ def _generator(self): self.fill_r0_container_from_zfile(event) yield self.data - @staticmethod def is_compatible(file_path): from .sst1meventsource import is_fits_in_header @@ -155,6 +152,7 @@ def fill_nectarcam_service_container_from_zfile(self): self.data.nectarcam.tels_with_data = [self.camera_config.telescope_id, ] svc_container = self.data.nectarcam.tel[self.camera_config.telescope_id].svc + svc_container.telescope_id = self.camera_config.telescope_id svc_container.cs_serial = self.camera_config.cs_serial svc_container.configuration_id = self.camera_config.configuration_id @@ -170,14 +168,10 @@ def fill_nectarcam_service_container_from_zfile(self): svc_container.idaq_version = self.camera_config.nectarcam.idaq_version svc_container.cdhs_version = self.camera_config.nectarcam.cdhs_version svc_container.algorithms = self.camera_config.nectarcam.algorithms - # svc_container.pre_proc_algorithms = camera_config.nectarcam.pre_proc_algorithms - - def fill_nectarcam_event_container_from_zfile(self, event): - event_container = self.data.nectarcam.tel[self.camera_config.telescope_id].evt event_container.configuration_id = event.configuration_id @@ -195,7 +189,6 @@ def fill_nectarcam_event_container_from_zfile(self, event): def fill_r0_camera_container_from_zfile(self, container, event): - container.num_samples = self.camera_config.num_samples container.trigger_time = event.trigger_time_s container.trigger_type = event.trigger_type @@ -213,11 +206,17 @@ def fill_r0_camera_container_from_zfile(self, container, event): .format(event.waveform.shape[0])) + container.pixel_status = np.zeros([self.n_camera_pixels]) + container.pixel_status[self.camera_config.expected_pixels_id] = \ + event.pixel_status + reshaped_waveform = np.array( event.waveform - ).reshape(n_gains, - self.camera_config.num_pixels, - container.num_samples) + ).reshape( + n_gains, + self.camera_config.num_pixels, + container.num_samples + ) # initialize the waveform container to zero container.waveform = np.zeros([n_gains, @@ -228,7 +227,6 @@ def fill_r0_camera_container_from_zfile(self, container, event): container.waveform[:, self.camera_config.expected_pixels_id, :] \ = reshaped_waveform - def fill_r0_container_from_zfile(self, event): container = self.data.r0 container.obs_id = -1 @@ -241,4 +239,3 @@ def fill_r0_container_from_zfile(self, event): r0_camera_container, event ) - diff --git a/examples/calc_flatfield.py b/examples/calc_flatfield.py new file mode 100644 index 00000000000..8cfca81d7e3 --- /dev/null +++ b/examples/calc_flatfield.py @@ -0,0 +1,103 @@ +""" +Extract flat field coefficients from flasher data files. +""" +from traitlets import Dict, List, Unicode + +from ctapipe.core import Provenance +from ctapipe.io import HDF5TableWriter +from ctapipe.core import Tool +from ctapipe.io import EventSourceFactory + +from ctapipe.image import ChargeExtractorFactory, WaveformCleanerFactory +from ctapipe.calib.camera.flatfield import FlatFieldFactory +from ctapipe.io.containers import MonitorDataContainer + + +class FlatFieldHDF5Writer(Tool): + name = "FlatFieldHDF5Writer" + description = "Generate a HDF5 file with flat field coefficients" + + output_file = Unicode( + 'flatfield.hdf5', + help='Name of the output flat field file file' + ).tag(config=True) + + aliases = Dict(dict( + input_file='EventSourceFactory.input_url', + max_events='EventSourceFactory.max_events', + allowed_tels='EventSourceFactory.allowed_tels', + charge_extractor='ChargeExtractorFactory.product', + window_width='ChargeExtractorFactory.window_width', + t0='ChargeExtractorFactory.t0', + window_shift='ChargeExtractorFactory.window_shift', + sig_amp_cut_HG='ChargeExtractorFactory.sig_amp_cut_HG', + sig_amp_cut_LG='ChargeExtractorFactory.sig_amp_cut_LG', + lwt='ChargeExtractorFactory.lwt', + cleaner='WaveformCleanerFactory.product', + cleaner_width='WaveformCleanerFactory.baseline_width', + cleaner_shift='WaveformCleanerFactory.window_shift', + generator='FlatFieldFactory.product', + tel_id='FlatFieldFactory.tel_id', + sample_duration='FlatFieldFactory.sample_duration', + sample_size='FlatFieldFactory.sample_size', + n_channels='FlatFieldFactory.n_channels', + )) + + classes = List([EventSourceFactory, + ChargeExtractorFactory, + WaveformCleanerFactory, + FlatFieldFactory, + MonitorDataContainer, + HDF5TableWriter + ]) + + def __init__(self, **kwargs): + super().__init__(**kwargs) + self.eventsource = None + self.flatfield = None + self.container = None + self.writer = None + + def setup(self): + kwargs = dict(config=self.config, tool=self) + self.eventsource = EventSourceFactory.produce(**kwargs) + self.flatfield = FlatFieldFactory.produce(**kwargs) + + self.container = MonitorDataContainer() + self.writer = HDF5TableWriter( + filename=self.output_file, group_name='flatfield', overwrite=True + ) + + def start(self): + '''Flat field coefficient calculator''' + + # initialize the flat fieLd containers + self.container.flatfield.tels_with_data.append(self.flatfield.tel_id) + + for count, event in enumerate(self.eventsource): + + ff_data = self.flatfield.calculate_relative_gain(event) + + if ff_data: + self.container.flatfield.tel[self.flatfield.tel_id] = ff_data + table_name = 'tel_' + str(self.flatfield.tel_id) + + self.log.info("write event in table: /flatfield/%s", + table_name) + self.writer.write(table_name, ff_data) + + def finish(self): + Provenance().add_output_file( + self.output_file, + role='mon.tel.flatfield' + ) + self.writer.close() + + +def main(): + exe = FlatFieldHDF5Writer() + exe.run() + + +if __name__ == '__main__': + main()