Merge pull request #206 from cta-observatory/simtel_photons

Add support in SimTelFile for files with full corsika photon information
cta-observatory · May 22, 2020 · d8c1abc · d8c1abc
2 parents 8cdcf84 + 0863e57
commit d8c1abc
Show file tree

Hide file tree

Showing 4 changed files with 48 additions and 9 deletions.
diff --git a/eventio/__init__.py b/eventio/__init__.py
@@ -6,7 +6,7 @@
 from .histograms import Histograms
 
 
-__version__ = '1.2.0'
+__version__ = '1.3.0'
 
 __all__ = [
     'EventIOFile',

diff --git a/eventio/simtel/simtelfile.py b/eventio/simtel/simtelfile.py
@@ -92,6 +92,8 @@ def __init__(self, path, allowed_telescopes=None, skip_calibration=False, zcat=T
         self.current_mc_event = None
         self.current_photoelectron_sum = None
         self.current_photoelectrons = {}
+        self.current_photons = {}
+        self.current_emitter = {}
         self.current_array_event = None
         self.current_calibration_event = None
         self.skip_calibration = skip_calibration
@@ -145,7 +147,10 @@ def next_low_level(self):
             )
 
         elif isinstance(o, iact.TelescopeData):
-            self.current_photoelectrons = parse_photoelectrons(o)
+            photons, emitter, photoelectrons = parse_telescope_data(o)
+            self.current_photons = photons
+            self.current_emitter = emitter
+            self.current_photoelectrons = photoelectrons
 
         elif isinstance(o, CameraMonitoring):
             self.camera_monitorings[o.telescope_id].update(o.parse())
@@ -205,6 +210,9 @@ def try_build_mc_event(self):
                 'event_id': self.current_mc_event_id,
                 'mc_shower': self.current_mc_shower,
                 'mc_event': self.current_mc_event,
+                'photons': self.current_photons,
+                'emitter': self.current_emitter,
+                'photoelectrons': self.current_photoelectrons,
             }
             self.current_mc_event = None
             return event_data
@@ -242,8 +250,10 @@ def try_build_event(self):
                 'telescope_events': self.current_array_event['telescope_events'],
                 'tracking_positions': self.current_array_event['tracking_positions'],
                 'trigger_information': self.current_array_event['trigger_information'],
-                'photoelectron_sums': self.current_photoelectron_sum,
+                'photons': self.current_photons,
+                'emitter': self.current_emitter,
                 'photoelectrons': self.current_photoelectrons,
+                'photoelectron_sums': self.current_photoelectron_sum,
             }
 
             event_data['camera_monitorings'] = {
@@ -349,15 +359,23 @@ def parse_array_event(array_event, allowed_telescopes=None):
     }
 
 
-def parse_photoelectrons(telescope_data):
+def parse_telescope_data(telescope_data):
+    ''' Parse the TelescopeData block with Cherenkov Photon information'''
     check_type(telescope_data, iact.TelescopeData)
 
+    photons = {}
+    emitter = {}
     photo_electrons = {}
     for o in telescope_data:
-        check_type(o, iact.PhotoElectrons)
-        photo_electrons[o.telescope_id] = o.parse()
-
-    return photo_electrons
+        if isinstance(o, iact.PhotoElectrons):
+            photo_electrons[o.telescope_id] = o.parse()
+        elif isinstance(o, iact.Photons):
+            p, e = o.parse()
+            photons[o.telescope_id] = p
+            if e is not None:
+                emitter[o.telescope_id] = e
+
+    return photons, emitter, photo_electrons
 
 
 def parse_telescope_event(telescope_event):

diff --git a/tests/resources/lst_with_photons.simtel.zst b/tests/resources/lst_with_photons.simtel.zst
diff --git a/tests/simtel/test_simtelfile.py b/tests/simtel/test_simtelfile.py
@@ -68,6 +68,8 @@ def test_show_event_is_not_empty_and_has_some_members_for_sure():
                 'tracking_positions',
                 'photoelectron_sums',
                 'photoelectrons',
+                'photons',
+                'emitter',
                 'camera_monitorings',
                 'laser_calibrations',
             }
@@ -128,7 +130,11 @@ def test_iterate_mc_events():
     expected = 200
     with SimTelFile(prod4_path) as f:
         for counter, event in enumerate(f.iter_mc_events(), start=1):
-            pass
+            assert event.keys() == {
+                'event_id',
+                'mc_shower', 'mc_event',
+                'photoelectrons', 'photons', 'emitter',
+            }
     assert counter == expected
 
 
@@ -192,3 +198,18 @@ def test_new_prod4():
         for e in f:
             i += 1
         assert i == 10
+
+
+def test_photons():
+    from eventio.iact.objects import Photons
+
+    with SimTelFile('tests/resources/lst_with_photons.simtel.zst') as f:
+        e = next(iter(f))
+
+        assert len(e['photons']) == 1
+        photons = e['photons'][0]
+        assert photons.dtype == Photons.long_dtype
+
+        # no emitter info in file
+        print(e['emitter'])
+        assert len(e['emitter']) == 0