From 271600c026824df63331a3c1e59670bb5920fba3 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Thu, 19 Oct 2023 15:13:15 +0200 Subject: [PATCH] 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``. --- ctapipe/containers.py | 1 + ctapipe/io/simteleventsource.py | 42 ++++++++++++++-------- ctapipe/io/tests/test_simteleventsource.py | 22 ++++++++++++ docs/changes/2411.features.rst | 3 ++ 4 files changed, 54 insertions(+), 14 deletions(-) create mode 100644 docs/changes/2411.features.rst diff --git a/ctapipe/containers.py b/ctapipe/containers.py index cb4282c5d4c..df8c8621b2d 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 632277a566c..70e8e93f81a 100644 --- a/ctapipe/io/simteleventsource.py +++ b/ctapipe/io/simteleventsource.py @@ -49,7 +49,7 @@ 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 @@ def _make_dummy_location(self): 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 @@ def _generate_events(self): 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 @@ def _parse_simulation_header(self): 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 @@ def _parse_simulation_header(self): 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 @@ def _fill_scheduling_and_observation_blocks(self): """ 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 @@ def _fill_scheduling_and_observation_blocks(self): } 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/ctapipe/io/tests/test_simteleventsource.py b/ctapipe/io/tests/test_simteleventsource.py index b4d3c0054ff..116a8a91f8a 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 diff --git a/docs/changes/2411.features.rst b/docs/changes/2411.features.rst new file mode 100644 index 00000000000..d81a88b977f --- /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.