Repository: cta-observatory/ctapipe
Pull request #2411 (2023-10-19 13:16:58)

Add option to override obs_id in SimTelEventSource
Simulations were performed that re-use the same CORSIKA run numbers, e.g. different pointing positions for the same particles in prod5 and also different particles for the same pointing.

To still enable unique obs_ids in merged files for these productions, an "override_obs_id" option is added to the ``SimTelEventSource`` that can be used to assign new, globally unique ids.

The original run number is added to the ``SimulationConfigContainer``.

We have discussed before to decouple simulation obs_ids from CORSIKA / simtel run numbers.

This is a first step in enabling this. This was triggered also because the LST MC production reuses run numbers in a way that completely prevented merging files for model training.

maxnoe: I need this now for really the case where I need to assign custom obs_ids. The LST simulations restart run numbers for each pointing and I need to merge ~10 pointing directions for model training.

I wouldn't just add the `4`, since `ctapipe` is used on not just `SDMC-DPPS` produced simtel files, in fact, no such files really exist as of now.

Patch:
diff --git a/ctapipe/containers.py b/ctapipe/containers.py
index cb4282c5..df8c8621 100644
--- a/ctapipe/containers.py
+++ b/ctapipe/containers.py
@@ -772,6 +772,7 @@ class SimulationConfigContainer(Container):
     Configuration parameters of the simulation
     """
 
+    run_number = Field(np.int32(-1), description="Original sim_telarray run number")
     corsika_version = Field(nan, description="CORSIKA version * 1000")
     simtel_version = Field(nan, description="sim_telarray version * 1000")
     energy_range_min = Field(
diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py
index 632277a5..70e8e93f 100644
--- a/ctapipe/io/simteleventsource.py
+++ b/ctapipe/io/simteleventsource.py
@@ -49,7 +49,7 @@ from ..containers import (
 from ..coordinates import CameraFrame, shower_impact_distance
 from ..core import Map
 from ..core.provenance import Provenance
-from ..core.traits import Bool, ComponentName, Float, Undefined, UseEnum
+from ..core.traits import Bool, ComponentName, Float, Integer, Undefined, UseEnum
 from ..instrument import (
     CameraDescription,
     CameraGeometry,
@@ -483,6 +483,12 @@ class SimTelEventSource(EventSource):
         ),
     ).tag(config=True)
 
+    override_obs_id = Integer(
+        default_value=None,
+        allow_none=True,
+        help="Use the given obs_id instead of the run number from sim_telarray",
+    ).tag(config=True)
+
     def __init__(self, input_url=Undefined, config=None, parent=None, **kwargs):
         """
         EventSource for simtelarray files using the pyeventio library.
@@ -701,11 +707,10 @@ class SimTelEventSource(EventSource):
             SimulationConfigContainer), but with the lat/lon set to (0,0), i.e.
             on "Null Island"
         """
-        obs_id = self.file_.header["run"]
         return EarthLocation(
             lon=0 * u.deg,
             lat=0 * u.deg,
-            height=self._simulation_config[obs_id].prod_site_alt,
+            height=self._simulation_config[self.obs_id].prod_site_alt,
         )
 
     @staticmethod
@@ -743,7 +748,7 @@ class SimTelEventSource(EventSource):
                 pseudo_event_id -= 1
                 event_id = pseudo_event_id
 
-            obs_id = self.file_.header["run"]
+            obs_id = self.obs_id
 
             trigger = self._fill_trigger_info(array_event)
             if trigger.event_type == EventType.SUBARRAY:
@@ -1007,12 +1012,14 @@ class SimTelEventSource(EventSource):
         point in time, this dictionary will always have
         length 1.
         """
