Merge pull request #295 from zhipzhang/split_always_mode

Add support for split-always mode event handling in IACTFile

Author: maxnoe

src/eventio/iact/__init__.py

@@ -18,6 +18,8 @@
     RunEnd,
     Longitudinal,
     InputCard,
+    TelescopeArrayHead,
+    TelescopeArrayEnd,
     AtmosphericProfile,
 )
 
@@ -43,6 +45,8 @@
     'RunEnd',
     'Longitudinal',
     'InputCard',
+    'TelescopeArrayHead',
+    'TelescopeArrayEnd',
     'AtmosphericProfile',
 ]
 
@@ -80,6 +84,18 @@ class IACTFile(EventIOFile):
       EventEnd
 
     RunEnd
+
+    For Split-Always mode, the structure is slightly different:
+    For each Event:
+      EventHeader
+      ArrayOffsets
+      Longitudinal (optional)
+      For each reuse:
+        TelescopeArrayHead
+        For each Telescopes:
+          Photons
+        TelescopeArrayEnd
+      EventEnd
     '''
 
     def __init__(self, path, zcat=True):
@@ -140,7 +156,7 @@ def __iter__(self):
 
             obj = next(self)
 
-            while not isinstance(obj, TelescopeData):
+            while not isinstance(obj, (TelescopeData, TelescopeArrayHead)):
                 if isinstance(obj, Longitudinal):
                     longitudinal = obj.parse()
 
@@ -152,22 +168,34 @@ def __iter__(self):
 
                 obj = next(self)
 
-            for reuse in range(n_reuses):
-
-                check_type(obj, TelescopeData)
-                telescope_data_obj = obj
+            split_always = isinstance(obj, TelescopeArrayHead)
+            if split_always:
+                obj = next(self)
 
+            for reuse in range(n_reuses):
                 photon_bunches = {}
                 emitter_bunches = {}
                 n_photons = {}
                 n_bunches = {}
-                for data in telescope_data_obj:
-                    if isinstance(data, Photons):
-                        photons, emitter = data.parse()
-                        photon_bunches[data.telescope] = photons
-                        emitter_bunches[data.telescope] = emitter
-                        n_photons[data.telescope] = data.n_photons
-                        n_bunches[data.telescope] = data.n_bunches
+                if split_always:
+                    while not isinstance(obj, TelescopeArrayEnd):
+                        check_type(obj, Photons)
+                        photons, emitter = obj.parse()
+                        photon_bunches[obj.telescope] = photons
+                        emitter_bunches[obj.telescope] = emitter
+                        n_photons[obj.telescope] = obj.n_photons
+                        n_bunches[obj.telescope] = obj.n_bunches
+                        obj = next(self)
+                else:
+                    check_type(obj, TelescopeData)
+                    telescope_data_obj = obj
+                    for data in telescope_data_obj:
+                        if isinstance(data, Photons):
+                            photons, emitter = data.parse()
+                            photon_bunches[data.telescope] = photons
+                            emitter_bunches[data.telescope] = emitter
+                            n_photons[data.telescope] = data.n_photons
+                            n_bunches[data.telescope] = data.n_bunches
 
                 if len(array_offsets.dtype) == 3:
                     weight = array_offsets[reuse]['weight']

src/eventio/iact/objects.py

@@ -32,6 +32,8 @@
     'RunEnd',
     'Longitudinal',
     'InputCard',
+    'MarkersBegin',
+    'MakersEnd',
 ]
 
 
@@ -412,6 +414,13 @@ def parse(self):
             input_card += b'\n'
         return input_card
 
+class TelescopeArrayHead(EventIOObject):
+    """This Object is just a marker to indcate the beginning of separate telescope data."""
+    eventio_type = 1213
+
+class TelescopeArrayEnd(EventIOObject):
+    """This Object is just a marker to indcate the end of separate telescope data."""
+    eventio_type = 1214
 
 class AtmosphericProfile(EventIOObject):
     ''' This Object contains the CORSIKA/atmext atmospheric profile '''

tests/iact/test_iact_file.py

@@ -7,6 +7,7 @@
 testfile_reuse = 'tests/resources/3_gammas_reuse_5.dat'
 testfile_two_telescopes = 'tests/resources/two_telescopes.dat'
 testfile_zstd = 'tests/resources/run102_gamma_za20deg_azm0deg-paranal-sst.corsika.zst'
+testfile_split_always = 'tests/resources/5_gamma_2_telescopes_split_always.zst'
 
 
 def test_file_open():
@@ -135,3 +136,12 @@ def test_event_with_reuse():
         for i, e in enumerate(f):
             assert e.event_number == i // 5 + 1
             assert e.reuse == (i % 5) + 1
+
+
+def test_event_split_always():
+    with eventio.IACTFile(testfile_split_always) as f:
+        event = next(iter(f))
+        assert event.header['zenith'] == approx(20.0/ 180 * 3.14159)
+        assert event.header['azimuth'] == approx(0.0)
+        assert event.header['total_energy'] == approx(50)
+        assert len(event.photon_bunches) == 2