Use eventio v0.17.1 (#966)

* start using eventio v0.17.0

* Use is_eventio

* Implement reading calibration events

* Call it skip_calibration_events

* Use pip for now for extra

* Use eventio 0.17.1

* Update version of eventio also in setup.py

ctapipe/io/containers.py

@@ -428,6 +428,7 @@ class TelescopePointingContainer(Container):
 class DataContainer(Container):
     """ Top-level container for all event information """
 
+    event_type = Field('data', "Event type")
     r0 = Field(R0Container(), "Raw Data")
     r1 = Field(R1Container(), "R1 Calibrated Data")
     dl0 = Field(DL0Container(), "DL0 Data Volume Reduced Data")

ctapipe/io/simteleventsource.py

@@ -6,20 +6,25 @@
 from astropy.coordinates import Angle
 from astropy.time import Time
 from ctapipe.instrument import TelescopeDescription, SubarrayDescription
-import struct
-import gzip
+from traitlets import Bool
+
 
 from eventio.simtel.simtelfile import SimTelFile
+from eventio.file_types import is_eventio
 
 __all__ = ['SimTelEventSource']
 
 
 class SimTelEventSource(EventSource):
+    skip_calibration_events = Bool(True, help='Skip calibration events').tag(config=True)
 
     def __init__(self, config=None, tool=None, **kwargs):
         super().__init__(config=config, tool=tool, **kwargs)
         self.metadata['is_simulation'] = True
-        self.file_ = SimTelFile(self.input_url)
+        self.file_ = SimTelFile(
+            self.input_url,
+            skip_calibration=self.skip_calibration_events
+        )
 
         self._subarray_info = self.prepare_subarray_info(
             self.file_.telescope_descriptions,
@@ -74,18 +79,7 @@ def prepare_subarray_info(telescope_descriptions, header):
 
     @staticmethod
     def is_compatible(file_path):
-        # read the first 4 bytes
-        with open(file_path, 'rb') as f:
-            marker_bytes = f.read(4)
-
-        # if file is gzip, read the first 4 bytes with gzip again
-        if marker_bytes[0] == 0x1f and marker_bytes[1] == 0x8b:
-            with gzip.open(file_path, 'rb') as f:
-                marker_bytes = f.read(4)
-
-        # check for the simtel magic marker
-        int_marker, = struct.unpack('I', marker_bytes)
-        return int_marker == 3558836791 or int_marker == 931798996
+        return is_eventio(file_path)
 
     def _generator(self):
         try:
@@ -106,8 +100,14 @@ def __generator(self):
         for counter, array_event in enumerate(self.file_):
             # next lines are just for debugging
             self.array_event = array_event
+            data.event_type = array_event['type']
+
+            # calibration events do not have an event id
+            if data.event_type == 'calibration':
+                event_id = -1
+            else:
+                event_id = array_event['event_id']
 
-            event_id = array_event['event_id']
             data.inst.subarray = self._subarray_info
 
             obs_id = self.file_.header['run']
@@ -134,67 +134,14 @@ def __generator(self):
                 data.dl0.tels_with_data = selected
 
             trigger_information = array_event['trigger_information']
-            mc_event = array_event['mc_event']
-            mc_shower = array_event['mc_shower']
 
             data.trig.tels_with_trigger = trigger_information['triggered_telescopes']
             time_s, time_ns = trigger_information['gps_time']
             data.trig.gps_time = Time(time_s * u.s, time_ns * u.ns,
                                       format='unix', scale='utc')
 
-            data.mc.energy = mc_shower['energy'] * u.TeV
-            data.mc.alt = Angle(mc_shower['altitude'], u.rad)
-            data.mc.az = Angle(mc_shower['azimuth'], u.rad)
-            data.mc.core_x = mc_event['xcore'] * u.m
-            data.mc.core_y = mc_event['ycore'] * u.m
-            first_int = mc_shower['h_first_int'] * u.m
-            data.mc.h_first_int = first_int
-            data.mc.x_max = mc_shower['xmax'] * u.g / (u.cm**2)
-            data.mc.shower_primary_id = mc_shower['primary_id']
-
-            # mc run header data
-            data.mcheader.run_array_direction = Angle(
-                self.file_.header['direction'] * u.rad
-            )
-            mc_run_head = self.file_.mc_run_headers[-1]
-            data.mcheader.corsika_version = mc_run_head['shower_prog_vers']
-            data.mcheader.simtel_version = mc_run_head['detector_prog_vers']
-            data.mcheader.energy_range_min = mc_run_head['E_range'][0] * u.TeV
-            data.mcheader.energy_range_max = mc_run_head['E_range'][1] * u.TeV
-            data.mcheader.prod_site_B_total = mc_run_head['B_total'] * u.uT
-            data.mcheader.prod_site_B_declination = Angle(
-                mc_run_head['B_declination'] * u.rad)
-            data.mcheader.prod_site_B_inclination = Angle(
-                mc_run_head['B_inclination'] * u.rad)
-            data.mcheader.prod_site_alt = mc_run_head['obsheight'] * u.m
-            data.mcheader.spectral_index = mc_run_head['spectral_index']
-            data.mcheader.shower_prog_start = mc_run_head['shower_prog_start']
-            data.mcheader.shower_prog_id = mc_run_head['shower_prog_id']
-            data.mcheader.detector_prog_start = mc_run_head['detector_prog_start']
-            data.mcheader.detector_prog_id = mc_run_head['detector_prog_id']
-            data.mcheader.num_showers = mc_run_head['num_showers']
-            data.mcheader.shower_reuse = mc_run_head['num_use']
-            data.mcheader.max_alt = mc_run_head['alt_range'][1] * u.rad
-            data.mcheader.min_alt = mc_run_head['alt_range'][0] * u.rad
-            data.mcheader.max_az = mc_run_head['az_range'][1] * u.rad
-            data.mcheader.min_az = mc_run_head['az_range'][0] * u.rad
-            data.mcheader.diffuse = mc_run_head['diffuse']
-            data.mcheader.max_viewcone_radius = mc_run_head['viewcone'][1] * u.deg
-            data.mcheader.min_viewcone_radius = mc_run_head['viewcone'][0] * u.deg
-            data.mcheader.max_scatter_range = mc_run_head['core_range'][1] * u.m
-            data.mcheader.min_scatter_range = mc_run_head['core_range'][0] * u.m
-            data.mcheader.core_pos_mode = mc_run_head['core_pos_mode']
-            data.mcheader.injection_height = mc_run_head['injection_height'] * u.m
-            data.mcheader.atmosphere = mc_run_head['atmosphere']
-            data.mcheader.corsika_iact_options = mc_run_head['corsika_iact_options']
-            data.mcheader.corsika_low_E_model = mc_run_head['corsika_low_E_model']
-            data.mcheader.corsika_high_E_model = mc_run_head['corsika_high_E_model']
-            data.mcheader.corsika_bunchsize = mc_run_head['corsika_bunchsize']
-            data.mcheader.corsika_wlen_min = mc_run_head['corsika_wlen_min'] * u.nm
-            data.mcheader.corsika_wlen_max = mc_run_head['corsika_wlen_max'] * u.nm
-            data.mcheader.corsika_low_E_detail = mc_run_head['corsika_low_E_detail']
-            data.mcheader.corsika_high_E_detail = mc_run_head['corsika_high_E_detail']
-
+            if data.event_type == 'data':
+                self.fill_mc_information(data, array_event)
 
             # this should be done in a nicer way to not re-allocate the
             # data each time (right now it's just deleted and garbage
@@ -228,7 +175,7 @@ def __generator(self):
                     data.r0.tel[tel_id].trig_pix_id = pixel_lists[0]['pixel_list']
 
                 pixel_settings = telescope_description['pixel_settings']
-                data.mc.tel[tel_id].reference_pulse_shape = pixel_settings['refshape'].astype('float64')
+                data.mc.tel[tel_id].reference_pulse_shape = pixel_settings['ref_shape'].astype('float64')
                 data.mc.tel[tel_id].meta['refstep'] = float(pixel_settings['ref_step'])
                 data.mc.tel[tel_id].time_slice = float(pixel_settings['time_slice'])
 
@@ -246,5 +193,60 @@ def __generator(self):
                 data.mc.tel[tel_id].altitude_cor = tracking_position.get('altitude_cor', 0)
             yield data
 
-
+    def fill_mc_information(self, data, array_event):
+        mc_event = array_event['mc_event']
+        mc_shower = array_event['mc_shower']
+
+        data.mc.energy = mc_shower['energy'] * u.TeV
+        data.mc.alt = Angle(mc_shower['altitude'], u.rad)
+        data.mc.az = Angle(mc_shower['azimuth'], u.rad)
+        data.mc.core_x = mc_event['xcore'] * u.m
+        data.mc.core_y = mc_event['ycore'] * u.m
+        first_int = mc_shower['h_first_int'] * u.m
+        data.mc.h_first_int = first_int
+        data.mc.x_max = mc_shower['xmax'] * u.g / (u.cm**2)
+        data.mc.shower_primary_id = mc_shower['primary_id']
+
+        # mc run header data
+        data.mcheader.run_array_direction = Angle(
+            self.file_.header['direction'] * u.rad
+        )
+        mc_run_head = self.file_.mc_run_headers[-1]
+        data.mcheader.corsika_version = mc_run_head['shower_prog_vers']
+        data.mcheader.simtel_version = mc_run_head['detector_prog_vers']
+        data.mcheader.energy_range_min = mc_run_head['E_range'][0] * u.TeV
+        data.mcheader.energy_range_max = mc_run_head['E_range'][1] * u.TeV
+        data.mcheader.prod_site_B_total = mc_run_head['B_total'] * u.uT
+        data.mcheader.prod_site_B_declination = Angle(
+            mc_run_head['B_declination'] * u.rad)
+        data.mcheader.prod_site_B_inclination = Angle(
+            mc_run_head['B_inclination'] * u.rad)
+        data.mcheader.prod_site_alt = mc_run_head['obsheight'] * u.m
+        data.mcheader.spectral_index = mc_run_head['spectral_index']
+        data.mcheader.shower_prog_start = mc_run_head['shower_prog_start']
+        data.mcheader.shower_prog_id = mc_run_head['shower_prog_id']
+        data.mcheader.detector_prog_start = mc_run_head['detector_prog_start']
+        data.mcheader.detector_prog_id = mc_run_head['detector_prog_id']
+        data.mcheader.num_showers = mc_run_head['n_showers']
+        data.mcheader.shower_reuse = mc_run_head['n_use']
+        data.mcheader.max_alt = mc_run_head['alt_range'][1] * u.rad
+        data.mcheader.min_alt = mc_run_head['alt_range'][0] * u.rad
+        data.mcheader.max_az = mc_run_head['az_range'][1] * u.rad
+        data.mcheader.min_az = mc_run_head['az_range'][0] * u.rad
+        data.mcheader.diffuse = mc_run_head['diffuse']
+        data.mcheader.max_viewcone_radius = mc_run_head['viewcone'][1] * u.deg
+        data.mcheader.min_viewcone_radius = mc_run_head['viewcone'][0] * u.deg
+        data.mcheader.max_scatter_range = mc_run_head['core_range'][1] * u.m
+        data.mcheader.min_scatter_range = mc_run_head['core_range'][0] * u.m
+        data.mcheader.core_pos_mode = mc_run_head['core_pos_mode']
+        data.mcheader.injection_height = mc_run_head['injection_height'] * u.m
+        data.mcheader.atmosphere = mc_run_head['atmosphere']
+        data.mcheader.corsika_iact_options = mc_run_head['corsika_iact_options']
+        data.mcheader.corsika_low_E_model = mc_run_head['corsika_low_E_model']
+        data.mcheader.corsika_high_E_model = mc_run_head['corsika_high_E_model']
+        data.mcheader.corsika_bunchsize = mc_run_head['corsika_bunchsize']
+        data.mcheader.corsika_wlen_min = mc_run_head['corsika_wlen_min'] * u.nm
+        data.mcheader.corsika_wlen_max = mc_run_head['corsika_wlen_max'] * u.nm
+        data.mcheader.corsika_low_E_detail = mc_run_head['corsika_low_E_detail']
+        data.mcheader.corsika_high_E_detail = mc_run_head['corsika_high_E_detail']
 

ctapipe/io/tests/test_simteleventsource.py

@@ -8,6 +8,7 @@
 
 gamma_test_large_path = get_dataset_path("gamma_test_large.simtel.gz")
 gamma_test_path = get_dataset_path("gamma_test.simtel.gz")
+calib_events_path = get_dataset_path('calib_events.simtel.gz')
 
 
 def compare_sources(input_url):
@@ -161,3 +162,44 @@ def test_additional_meta_data_from_mc_header():
             value.to_value(expectation.unit),
             expectation.to_value(expectation.unit)
         )
+
+
+def test_hessio_file_reader():
+    dataset = gamma_test_path
+    with SimTelEventSource(input_url=dataset) as reader:
+        assert reader.is_compatible(dataset)
+        assert not reader.is_stream
+        for event in reader:
+            if event.count == 0:
+                assert event.r0.tels_with_data == {38, 47}
+            elif event.count == 1:
+                assert event.r0.tels_with_data == {11, 21, 24, 26, 61, 63, 118,
+                                                   119}
+            else:
+                break
+        for event in reader:
+            # Check generator has restarted from beginning
+            assert event.count == 0
+            break
+
+    # test that max_events works:
+    max_events = 5
+    with SimTelEventSource(input_url=dataset, max_events=max_events) as reader:
+        count = 0
+        for _ in reader:
+            count += 1
+        assert count == max_events
+
+    # test that the allowed_tels mask works:
+    with SimTelEventSource(input_url=dataset, allowed_tels={3, 4}) as reader:
+        for event in reader:
+            assert event.r0.tels_with_data.issubset(reader.allowed_tels)
+
+
+def test_calibration_events():
+    with SimTelEventSource(
+            input_url=calib_events_path,
+            skip_calibration_events=False,
+    ) as reader:
+        for e in reader:
+            pass

environment.yml

@@ -4,7 +4,6 @@ channels:
   - default
   - cta-observatory
 dependencies:
-  - ctapipe-extra=0.2.16
   - zeromq
   - pyzmq>=17  # needed for correct function of notebooks on OSX
   - astropy
@@ -42,6 +41,7 @@ dependencies:
   - conda-forge::nbsphinx
   - pip:
       - pytest_runner
-      - eventio==0.11.0
+      - eventio==0.17.1
+      - https://github.com/cta-observatory/ctapipe-extra/archive/v0.2.17.tar.gz
 
 

py3.6_env.yaml

@@ -3,7 +3,6 @@ channels:
   - default
   - cta-observatory
 dependencies:
-  - ctapipe-extra=0.2.16
   - python=3.6 # nail the python version, so conda does not try upgrading / dowgrading
   - zeromq
   - pyzmq>=17  # needed for correct function of notebooks on OSX
@@ -42,5 +41,6 @@ dependencies:
   - conda-forge::nbsphinx
   - pip:
       - pytest_runner
-      - eventio==0.11.0
+      - eventio==0.17.1
       - https://github.com/cta-observatory/pyhessio/archive/v2.1.1.tar.gz
+      - https://github.com/cta-observatory/ctapipe-extra/archive/v0.2.17.tar.gz

py3.7_env.yaml

@@ -3,7 +3,6 @@ channels:
   - default
   - cta-observatory
 dependencies:
-  - ctapipe-extra=0.2.16
   - python=3.7
   - zeromq
   - pyzmq>=17  # needed for correct function of notebooks on OSX
@@ -42,5 +41,6 @@ dependencies:
   - conda-forge::nbsphinx
   - pip:
       - pytest_runner
-      - eventio==0.11.0
+      - eventio==0.17.1
       - https://github.com/cta-observatory/pyhessio/archive/v2.1.1.tar.gz
+      - https://github.com/cta-observatory/ctapipe-extra/archive/v0.2.17.tar.gz

setup.py

@@ -65,7 +65,7 @@
           'pandas',
           'bokeh>=1.0.1',
           'scikit-learn',
-          'eventio==0.11.0',
+          'eventio==0.17.1',
       ],
       tests_require=[
           'pytest',