-        assert len(self.obs_ids) == 1
-        obs_id = self.obs_ids[0]
+        obs_id = self.obs_id
         # With only one run, we can take the first entry:
         mc_run_head = self.file_.mc_run_headers[-1]
 
         simulation_config = SimulationConfigContainer(
+            run_number=self.file_.header["run"],
+            detector_prog_start=mc_run_head["detector_prog_start"],
+            detector_prog_id=mc_run_head["detector_prog_id"],
             corsika_version=mc_run_head["shower_prog_vers"],
             simtel_version=mc_run_head["detector_prog_vers"],
             energy_range_min=mc_run_head["E_range"][0] * u.TeV,
@@ -1024,8 +1031,6 @@ class SimTelEventSource(EventSource):
             spectral_index=mc_run_head["spectral_index"],
             shower_prog_start=mc_run_head["shower_prog_start"],
             shower_prog_id=mc_run_head["shower_prog_id"],
-            detector_prog_start=mc_run_head["detector_prog_start"],
-            detector_prog_id=mc_run_head["detector_prog_id"],
             n_showers=mc_run_head["n_showers"],
             shower_reuse=mc_run_head["n_use"],
             max_alt=_clip_altitude_if_close(mc_run_head["alt_range"][1]) * u.rad,
@@ -1056,11 +1061,20 @@ class SimTelEventSource(EventSource):
         """
 
         az, alt = self.file_.header["direction"]
-        obs_id = self.file_.header["run"]
+
+        # this event source always contains only a single OB, so we can
+        # also assign a single obs_id
+        if self.override_obs_id is not None:
+            self.obs_id = self.override_obs_id
+        else:
+            self.obs_id = self.file_.header["run"]
+
+        # simulations at the moment do not have SBs, use OB id
+        self.sb_id = self.obs_id
 
         sb_dict = {
-            obs_id: SchedulingBlockContainer(
-                sb_id=np.uint64(obs_id),  # simulations have no SBs, use OB id
+            self.sb_id: SchedulingBlockContainer(
+                sb_id=np.uint64(self.sb_id),
                 sb_type=SchedulingBlockType.OBSERVATION,
                 producer_id="simulation",
                 observing_mode=ObservingMode.UNKNOWN,
@@ -1069,9 +1083,9 @@ class SimTelEventSource(EventSource):
         }
 
         ob_dict = {
-            obs_id: ObservationBlockContainer(
-                obs_id=np.uint64(obs_id),
-                sb_id=np.uint64(obs_id),  # see comment above
+            self.obs_id: ObservationBlockContainer(
+                obs_id=np.uint64(self.obs_id),
+                sb_id=np.uint64(self.sb_id),
                 producer_id="simulation",
                 state=ObservationBlockState.COMPLETED_SUCCEDED,
                 subarray_pointing_lat=alt * u.rad,
diff --git a/docs/changes/2411.features.rst b/docs/changes/2411.features.rst
new file mode 100644
index 00000000..d81a88b9
--- /dev/null
+++ b/docs/changes/2411.features.rst
@@ -0,0 +1,3 @@
+Add option ``override_obs_id`` to ``SimTelEventSource`` which allows
+assigning new, unique ``obs_ids`` in case productions re-use CORSIKA run
+numbers.


Test patch:
diff --git a/ctapipe/io/tests/test_simteleventsource.py b/ctapipe/io/tests/test_simteleventsource.py
index b4d3c005..116a8a91 100644
--- a/ctapipe/io/tests/test_simteleventsource.py
+++ b/ctapipe/io/tests/test_simteleventsource.py
@@ -604,3 +604,25 @@ def test_starting_grammage():
     with SimTelEventSource(path, focal_length_choice="EQUIVALENT") as source:
         e = next(iter(source))
         assert e.simulation.shower.starting_grammage == 580 * u.g / u.cm**2
+
+
+@pytest.mark.parametrize("override_obs_id,expected_obs_id", [(None, 1), (5, 5)])
+def test_override_obs_id(override_obs_id, expected_obs_id, prod5_gamma_simtel_path):
+    """Test for the override_obs_id option"""
+    original_run_number = 1
+
+    with SimTelEventSource(
+        prod5_gamma_simtel_path, override_obs_id=override_obs_id
+    ) as s:
+        assert s.obs_id == expected_obs_id
+        assert s.obs_ids == [expected_obs_id]
+
+        assert s.simulation_config.keys() == {expected_obs_id}
+        assert s.observation_blocks.keys() == {expected_obs_id}
+        assert s.scheduling_blocks.keys() == {expected_obs_id}
+
+        # this should alway be the original run number
+        assert s.simulation_config[s.obs_id].run_number == original_run_number
+
+        for e in s:
+            assert e.index.obs_id == expected_obs_id