Repository: cta-observatory/ctapipe (python)

Issue:
StereoTrigger does only work when tels_with_trigger is a python list, real sources use numpy arrays
**Describe the bug**

The code assumed and the unit tests made it so that the `StereoTrigger` works only when `tels_with_trigger` is a python list,
however the event source implementations provide a `np.ndarray` with dtype `uint16`.


**To Reproduce**
Try processing a prod3 or prod6 divergent pointing file where there might be single LSTs with the stereo trigger set to simulate the LST hardware stereo trigger:

```
❯ ctapipe-process -i dataset://gamma_divergent_LaPalma_baseline_20Zd_180Az_prod3_test.simtel.gz \
   -o /tmp/divergent.dl1.h5 \
    --progress --overwrite \
   --SimTelEventSource.focal_length_choice=EQUIVALENT \
   --SoftwareTrigger.min_telescopes_of_type="[('type', '*', 0), ('type', 'LST*', 2)]"

/home/mnoethe/.local/conda/envs/cta-dev/lib/python3.9/site-packages/traitlets/traitlets.py:2930: FutureWarning: --SoftwareTrigger.min_telescopes_of_type=[('type', '*', 0), ('type', 'LST*', 2)] for containers is deprecated in traitlets 5.0. You can pass `--SoftwareTrigger.min_telescopes_of_type item` ... multiple times to add items to a list.
  warn(
2023-05-08 12:08:24,527 WARNING [ctapipe.core.telescope_component] (telescope_component.attach_subarray): TelescopeParameter type argument 'SST_1M_*' did not match any known telescope types
SimTelEventSource: 1ev [00:02,  2.17s/ev]
2023-05-08 12:08:26,718 ERROR [ctapipe.ctapipe-process] (tool.run): Caught unexpected exception: 'numpy.ndarray' object has no attribute 'remove'
Traceback (most recent call last):
  File "/home/mnoethe/Uni/CTA/ctapipe/ctapipe/core/tool.py", line 426, in run
    self.start()
  File "/home/mnoethe/Uni/CTA/ctapipe/ctapipe/tools/process.py", line 305, in start
    if not self.software_trigger(event):
  File "/home/mnoethe/Uni/CTA/ctapipe/ctapipe/instrument/trigger.py", line 99, in __call__
    event.trigger.tels_with_trigger.remove(tel_id)
AttributeError: 'numpy.ndarray' object has no attribute 'remove'
```


**Expected behavior**

Single LSTs are removed, events with less than 2 telescopes are completely discarded.

**Supporting information**



**Additional context**


This is needed for the prod6 divergent processing, since due to an error in the simtel configuration, the LST hardware stereo trigger was not actually already included in the simulation.

Pull request:
Fix SoftwareTrigger for actual simtel files
Fixes #2321 

Patch:
diff --git a/ctapipe/instrument/trigger.py b/ctapipe/instrument/trigger.py
index 31120160..f53d2cc1 100644
--- a/ctapipe/instrument/trigger.py
+++ b/ctapipe/instrument/trigger.py
@@ -1,3 +1,5 @@
+import numpy as np
+
 from ctapipe.containers import ArrayEventContainer
 from ctapipe.core import TelescopeComponent
 from ctapipe.core.traits import Integer, IntTelescopeParameter
