Repository: cta-observatory/pyeventio

Issue:
Support ObservationLevelParticles object
When I read a simtel file that was produced with `IACT STORE-PARTICLES YES` I see the following warning:
```
<snip>/lib/python3.12/site-packages/eventio/simtel/simtelfile.py:252: UnknownObjectWarning: object type
encountered, which is no handled at the moment: ObservationLevelParticles(n_particles=17740)
  warnings.warn(
```
This seems to be because, although `iact.Photons` knows how to parse a photon table that contains observation-level particles, the `SimTelFile` reader doesn't know when to expect this particular Photons block (it comes right after the `MCShower` block, and only once per `MCShower` regardless of how many `MCEvent`) or what to do with it.

I would like the observation-level particles to be exposed in the event data structure. (Eventually I'd like them to be available from the ctapipe data structure as well, but that's another story.)

There is a bit of a question as to where the particle table should go. Because it is associated with an MC shower, it might make sense to insert it into the `mc_shower` structure. However, the particle positions are in a core coordinate system, so the user must add the core position to offset the particle positions in each event to the array coordinate system. For this reason I wonder if it might be more convenient for the user if the particle table is placed in the `mc_event` structure and already has the core position added to the positions to properly offset them.

Here's an example simulation file that contains the particle table: [gamma_20deg_0deg_run12___phase2d3_desert-ws0-nsb1.00.simhess.gz](https://github.com/user-attachments/files/19064613/gamma_20deg_0deg_run12___phase2d3_desert-ws0-nsb1.00.simhess.gz) (220 kB) – I wasn't sure about which types of simulation configuration are okay to share publicly, so this simulation uses HESS configurations that I believe are public.

---

I'm happy to take on this issue and create a PR, if you want to support this use case.

Discussion:
@maxnoe: Hi,

Sure, in general we want to be able to read all information from simtel files.

Regarding the question where to put it, I think it should just be a new field on the event. mc_shower and mc_event directly relate to the corresponding objects, so I think adding it there would be wrong.

Just create a new field in the event for the `observation_level_particles`@maxnoe: Could you maybe also fix the typo in the warning while you are at it? That would be great, thanks @maxnoe: >  I wasn't sure about which types of simulation configuration are okay to share publicly

The CTAO ones are published on zenodo: https://zenodo.org/records/6218687

Solution patch:
diff --git a/src/eventio/simtel/simtelfile.py b/src/eventio/simtel/simtelfile.py
index 063fb91..ae4b87e 100644
--- a/src/eventio/simtel/simtelfile.py
+++ b/src/eventio/simtel/simtelfile.py
@@ -112,6 +112,7 @@ class SimTelFile(EventIOFile):
         self.current_photoelectrons = {}
         self.current_photons = {}
         self.current_emitter = {}
+        self.current_obslev_particles = None
         self.current_array_event = None
         self.current_calibration_event = None
         self.current_calibration_event_id = None
@@ -161,6 +162,9 @@ class SimTelFile(EventIOFile):
             self.current_mc_shower = o.parse()
             self.current_mc_shower_id = o.header.id
 
+        elif isinstance(o, iact.Photons):
+            self.current_obslev_particles = o.parse()
+
         elif isinstance(o, ArrayEvent):
             self.current_array_event = parse_array_event(
                 o,
@@ -249,8 +253,8 @@ class SimTelFile(EventIOFile):
             self.atmospheric_profiles.append(o.parse())
         else:
             warnings.warn(
-                'object type encountered, which is no handled'
-                ' at the moment: {}'.format(o),
+                'encountered object of type "{}" which is not handled'
+                ' at the moment'.format(o),
                 UnknownObjectWarning,
             )
 
@@ -270,6 +274,7 @@ class SimTelFile(EventIOFile):
                 'event_id': self.current_mc_event_id,
                 'mc_shower': self.current_mc_shower,
                 'mc_event': self.current_mc_event,
+                'observation_level_particles': self.current_obslev_particles,
             }
             # if next object is TelescopeData, it belongs to this event
             if isinstance(self.peek(), iact.TelescopeData):
