Implement type 2033, PixelMonitoring, fixes #243

eventio/simtel/__init__.py

@@ -27,6 +27,7 @@
     MCShower,
     PixelCalibration,
     PixelList,
+    PixelMonitoring,
     PixelSettings,
     PixelTiming,
     PixelTriggerTimes,
@@ -68,6 +69,7 @@
     'MCShower',
     'PixelCalibration',
     'PixelList',
+    'PixelMonitoring',
     'PixelSettings',
     'PixelTiming',
     'PixelTriggerTimes',

eventio/simtel/objects.py

@@ -1759,6 +1759,39 @@ def parse(self):
         return trigger_times
 
 
+class PixelMonitoring(TelescopeObject):
+    eventio_type = 2033
+
+    def parse(self):
+        assert_exact_version(self, supported_version=0)
+        byte_stream = BytesIO(self.read())
+
+        flags = read_int(byte_stream)
+        n_pixels = read_int(byte_stream)
+        n_gains = read_short(byte_stream)
+
+        data = {"flags": flags, "n_pixels": n_pixels, "n_gains": n_gains}
+        if flags & 0b0000_0001 != 0:
+            data['nsb_rate'] = read_array(byte_stream, np.float32, n_pixels)
+        if flags & 0b0000_0010 != 0:
+            data['qe_rel'] = read_array(byte_stream, np.float32, n_pixels)
+        if flags & 0b0000_0100 != 0:
+            data['gain_rel'] = read_array(byte_stream, np.float32, n_pixels)
+        if flags & 0b0000_1000 != 0:
+            data['hv_rel'] = read_array(byte_stream, np.float32, n_pixels)
+        if flags & 0b0001_0000 != 0:
+            data['current'] = read_array(byte_stream, np.float32, n_pixels)
+        if flags & 0b0010_0000 != 0:
+            data['fadc_amp_hg'] = read_array(byte_stream, np.float32, n_pixels)
+        if flags & 0b0100_0000 != 0 and n_gains > 1:
+            data['fadc_amp_lg'] = read_array(byte_stream, np.float32, n_pixels)
+        if flags & 0b1000_0000 != 0:
+            data['disabled'] = read_array(byte_stream, np.bool_, n_pixels)
+
+        return data
+
+
+
 def merge_structured_arrays_into_dict(arrays):
     result = dict()
     for array in arrays:

eventio/simtel/simtelfile.py

@@ -35,6 +35,7 @@
     PixelSettings,
     PixelTiming,
     PixelTriggerTimes,
+    PixelMonitoring,
     PointingCorrection,
     RunHeader,
     TelescopeEvent,
@@ -95,6 +96,7 @@ def __init__(self, path, allowed_telescopes=None, skip_calibration=False, zcat=T
         self.telescope_meta = {}
         self.global_meta = {}
         self.telescope_descriptions = defaultdict(dict)
+        self.pixel_monitorings = defaultdict(dict)
         self.camera_monitorings = defaultdict(dict)
         self.laser_calibrations = defaultdict(dict)
         self.current_mc_shower = None
@@ -172,6 +174,9 @@ def next_low_level(self):
         elif isinstance(o, LaserCalibration):
             self.laser_calibrations[o.telescope_id].update(o.parse())
 
+        elif isinstance(o, PixelMonitoring):
+            self.pixel_monitorings[o.telescope_id].update(o.parse())
+
         elif isinstance(o, telescope_descriptions_types):
             key = camel_to_snake(o.__class__.__name__)
             self.telescope_descriptions[o.telescope_id][key] = o.parse()
@@ -292,6 +297,11 @@ def try_build_event(self):
                 for telescope_id in self.current_array_event['telescope_events'].keys()
             }
 
+            event_data['pixel_monitorings'] = {
+                telescope_id: copy(self.pixel_monitorings[telescope_id])
+                for telescope_id in self.current_array_event['telescope_events'].keys()
+            }
+
             self.current_array_event = None
 
             return event_data

tests/simtel/test_simtel_objects.py

@@ -724,3 +724,39 @@ def test_2032():
             assert 'n_times' in d
             assert 'pixel_ids' in d
             assert 'trigger_times' in d
+
+
+def test_2033():
+    from eventio.simtel.objects import PixelMonitoring
+
+    with EventIOFile('tests/resources/type2033.simtel.zst') as f:
+        objects = yield_n_and_assert(f, PixelMonitoring, n=79)
+        tel_id = 0
+        for tel_id, o in enumerate(objects, start=1):
+            assert o.header.id == tel_id
+            data = o.parse()
+            if data['n_gains'] == 2:
+                assert data.keys() == {
+                    'flags', 'n_pixels', 'n_gains', 'nsb_rate', 'qe_rel',
+                    'gain_rel', 'hv_rel', 'current', 'fadc_amp_hg',
+                    'fadc_amp_lg', 'disabled',
+                }
+            else:
+                assert data.keys() == {
+                    'flags', 'n_pixels', 'n_gains', 'nsb_rate', 'qe_rel',
+                    'gain_rel', 'hv_rel', 'current', 'fadc_amp_hg', 'disabled',
+                }
+
+
+            assert (data['qe_rel'] > 0).any()
+            # file was simulated with random disabled pixels
+            assert np.count_nonzero(data['disabled']) > 0
+            # disabled pixels should have no qe
+            assert (data['qe_rel'][data['disabled']] == 0.0).all()
+            # amplification in hg should be higher than lg
+            valid = ~data['disabled']
+
+            if data['n_gains'] == 2:
+                assert (data['fadc_amp_hg'][valid] > data['fadc_amp_lg'][valid]).all()
+
+        assert tel_id == 79

tests/simtel/test_simtelfile.py

@@ -12,6 +12,7 @@
 calib_path = 'tests/resources/calib_events.simtel.gz'
 frankenstein_path = 'tests/resources/gamma_merged.simtel.gz'
 history_meta_path = 'tests/resources/history_meta_75.simtel.zst'
+pixel_monitoring_path = 'tests/resources/type2033.simtel.zst'
 
 
 test_paths = [prod2_path, prod3_path, prod4_path]
@@ -73,6 +74,7 @@ def test_show_event_is_not_empty_and_has_some_members_for_sure():
                 'emitter',
                 'camera_monitorings',
                 'laser_calibrations',
+                'pixel_monitorings',
             }
 
             telescope_events = event['telescope_events']
@@ -237,3 +239,21 @@ def test_history_meta():
         assert isinstance(f.telescope_meta, dict)
         assert len(f.global_meta) == 11
         assert len(f.telescope_meta) == 19
+
+
+def test_type_2033():
+    with SimTelFile(pixel_monitoring_path) as f:
+        e = next(iter(f))
+        assert e['pixel_monitorings'].keys() == e['telescope_events'].keys()
+
+        for pixel_monitorings in e['pixel_monitorings'].values():
+            if pixel_monitorings['n_gains'] == 2:
+                assert pixel_monitorings.keys() == {
+                    'flags', 'n_pixels', 'n_gains', 'nsb_rate', 'qe_rel', 'gain_rel',
+                    'hv_rel', 'current', 'fadc_amp_hg', 'fadc_amp_lg', 'disabled',
+                }
+            else:
+                assert pixel_monitorings.keys() == {
+                    'flags', 'n_pixels', 'n_gains', 'nsb_rate', 'qe_rel', 'gain_rel',
+                    'hv_rel', 'current', 'fadc_amp_hg', 'disabled',
+                }