diff --git a/examples/simple_event_writer.json b/examples/simple_event_writer.json new file mode 100644 index 00000000000..3dd4d7f7334 --- /dev/null +++ b/examples/simple_event_writer.json @@ -0,0 +1,9 @@ +{ + "Analysis": { + "allowed_tels": [1, 2, 3, 4] + }, + "Preselect": { + "image_amplitude": {"min": 50}, + "n_pixel": {"min": 3} + } +} \ No newline at end of file diff --git a/examples/simple_event_writer.py b/examples/simple_event_writer.py new file mode 100755 index 00000000000..51975d3ac22 --- /dev/null +++ b/examples/simple_event_writer.py @@ -0,0 +1,115 @@ +#!/usr/bin/env python3 +""" +This example uses a configuration file in JSON format to +process the events and apply pre-selection cuts to the images +(charge and number of pixels). +An HDF5 file is written with image MC and moment parameters +(e.g. length, width, image amplitude, etc.). +""" + +import numpy as np +from tqdm import tqdm + +from ctapipe.core import Tool +from ctapipe.core.traits import Unicode, List, Dict, Bool +from ctapipe.io import EventSourceFactory, HDF5TableWriter +from ctapipe.calib import CameraCalibrator +from ctapipe.utils.CutFlow import CutFlow +from ctapipe.image import hillas_parameters, tailcuts_clean + + +class SimpleEventWriter(Tool): + name = 'ctapipe-simple-event-writer' + description = Unicode(__doc__) + + infile = Unicode(help='input file to read', default='').tag(config=True) + outfile = Unicode(help='output file name', default_value='output.h5').tag(config=True) + progress = Bool(help='display progress bar', default_value=True).tag(config=True) + + aliases = Dict({ + 'infile': 'EventSourceFactory.input_url', + 'outfile': 'SimpleEventWriter.outfile', + 'max-events': 'EventSourceFactory.max_events', + 'progress': 'SimpleEventWriter.progress' + }) + classes = List([EventSourceFactory, CameraCalibrator, CutFlow]) + + def setup(self): + self.log.info('Configure EventSourceFactory...') + + self.event_source = EventSourceFactory.produce( + config=self.config, tool=self, product='HESSIOEventSource' + ) + self.event_source.allowed_tels = self.config['Analysis']['allowed_tels'] + + self.calibrator = CameraCalibrator( + config=self.config, tool=self, eventsource=self.event_source + ) + + self.writer = HDF5TableWriter( + filename=self.outfile, group_name='image_infos', overwrite=True + ) + + # Define Pre-selection for images + preselcuts = self.config['Preselect'] + self.image_cutflow = CutFlow('Image preselection') + self.image_cutflow.set_cuts(dict( + no_sel=None, + n_pixel=lambda s: np.count_nonzero(s) < preselcuts['n_pixel']['min'], + image_amplitude=lambda q: q < preselcuts['image_amplitude']['min'] + )) + + # Define Pre-selection for events + self.event_cutflow = CutFlow('Event preselection') + self.event_cutflow.set_cuts(dict( + no_sel=None + )) + + def start(self): + self.log.info('Loop on events...') + + for event in tqdm( + self.event_source, + desc='EventWriter', + total=self.event_source.max_events, + disable=~self.progress): + + self.event_cutflow.count('no_sel') + self.calibrator.calibrate(event) + + for tel_id in event.dl0.tels_with_data: + self.image_cutflow.count('no_sel') + + camera = event.inst.subarray.tel[tel_id].camera + dl1_tel = event.dl1.tel[tel_id] + + # Image cleaning + image = dl1_tel.image[0] # Waiting for automatic gain selection + mask = tailcuts_clean(camera, image, picture_thresh=10, boundary_thresh=5) + cleaned = image.copy() + cleaned[~mask] = 0 + + # Preselection cuts + if self.image_cutflow.cut('n_pixel', cleaned): + continue + if self.image_cutflow.cut('image_amplitude', np.sum(cleaned)): + continue + + # Image parametrisation + params = hillas_parameters(camera, cleaned, container=True) + + # Save Ids, MC infos and Hillas informations + self.writer.write(camera.cam_id, [event.r0, event.mc, params]) + + def finish(self): + self.log.info('End of job.') + + self.image_cutflow() + self.event_cutflow() + self.writer.close() + + +if __name__ == '__main__': + tool = SimpleEventWriter() + tool.run() +