Implement stereo trigger for subarray triggers

ctapipe/instrument/__init__.py

@@ -4,6 +4,7 @@
 from .optics import FocalLengthKind, OpticsDescription, ReflectorShape, SizeType
 from .subarray import SubarrayDescription, UnknownTelescopeID
 from .telescope import TelescopeDescription
+from .trigger import SoftwareTrigger
 
 __all__ = [
     "CameraDescription",
@@ -19,4 +20,5 @@
     "FocalLengthKind",
     "ReflectorShape",
     "SizeType",
+    "SoftwareTrigger",
 ]

ctapipe/instrument/tests/test_trigger.py

@@ -0,0 +1,52 @@
+from ctapipe.containers import ArrayEventContainer
+
+
+def test_software_trigger(subarray_prod5_paranal):
+    from ctapipe.instrument.trigger import SoftwareTrigger
+
+    subarray = subarray_prod5_paranal
+    trigger = SoftwareTrigger(
+        subarray=subarray,
+        min_telescopes=2,
+        min_telescopes_of_type=[
+            ("type", "*", 0),
+            ("type", "LST*", 2),
+        ],
+    )
+
+    # only one telescope, no SWAT
+    event = ArrayEventContainer()
+    event.trigger.tels_with_trigger = [5]
+    assert trigger(event) == False
+    assert event.trigger.tels_with_trigger == []
+
+    # 1 LST + 1 MST, 1 LST would not have triggered LST hardware trigger
+    # and after LST is removed, we only have 1 telescope, so no SWAT either
+    event = ArrayEventContainer()
+    event.trigger.tels_with_trigger = [1, 6]
+    assert trigger(event) == False
+    assert event.trigger.tels_with_trigger == []
+
+    # two MSTs and 1 LST, -> remove single LST
+    event = ArrayEventContainer()
+    event.trigger.tels_with_trigger = [1, 5, 6]
+    assert trigger(event) == True
+    assert event.trigger.tels_with_trigger == [5, 6]
+
+    # two MSTs, nothing to change
+    event = ArrayEventContainer()
+    event.trigger.tels_with_trigger = [5, 6]
+    assert trigger(event) == True
+    assert event.trigger.tels_with_trigger == [5, 6]
+
+    # three LSTs, nothing to change
+    event = ArrayEventContainer()
+    event.trigger.tels_with_trigger = [1, 2, 3]
+    assert trigger(event) == True
+    assert event.trigger.tels_with_trigger == [1, 2, 3]
+
+    # thee LSTs, plus MSTs, nothing to change
+    event = ArrayEventContainer()
+    event.trigger.tels_with_trigger = [1, 2, 3, 5, 6, 7]
+    assert trigger(event) == True
+    assert event.trigger.tels_with_trigger == [1, 2, 3, 5, 6, 7]

ctapipe/instrument/trigger.py

