Merge pull request #2735 from cta-observatory/simtel_true_image

Improve handling of missing true images

Author: kosack

docs/changes/2735.bugfix.rst

@@ -0,0 +1,4 @@
+Improve heuristic in ``SimTelEventSource`` when true Cherenkov
+images are expected to be present in the input file to fix missing
+true images in output files when the first event of a run
+had missing true images.

src/ctapipe/io/simteleventsource.py

@@ -121,6 +121,57 @@ def _clip_altitude_if_close(altitude):
     return altitude
 
 
+def _split_history(history):
+    tel_history = {0: []}
+    tel = 0
+
+    for timestamp, line in history:
+        line = line.decode("utf-8").strip()
+
+        if line.startswith("# Telescope-specific configuration follows"):
+            tel += 1
+            tel_history[tel] = []
+
+        tel_history[tel].append((timestamp, line))
+
+    global_history = tel_history.pop(0)
+    return global_history, tel_history
+
+
+def _has_true_image(history):
+    global_history, tel_history = _split_history(history)
+
+    # merge_simtel does not include photoelectrons, so always missing
+    if "merge_simtel" in global_history[0][1]:
+        return False
+
+    # first, check for global default
+    save_photons = 0
+    global_store_pe = -1
+    for _, line in global_history:
+        if line.lower().startswith("save_photons"):
+            save_photons = int(line.split()[1])
+
+        if line.lower().startswith("store_photoelectrons"):
+            global_store_pe = int(line.split()[1])
+
+    if (save_photons & 0b10) > 0:
+        return True
+
+    if global_store_pe > 0:
+        return True
+
+    # if global default is not set, we still might have true images
+    # via the telescope configuration
+    store_photoelectrons = {}
+    for tel, config in tel_history.items():
+        for _, line in config:
+            if "store_photoelectrons" in line.lower():
+                store_photoelectrons[tel] = int(line.split()[1])
+
+    return any(v > 0 for v in store_photoelectrons.values())
+
+
 @enum.unique
 class MirrorClass(enum.Enum):
     """Enum for the sim_telarray MIRROR_CLASS values"""
@@ -575,7 +626,13 @@ def __init__(self, input_url=Undefined, config=None, parent=None, **kwargs):
             self.file_, kind=self.atmosphere_profile_choice
         )
 
-        self._has_true_image = None
+        try:
+            self._has_true_image = _has_true_image(self.file_.history)
+        except Exception:
+            self.log.exception(
+                "Error trying to determine whether file has true_image, assuming yes"
+            )
+            self._has_true_image = True
         self.log.debug(f"Using gain selector {self.gain_selector}")
 
     def __exit__(self, exc_type, exc_val, exc_tb):
@@ -895,9 +952,6 @@ def _generate_events(self):
                     self.telescope_indices_original[tel_id]
                 ]
 
-                if self._has_true_image is None:
-                    self._has_true_image = true_image is not None
-
                 if self._has_true_image and true_image is None:
                     if true_image_sum > MISSING_IMAGE_BRIGHTNESS_LIMIT:
                         self.log.warning(

src/ctapipe/io/tests/test_simteleventsource.py

@@ -688,3 +688,28 @@ def test_prod6_issues():
 
     assert strange_trigger_events.keys() == events_checked_trigger
     assert missing_true_images.keys() == events_checked_image
+
+
+@pytest.mark.parametrize(
+    ("path", "expected", "kwargs"),
+    [
+        ("dataset://gamma_prod5.simtel.zst", True, {}),
+        ("dataset://test_pe_first_missing.simtel.zst", True, {}),
+        (gamma_test_large_path, False, {"focal_length_choice": "EQUIVALENT"}),
+    ],
+)
+def test_has_true_image(path, expected, kwargs):
+    with SimTelEventSource(path, **kwargs) as source:
+        assert source._has_true_image == expected
+
+
+def test_has_true_image_first_missing():
+    """Check that we have true images even if the first event is missing it"""
+    input_url = "dataset://test_pe_first_missing.simtel.zst"
+
+    with SimTelEventSource(input_url) as source:
+        n_found = 0
+        for event in source:
+            assert event.simulation.tel[1].true_image is not None
+            n_found += 1
+        assert n_found == 2