Merge pull request #153 from cta-observatory/simtel_2032

Simtel 2032, fixes #115

eventio/simtel/objects.py

@@ -1709,6 +1709,21 @@ class FSPhot(EventIOObject):
 class PixelTriggerTimes(TelescopeObject):
     eventio_type = 2032
 
+    def parse(self):
+        assert_exact_version(self, supported_version=0)
+        byte_stream = BytesIO(self.read())
+
+        trigger_times = {}
+        trigger_times['time_step'] = read_float(byte_stream)
+        n_times = read_utf8_like_signed_int(byte_stream)
+        trigger_times['n_times'] = n_times
+
+        data = byte_stream.read()
+        trigger_times['pixel_ids'], length = varint_array(data, n_times, offset=0)
+        trigger_times['trigger_times'], length = varint_array(data, n_times, offset=length)
+
+        return trigger_times
+
 
 def merge_structured_arrays_into_dict(arrays):
     result = dict()

eventio/simtel/simtelfile.py

@@ -19,7 +19,6 @@
     CameraOrganization,
     CameraSettings,
     CameraSoftwareSettings,
-    TriggerInformation,
     DisabledPixels,
     DriveSettings,
     History,
@@ -32,11 +31,13 @@
     PixelList,
     PixelSettings,
     PixelTiming,
+    PixelTriggerTimes,
     PointingCorrection,
     RunHeader,
     TelescopeEvent,
     TelescopeEventHeader,
     TrackingPosition,
+    TriggerInformation,
 )
 
 
@@ -361,4 +362,7 @@ def parse_telescope_event(telescope_event):
         elif isinstance(o, PixelList):
             event['pixel_lists'][o.code] = o.parse()
 
+        elif isinstance(o, PixelTriggerTimes):
+            event['pixel_trigger_times'] = o.parse()
+
     return event

examples/trigger_times.py

@@ -0,0 +1,62 @@
+# coding: utf-8
+from ctapipe.visualization import CameraDisplay
+from ctapipe.instrument import CameraGeometry
+from eventio import SimTelFile
+import astropy.units as u
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+
+f = SimTelFile('tests/resources/gamma_20deg_0deg_run103___cta-prod4-sst-astri_desert-2150m-Paranal-sst-astri.simtel.gz')
+cam = f.telescope_descriptions[1]['camera_settings']
+geom = CameraGeometry(
+    'astri',
+    np.arange(cam['n_pixels']),
+    cam['pixel_x'] * u.m,
+    cam['pixel_y'] * u.m,
+    cam['pixel_area']* u.m**2,
+    pix_type='rectangular',
+)
+
+
+it = iter(f)
+
+plt.ion()
+
+fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
+
+d1 = CameraDisplay(geom, ax=ax1)
+d2 = CameraDisplay(geom, ax=ax2)
+
+d1.add_colorbar(ax=ax1)
+d2.add_colorbar(ax=ax2)
+ax1.set_title('ADCSum')
+ax2.set_title('Pixel Trigger Times')
+
+x = geom.pix_x.to_value(u.m)
+y = geom.pix_y.to_value(u.m)
+for ax in (ax1, ax2):
+    ax.set_xlim(1.01 * x.min(), 1.01 * x.max())
+    ax.set_ylim(1.01 * y.min(), 1.01 * y.max())
+
+cmap = plt.get_cmap('viridis')
+cmap.set_bad('gray')
+d2.cmap = cmap
+
+fig.show()
+fig.tight_layout()
+
+
+for event in f:
+    for t in event['telescope_events'].values():
+
+        d1.image = t['adc_sums'][0]
+
+        trig = np.full(geom.n_pixels, np.nan)
+        trig[t['pixel_trigger_times']['pixel_ids']] = t['pixel_trigger_times']['trigger_times']
+        trig = np.ma.array(trig, mask=np.isnan(trig))
+
+        d2.image = trig
+
+        input('Enter for next telecope event')

tests/simtel/test_simtel_objects.py

@@ -638,6 +638,13 @@ def test_2031():
     assert False
 
 
-@pytest.mark.xfail
 def test_2032():
-    assert False
+    from eventio.simtel.objects import PixelTriggerTimes
+
+    with EventIOFile(prod4b_astri_file) as f:
+        for i, o in enumerate(yield_n_and_assert(f, PixelTriggerTimes, n=3)):
+            d = parse_and_assert_consumption(o, limit=2)
+            print(d)
+            assert 'n_times' in d
+            assert 'pixel_ids' in d
+            assert 'trigger_times' in d

tests/simtel/test_simtelfile.py

@@ -4,6 +4,7 @@
 prod2_path = 'tests/resources/gamma_test.simtel.gz'
 prod3_path = 'tests/resources/gamma_test_large_truncated.simtel.gz'
 prod4_path = 'tests/resources/gamma_20deg_0deg_run102___cta-prod4-sst-1m_desert-2150m-Paranal-sst-1m.simtel.gz'
+prod4_astri_path = 'tests/resources/gamma_20deg_0deg_run103___cta-prod4-sst-astri_desert-2150m-Paranal-sst-astri.simtel.gz'
 
 # using a zstd file ensures SimTelFile is not seeking back, when reading
 # a file
@@ -147,6 +148,14 @@ def test_allowed_tels():
     assert n_read == 3
 
 
+def test_pixel_trigger_times():
+    # astri files must have trigger times
+    with SimTelFile(prod4_astri_path) as f:
+        for counter, event in enumerate(f, start=1):
+            for telescope_event in event['telescope_events'].values():
+                assert 'pixel_trigger_times' in telescope_event
+
+
 def test_calibration_events():
     with SimTelFile(calib_path) as f:
         i = 0