@@ -313,6 +318,7 @@ class SimTelFile(EventIOFile):
                 'event_id': event_id,
                 'mc_shower': None,
                 'mc_event': None,
+                'observation_level_particles': None,
                 'telescope_events': self.current_array_event['telescope_events'],
                 'tracking_positions': self.current_array_event['tracking_positions'],
                 'trigger_information': self.current_array_event['trigger_information'],
@@ -325,6 +331,7 @@ class SimTelFile(EventIOFile):
             if self.current_mc_event_id == event_id:
                 event_data['mc_shower'] = self.current_mc_shower
                 event_data['mc_event'] =  self.current_mc_event
+                event_data['observation_level_particles'] = self.current_obslev_particles
 
             if self.current_telescope_data_event_id == event_id:
                 event_data['photons'] = self.current_photons


Test patch:
diff --git a/tests/iact/test_iact_objects.py b/tests/iact/test_iact_objects.py
index e5a431b..c79a9b3 100644
--- a/tests/iact/test_iact_objects.py
+++ b/tests/iact/test_iact_objects.py
@@ -152,6 +152,30 @@ def test_emitter_bunches():
         # second lightest particle should be a muon
         assert np.isclose(np.unique(emitter['mass'])[1], 0.105, rtol=0.1)
 
+def test_obslev_particle_bunches():
+    from eventio.iact import Photons
+
+    columns = ('x', 'y', 'cx', 'cy', 'time', 'momentum', 'weight', 'particle_id')
+    with eventio.EventIOFile(test_emitter_file) as f:
+        # observation level particles object should be 5th in the file
+        for i in range(6):
+            obj = next(f)
+
+        assert isinstance(obj, Photons)
+        particles = obj.parse()
+
+        assert particles.dtype.names == columns
+        assert len(particles) == 14
+        assert np.all(particles['weight'] == np.float32(1))
+
+        sort_particles = np.sort(particles, order = 'momentum')
+
+        # lowest momentum particle should be a photon
+        assert np.floor(sort_particles[0]['particle_id']/1000) == np.float32(1)
+
+        # highest momentum particle should be a muon
+        assert np.floor(sort_particles[-1]['particle_id']/1000) == np.float32(5)
+
 
 def test_atmospheric_profile():
     from eventio.iact.objects import AtmosphericProfile
diff --git a/tests/simtel/test_simtelfile.py b/tests/simtel/test_simtelfile.py
index a4a2d54..3689315 100644
--- a/tests/simtel/test_simtelfile.py
+++ b/tests/simtel/test_simtelfile.py
@@ -33,7 +33,7 @@ def test_at_least_one_event_found():
         assert one_found, path
 
 
-def test_show_we_get_a_dict_with_hower_and_event():
+def test_show_we_get_a_dict_with_shower_and_event():
     for path in test_paths:
         for event in SimTelFile(path):
             assert 'mc_shower' in event
@@ -66,6 +66,7 @@ def test_show_event_is_not_empty_and_has_some_members_for_sure():
                 'event_id',
                 'mc_shower',
                 'mc_event',
+                'observation_level_particles',
                 'telescope_events',
                 'trigger_information',
                 'tracking_positions',
@@ -145,6 +146,7 @@ def test_iterate_mc_events():
             assert event.keys() == {
                 'event_id',
                 'mc_shower', 'mc_event',
+                'observation_level_particles',
                 'photons', 'photoelectrons', 'emitter'
             }
 
@@ -243,6 +245,24 @@ def test_missing_photons():
         assert e['emitter'] == {}
 
 
+def test_particles():
+    from eventio.iact.objects import Photons
+
+    with SimTelFile('tests/resources/gamma_with_obslevel_particles.simtel.zst') as f:
+        e = next(iter(f))
+
+        assert len(e['observation_level_particles']) == 64
+        particle = e['observation_level_particles'][0]
+        assert particle.dtype == Photons.particle_dtype
+
+
+def test_missing_particles():
+    with SimTelFile('tests/resources/gamma_test.simtel.gz') as f:
+        e = next(iter(f))
+
+        assert e['observation_level_particles'] == None
+
+
 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)):