Merge pull request #252 from cta-observatory/calib_true_pe

Implement reading of photoelectrons for calibration events, fixes #250

.github/workflows/ci.yml

@@ -30,7 +30,6 @@ jobs:
           python --version
           pip install codecov pytest-cov
           pip install -U pytest docutils pygments Cython matplotlib
-          pip install https://github.com/cta-observatory/pyhessio/archive/v2.1.1.tar.gz
           pip install -e .
           pip freeze
 

.gitignore

@@ -1,4 +1,4 @@
-eventio/_version.py
+src/eventio/_version.py
 # Cython extensions
 header.c
 var_int.c

src/eventio/simtel/objects.py

@@ -1797,6 +1797,13 @@ def parse(self):
         return data
 
 
+class CalibrationPhotoelectrons(EventIOObject):
+    eventio_type = 2034
+
+    def __str__(self):
+        return f'{self.__class__.__name__}(type={self.header.id})'
+
+
 
 def merge_structured_arrays_into_dict(arrays):
     result = dict()

src/eventio/simtel/simtelfile.py

@@ -42,6 +42,7 @@
     TelescopeEventHeader,
     TrackingPosition,
     TriggerInformation,
+    CalibrationPhotoelectrons,
 )
 
 
@@ -112,6 +113,7 @@ def __init__(self, path, allowed_telescopes=None, skip_calibration=False, zcat=T
         self.current_array_event = None
         self.current_calibration_event = None
         self.current_calibration_event_id = None
+        self.current_calibration_pe = {}
         self.skip_calibration = skip_calibration
 
         # read the header:
@@ -208,6 +210,26 @@ def next_low_level(self):
                 # duplicated event_ids
                 self.current_calibration_event_id = -array_event.header.id
                 self.current_calibration_event['calibration_type'] = o.type
+        elif isinstance(o, CalibrationPhotoelectrons):
+            telescope_data = next(o)
+            if not isinstance(telescope_data, iact.TelescopeData):
+                warnings.warn(
+                    f"Unexpected sub-object: {telescope_data} in {o}, ignoring",
+                    UnknownObjectWarning
+                )
+                return
+
+            self.current_calibration_pe = {}
+            for photoelectrons in telescope_data:
+                if not isinstance(photoelectrons, iact.PhotoElectrons):
+                    warnings.warn(
+                        f"Unexpected sub-object: {photoelectrons} in {telescope_data}, ignoring",
+                        UnknownObjectWarning
+                    )
+
+                tel_id = photoelectrons.telescope_id
+                self.current_calibration_pe[tel_id] = photoelectrons.parse()
+
 
         elif isinstance(o, History):
             for sub in o:
@@ -344,6 +366,7 @@ def try_build_event(self):
                 'tracking_positions': event['tracking_positions'],
                 'trigger_information': event['trigger_information'],
                 'calibration_type': event['calibration_type'],
+                'photoelectrons': self.current_calibration_pe,
             }
 
             event_data['camera_monitorings'] = {

tests/simtel/test_simtel_objects.py

@@ -788,3 +788,20 @@ def test_2033():
                 assert (data['fadc_amp_hg'][valid] > data['fadc_amp_lg'][valid]).all()
 
         assert tel_id == 79
+
+
+def test_2034():
+    from eventio.simtel.objects import CalibrationPhotoelectrons
+    from eventio.iact.objects import TelescopeData, PhotoElectrons
+
+    with EventIOFile('tests/resources/calib_true_pe.simtel.zst') as f:
+        # file contains two calibration events
+        objects = yield_n_and_assert(f, CalibrationPhotoelectrons, n=2)
+
+        for object in objects:
+            # CalibrationEventPhotoelectrons should be a container object
+            # holding a corsika array with photo electrons
+            telescope_data = next(object)
+            assert isinstance(telescope_data, TelescopeData)
+            photo_electrons = next(telescope_data)
+            assert isinstance(photo_electrons, PhotoElectrons)

tests/simtel/test_simtelfile.py

@@ -1,5 +1,6 @@
 from pytest import importorskip
 from eventio.simtel import SimTelFile
+import numpy as np
 
 prod2_path = 'tests/resources/gamma_test.simtel.gz'
 prod3_path = 'tests/resources/gamma_test_large_truncated.simtel.gz'
@@ -242,6 +243,14 @@ def test_missing_photons():
         assert e['emitter'] == {}
 
 
+def test_calibration_photoelectrons():
+    with SimTelFile('tests/resources/calib_true_pe.simtel.zst') as f:
+        for e, expected_pe in zip(f, (1.7, 20)):
+            assert 'photoelectrons' in e
+            assert 0 in e['photoelectrons']
+            true_image = e['photoelectrons'][0]['photoelectrons']
+            assert np.isclose(np.mean(true_image), expected_pe, rtol=0.05)
+
 def test_history_meta():
     with SimTelFile(history_meta_path) as f:
         assert isinstance(f.global_meta, dict)