@@ -76,6 +78,7 @@ class SoftwareTrigger(TelescopeComponent):
             Whether or not this event would have triggered the stereo trigger
         """
 
+        tels_removed = set()
         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]
@@ -87,6 +90,7 @@ class SoftwareTrigger(TelescopeComponent):
             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(
@@ -96,7 +100,7 @@ class SoftwareTrigger(TelescopeComponent):
                     )
 
                     # remove from tels_with_trigger
-                    event.trigger.tels_with_trigger.remove(tel_id)
+                    tels_removed.add(tel_id)
 
                     # remove any related data
                     for container in ("trigger", "r0", "r1", "dl0", "dl1", "dl2"):
@@ -104,6 +108,13 @@ class SoftwareTrigger(TelescopeComponent):
                         if tel_id in tel_map:
                             del tel_map[tel_id]
 
+        if len(tels_removed) > 0:
+            # convert to array with correct dtype to have setdiff1d work correctly
+            tels_removed = np.fromiter(tels_removed, np.uint16, len(tels_removed))
+            event.trigger.tels_with_trigger = np.setdiff1d(
+                event.trigger.tels_with_trigger, tels_removed, assume_unique=True
+            )
+
         if len(event.trigger.tels_with_trigger) < self.min_telescopes:
             event.trigger.tels_with_trigger = []
             for container in ("trigger", "r0", "r1", "dl0", "dl1", "dl2"):
diff --git a/docs/changes/2320.bugfix.rst b/docs/changes/2320.bugfix.rst
new file mode 100644
index 00000000..384fb167
--- /dev/null
+++ b/docs/changes/2320.bugfix.rst
@@ -0,0 +1,2 @@
+Fix StereoTrigger assuming the wrong data type for ``tels_with_trigger``, resulting in
+it not working for actual events read from an EventSource.


Test patch:
diff --git a/ctapipe/instrument/tests/test_trigger.py b/ctapipe/instrument/tests/test_trigger.py
index 3371cdfa..ed6dec91 100644
--- a/ctapipe/instrument/tests/test_trigger.py
+++ b/ctapipe/instrument/tests/test_trigger.py
@@ -1,7 +1,13 @@
+import numpy as np
+import pytest
+from numpy.testing import assert_equal
+
 from ctapipe.containers import ArrayEventContainer
+from ctapipe.io import EventSource
 
 
-def test_software_trigger(subarray_prod5_paranal):
+@pytest.mark.parametrize("data_type", (list, np.array))
+def test_software_trigger(subarray_prod5_paranal, data_type):
     from ctapipe.instrument.trigger import SoftwareTrigger
 
     subarray = subarray_prod5_paranal
@@ -16,37 +22,83 @@ def test_software_trigger(subarray_prod5_paranal):
 
     # only one telescope, no SWAT
     event = ArrayEventContainer()
-    event.trigger.tels_with_trigger = [5]
+    event.trigger.tels_with_trigger = data_type([5])
     assert trigger(event) == False
-    assert event.trigger.tels_with_trigger == []
+    assert_equal(event.trigger.tels_with_trigger, data_type([]))
 
     # 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]
+    event.trigger.tels_with_trigger = data_type([1, 6])
     assert trigger(event) == False
-    assert event.trigger.tels_with_trigger == []
+    assert_equal(event.trigger.tels_with_trigger, data_type([]))
 
     # two MSTs and 1 LST, -> remove single LST
     event = ArrayEventContainer()
-    event.trigger.tels_with_trigger = [1, 5, 6]
+    event.trigger.tels_with_trigger = data_type([1, 5, 6])
     assert trigger(event) == True
-    assert event.trigger.tels_with_trigger == [5, 6]
+    assert_equal(event.trigger.tels_with_trigger, data_type([5, 6]))
 
     # two MSTs, nothing to change
     event = ArrayEventContainer()
-    event.trigger.tels_with_trigger = [5, 6]
+    event.trigger.tels_with_trigger = data_type([5, 6])
     assert trigger(event) == True
-    assert event.trigger.tels_with_trigger == [5, 6]
+    assert_equal(event.trigger.tels_with_trigger, data_type([5, 6]))
 
     # three LSTs, nothing to change
     event = ArrayEventContainer()
-    event.trigger.tels_with_trigger = [1, 2, 3]
+    event.trigger.tels_with_trigger = data_type([1, 2, 3])
     assert trigger(event) == True
-    assert event.trigger.tels_with_trigger == [1, 2, 3]
+    assert_equal(event.trigger.tels_with_trigger, data_type([1, 2, 3]))
 
     # thee LSTs, plus MSTs, nothing to change
     event = ArrayEventContainer()
-    event.trigger.tels_with_trigger = [1, 2, 3, 5, 6, 7]
+    event.trigger.tels_with_trigger = data_type([1, 2, 3, 5, 6, 7])
     assert trigger(event) == True
-    assert event.trigger.tels_with_trigger == [1, 2, 3, 5, 6, 7]
+    assert_equal(event.trigger.tels_with_trigger, data_type([1, 2, 3, 5, 6, 7]))
+
+
+@pytest.mark.parametrize("allowed_tels", (None, list(range(1, 20))))
+def test_software_trigger_simtel(allowed_tels):
+    from ctapipe.instrument.trigger import SoftwareTrigger
+
+    path = "dataset://gamma_divergent_LaPalma_baseline_20Zd_180Az_prod3_test.simtel.gz"
+
+    expected = [
+        [12, 16],
+        [],
+        [1, 2, 3, 4],
+        [1, 4],
+        [],
+        [1, 3],
+        [1, 2, 3, 4, 5, 6, 7, 12, 15, 16, 17, 18],
+        [13, 14],
+        [],
+        [2, 3, 7, 12],
+        [1, 2, 5, 17],
+        [],
+        [13, 19],
+        [],
+        [],
+        [1, 2, 4, 5, 11, 18],
+        [17, 18],
+        [7, 12],
+        [],
+    ]
+
+    with EventSource(
+        path, focal_length_choice="EQUIVALENT", allowed_tels=allowed_tels
+    ) as source:
+
+        trigger = SoftwareTrigger(
+            subarray=source.subarray,
+            min_telescopes=2,
+            min_telescopes_of_type=[
+                ("type", "*", 0),
+                ("type", "LST*", 2),
+            ],
+        )
+
+        for e, expected_tels in zip(source, expected):
+            trigger(e)
+            assert_equal(e.trigger.tels_with_trigger, expected_tels)