@@ -0,0 +1,113 @@
+from ctapipe.containers import ArrayEventContainer
+from ctapipe.core import TelescopeComponent
+from ctapipe.core.traits import Integer, IntTelescopeParameter
+
+
+class SoftwareTrigger(TelescopeComponent):
+    """
+    A stereo trigger that can remove telescope events from subarray events.
+
+    This class is needed to correctly handle super-arrays as simulated for
+    CTA and still handle the LST hardware stereo trigger and the normal stereo
+    trigger correctly.
+
+    When selecting subarrays from simulations that contain many more telescopes,
+    as is done in all major CTA productions to date, the stereo trigger is not
+    correctly simulated as in that after selecting a realistic subarray, events
+    are still in the data stream where only one telescope of the selected subarray
+    triggered, which would in reality not trigger the stereo trigger.
+
+    An additional complexity is the LST hardware stereo trigger, that forces that
+    an array event has always to contain no or at least two LST telescope events.
+
+    This means that after selectig a subarray, we need to:
+    - Remove LST telescope events from the subarray if only one LST triggered
+    - Ignore events with only 1 telescope after this has been applied
+
+    With the default settings, this class is a no-op. To get the correct behavior
+    for CTA simulations, use the following configuration:
+
+    ..
+        SoftwareTrigger:
+            min_telescopes: 2
+            min_telescopes_of_type:
+                - ["type", "*", 0]
+                - ["type", "LST*", 2]
+    """
+
+    min_telescopes = Integer(
+        default_value=1,
+        help=(
+            "Minimum number of telescopes required globally."
+            " Events with fewer telescopes will be filtered out completely."
+        ),
+    ).tag(config=True)
+
+    min_telescopes_of_type = IntTelescopeParameter(
+        default_value=0,
+        help=(
+            "Minimum number of telescopes required for a specific type."
+            " In events with fewer telescopes of that type"
+            " , those telescopes will be removed from the array event."
+            " This might result in the event not fullfilling ``min_telescopes`` anymore"
+            " and thus being filtered completely."
+        ),
+    ).tag(config=True)
+
+    def __init__(self, subarray, *args, **kwargs):
+        super().__init__(subarray, *args, **kwargs)
+
+        self._ids_by_type = {
+            str(type): set(self.subarray.get_tel_ids_for_type(type))
+            for type in self.subarray.telescope_types
+        }
+
+    def __call__(self, event: ArrayEventContainer) -> bool:
+        """
+        Remove telescope events that have not the required number of telescopes of
+        a given type from the subarray event and decide if the event would
+        have triggered the stereo trigger.
+
+        Data is cleared from events that did not trigger.
+
+        Returns
+        -------
+        triggered : bool
+            Whether or not this event would have triggered the stereo trigger
+        """
+
+        for tel_type in self.subarray.telescope_types:
+            tel_type_str = str(tel_type)
+            min_tels = self.min_telescopes_of_type.tel[tel_type_str]
+
+            # no need to check telescopes for which we have no min requirement
+            if min_tels == 0:
+                continue
+
+            tels_with_trigger = set(event.trigger.tels_with_trigger)
+            tel_ids = self._ids_by_type[tel_type_str]
+            tels_in_event = tels_with_trigger.intersection(tel_ids)
+            if len(tels_in_event) < min_tels:
+                for tel_id in tels_in_event:
+                    self.log.debug(
+                        "Removing tel_id %d of type %s from event due to type requirement",
+                        tel_id,
+                        tel_type_str,
+                    )
+
+                    # remove from tels_with_trigger
+                    event.trigger.tels_with_trigger.remove(tel_id)
+
+                    # remove any related data
+                    for container in ("trigger", "r0", "r1", "dl0", "dl1", "dl2"):
+                        tel_map = getattr(event, container).tel
+                        if tel_id in tel_map:
+                            del tel_map[tel_id]
+
+        if len(event.trigger.tels_with_trigger) < self.min_telescopes:
+            event.trigger.tels_with_trigger = []
+            for container in ("trigger", "r0", "r1", "dl0", "dl1", "dl2"):
+                tel_map = getattr(event, container).tel
+                tel_map.clear()
+            return False
+        return True

ctapipe/tools/process.py

@@ -11,6 +11,7 @@
 from ..core.traits import Bool, classes_with_traits, flag
 from ..image import ImageCleaner, ImageModifier, ImageProcessor
 from ..image.extractor import ImageExtractor
+from ..instrument import SoftwareTrigger
 from ..io import (
     DataLevel,
     DataWriter,
@@ -144,6 +145,7 @@ class ProcessorTool(Tool):
             ShowerProcessor,
             metadata.Instrument,
             metadata.Contact,
+            SoftwareTrigger,
         ]
         + classes_with_traits(EventSource)
         + classes_with_traits(ImageCleaner)
@@ -169,17 +171,11 @@ def setup(self):
             )
             sys.exit(1)
 
-        self.calibrate = CameraCalibrator(
-            parent=self, subarray=self.event_source.subarray
-        )
-        self.process_images = ImageProcessor(
-            subarray=self.event_source.subarray, parent=self
-        )
-
-        self.process_shower = ShowerProcessor(
-            subarray=self.event_source.subarray, parent=self
-        )
-
+        subarray = self.event_source.subarray
+        self.software_trigger = SoftwareTrigger(parent=self, subarray=subarray)
+        self.calibrate = CameraCalibrator(parent=self, subarray=subarray)
+        self.process_images = ImageProcessor(subarray=subarray, parent=self)
+        self.process_shower = ShowerProcessor(subarray=subarray, parent=self)
         self.write = DataWriter(event_source=self.event_source, parent=self)
 
         # add ml reco classes if model paths were supplied via cli and not already configured
@@ -196,7 +192,7 @@ def setup(self):
                 reconstructor = Reconstructor.from_name(
                     name,
                     parent=self.process_shower,
-                    subarray=self.event_source.subarray,
+                    subarray=subarray,
                 )
                 self.process_shower.reconstructors.append(reconstructor)
                 self.process_shower.reconstructor_types.append(name)
@@ -298,6 +294,12 @@ def start(self):
             if not self.event_type_filter(event):
                 continue
 
+            if not self.software_trigger(event):
+                self.log.debug(
+                    "Skipping event %i due to software trigger", event.index.event_id
+                )
+                continue
+
             if self.should_calibrate:
                 self.calibrate(event)