Repository: cta-observatory/ctapipe
Pull request #1929 (2022-06-10 07:21:44)

Support reading files created with `extract_calibevent`
Test file to be added.


Patch:
diff --git a/ctapipe/coordinates/__init__.py b/ctapipe/coordinates/__init__.py
index d32680ed..5a210b55 100644
--- a/ctapipe/coordinates/__init__.py
+++ b/ctapipe/coordinates/__init__.py
@@ -10,12 +10,7 @@ from astropy.coordinates import (
 import warnings
 from .telescope_frame import TelescopeFrame
 from .nominal_frame import NominalFrame
-from .ground_frames import (
-    GroundFrame,
-    TiltedGroundFrame,
-    project_to_ground,
-    EastingNorthingFrame,
-)
+from .ground_frames import GroundFrame, TiltedGroundFrame, project_to_ground
 from .camera_frame import CameraFrame, EngineeringCameraFrame
 
 
@@ -26,7 +21,6 @@ __all__ = [
     "NominalFrame",
     "GroundFrame",
     "TiltedGroundFrame",
-    "EastingNorthingFrame",
     "MissingFrameAttributeWarning",
     "project_to_ground",
 ]
diff --git a/ctapipe/coordinates/ground_frames.py b/ctapipe/coordinates/ground_frames.py
index 8b19bb4d..7285607b 100644
--- a/ctapipe/coordinates/ground_frames.py
+++ b/ctapipe/coordinates/ground_frames.py
@@ -14,25 +14,21 @@ TODO:
 - Tests Tests Tests!
 """
 import astropy.units as u
-from astropy.units.quantity import Quantity
 import numpy as np
 from astropy.coordinates import (
-    AltAz,
     BaseCoordinateFrame,
     CartesianRepresentation,
-    CoordinateAttribute,
     FunctionTransform,
-    RepresentationMapping,
-    frame_transform_graph,
-    AffineTransform,
+    CoordinateAttribute,
+    AltAz,
 )
+from astropy.coordinates import frame_transform_graph
 from numpy import cos, sin
 
 __all__ = [
     "GroundFrame",
     "TiltedGroundFrame",
     "project_to_ground",
-    "EastingNorthingFrame",
 ]
 
 
@@ -51,25 +47,6 @@ class GroundFrame(BaseCoordinateFrame):
     default_representation = CartesianRepresentation
 
 
-class EastingNorthingFrame(BaseCoordinateFrame):
-    """GroundFrame but in standard Easting/Northing coordinates instead of
-    SimTel/Corsika conventions
-
-    Frame attributes: None
-
-    """
-
-    default_representation = CartesianRepresentation
-
-    frame_specific_representation_info = {
-        CartesianRepresentation: [
-            RepresentationMapping("x", "easting"),
-            RepresentationMapping("y", "northing"),
-            RepresentationMapping("z", "height"),
-        ]
-    }
-
-
 class TiltedGroundFrame(BaseCoordinateFrame):
     """Tilted ground coordinate frame.  The tilted ground coordinate frame
     is a cartesian system describing the 2 dimensional projected
@@ -223,34 +200,3 @@ def project_to_ground(tilt_system):
     y_projected = y_initial - trans[2][1] * z_initial / trans[2][2]
 
     return GroundFrame(x=x_projected * unit, y=y_projected * unit, z=0 * unit)
-
-
-@frame_transform_graph.transform(FunctionTransform, GroundFrame, GroundFrame)
-def ground_to_ground(ground_coords, ground_frame):
-    """Null transform for going from ground to ground, since there are no
-    attributes of the GroundSystem"""
-    return ground_coords
-
-
-# Matrices for transforming between GroundFrame and EastingNorthingFrame
-NO_OFFSET = CartesianRepresentation(Quantity([0, 0, 0], u.m))
-GROUND_TO_EASTNORTH = np.asarray([[0, -1, 0], [1, 0, 0], [0, 0, 1]])
-
-
-@frame_transform_graph.transform(AffineTransform, GroundFrame, EastingNorthingFrame)
-def ground_to_easting_northing(ground_coords, eastnorth_frame):
-    """
-    convert GroundFrame points into eastings/northings for plotting purposes
-
-    """
-
-    return GROUND_TO_EASTNORTH, NO_OFFSET
-
-
-@frame_transform_graph.transform(AffineTransform, EastingNorthingFrame, GroundFrame)
-def easting_northing_to_ground(eastnorth_coords, ground_frame):
-    """
-    convert  eastings/northings back to GroundFrame
-
-    """
-    return GROUND_TO_EASTNORTH.T, NO_OFFSET
diff --git a/ctapipe/instrument/__init__.py b/ctapipe/instrument/__init__.py
index 561e992e..79442154 100644
--- a/ctapipe/instrument/__init__.py
+++ b/ctapipe/instrument/__init__.py
@@ -2,7 +2,7 @@ from .camera import CameraDescription, CameraGeometry, CameraReadout, PixelShape
 from .atmosphere import get_atmosphere_profile_functions
 from .telescope import TelescopeDescription
 from .optics import OpticsDescription
-from .subarray import SubarrayDescription, UnknownTelescopeID
+from .subarray import SubarrayDescription
 from .guess import guess_telescope
 
 
@@ -11,10 +11,9 @@ __all__ = [
     "CameraGeometry",
     "CameraReadout",
     "get_atmosphere_profile_functions",
-    "guess_telescope",
+    "TelescopeDescription",
     "OpticsDescription",
-    "PixelShape",
     "SubarrayDescription",
-    "TelescopeDescription",
-    "UnknownTelescopeID",
+    "guess_telescope",
+    "PixelShape",
 ]
diff --git a/ctapipe/instrument/subarray.py b/ctapipe/instrument/subarray.py
index 24ac847b..557040b4 100644
--- a/ctapipe/instrument/subarray.py
+++ b/ctapipe/instrument/subarray.py
@@ -1,30 +1,28 @@
 """
 Description of Arrays or Subarrays of telescopes
 """
-import warnings
-from contextlib import ExitStack
-from copy import copy
-from itertools import groupby
 from typing import Dict, List, Union
+from contextlib import ExitStack
+import warnings
 
 import numpy as np
-import tables
 from astropy import units as u
 from astropy.coordinates import SkyCoord
 from astropy.table import QTable, Table
 from astropy.utils import lazyproperty
+import tables
+from copy import copy
+from itertools import groupby
 
-from .. import __version__ as CTAPIPE_VERSION
-from ..coordinates import CameraFrame, GroundFrame
-from .camera import CameraDescription, CameraGeometry, CameraReadout
-from .optics import OpticsDescription
-from .telescope import TelescopeDescription
+import ctapipe
 
-__all__ = ["SubarrayDescription"]
+from ..coordinates import GroundFrame, CameraFrame
+from .telescope import TelescopeDescription
+from .camera import CameraDescription, CameraReadout, CameraGeometry
+from .optics import OpticsDescription
 
 
-class UnknownTelescopeID(KeyError):
-    """Raised when an unknown telescope id is encountered"""
+__all__ = ["SubarrayDescription"]
 
 
 def _group_consecutives(sequence):
@@ -33,7 +31,6 @@ def _group_consecutives(sequence):
 
     from https://codereview.stackexchange.com/questions/214820/codewars-range-extraction
     """
-    sequence = sorted(sequence)
     for _, g in groupby(enumerate(sequence), lambda i_x: i_x[0] - i_x[1]):
         r = [x for _, x in g]
         if len(r) > 2:
@@ -93,7 +90,7 @@ class SubarrayDescription:
 
     @property
     def tel(self):
-        """for backward compatibility"""
+        """ for backward compatibility"""
         return self.tels
 
     @property
@@ -135,7 +132,7 @@ class SubarrayDescription:
 
     @lazyproperty
     def tel_coords(self):
-        """returns telescope positions as astropy.coordinates.SkyCoord"""
+        """ returns telescope positions as astropy.coordinates.SkyCoord"""
 
         pos_x = np.array([p[0].to("m").value for p in self.positions.values()]) * u.m
         pos_y = np.array([p[1].to("m").value for p in self.positions.values()]) * u.m
@@ -145,7 +142,7 @@ class SubarrayDescription:
 
     @lazyproperty
     def tel_ids(self):
-        """telescope IDs as an array"""
+        """ telescope IDs as an array"""
         return np.array(list(self.tel.keys()))
 
     @lazyproperty
@@ -241,7 +238,7 @@ class SubarrayDescription:
         meta = {
             "ORIGIN": "ctapipe.instrument.SubarrayDescription",
             "SUBARRAY": self.name,
-            "SOFT_VER": CTAPIPE_VERSION,
+            "SOFT_VER": ctapipe.__version__,
             "TAB_TYPE": kind,
         }
 
@@ -280,8 +277,8 @@ class SubarrayDescription:
             unique_types = self.telescope_types
 
             mirror_area = u.Quantity(
-                [t.optics.mirror_area.to_value(u.m**2) for t in unique_types],
-                u.m**2,
+                [t.optics.mirror_area.to_value(u.m ** 2) for t in unique_types],
+                u.m ** 2,
             )
             focal_length = u.Quantity(
                 [t.optics.equivalent_focal_length.to_value(u.m) for t in unique_types],
@@ -319,15 +316,11 @@ class SubarrayDescription:
         SubarrayDescription
         """
 
-        unknown_tel_ids = set(tel_ids).difference(self.tel.keys())
-        if len(unknown_tel_ids) > 0:
-            known = _range_extraction(self.tel.keys())
-            raise UnknownTelescopeID(f"{unknown_tel_ids}, known telescopes: {known}")
-
         tel_positions = {tid: self.positions[tid] for tid in tel_ids}
         tel_descriptions = {tid: self.tel[tid] for tid in tel_ids}
 
         if not name:
+            tel_ids = sorted(tel_ids)
             name = self.name + "_" + _range_extraction(tel_ids)
 
         newsub = SubarrayDescription(
@@ -339,29 +332,47 @@ class SubarrayDescription:
         """
         Draw a quick matplotlib plot of the array
         """
-        from ctapipe.coordinates.ground_frames import EastingNorthingFrame
-        from ctapipe.visualization import ArrayDisplay
         from matplotlib import pyplot as plt
+        from astropy.visualization import quantity_support
+
+        types = set(self.tels.values())
+        tab = self.to_table()
 
         plt.figure(figsize=(8, 8))
-        ad = ArrayDisplay(subarray=self, frame=EastingNorthingFrame(), tel_scale=0.75)
-        ad.add_labels()
-        plt.title(self.name)
-        plt.tight_layout()
+
+        with quantity_support():
+            for tel_type in types:
+                tels = tab[tab["tel_description"] == str(tel_type)]["tel_id"]
+                sub = self.select_subarray(tels, name=tel_type)
+                tel_coords = sub.tel_coords
+                radius = np.array(
+                    [
+                        np.sqrt(tel.optics.mirror_area / np.pi).value
+                        for tel in sub.tels.values()
+                    ]
+                )
+
+                plt.scatter(
+                    tel_coords.x, tel_coords.y, s=radius * 8, alpha=0.5, label=tel_type
+                )
+
+            plt.legend(loc="best")
+            plt.title(self.name)
+            plt.tight_layout()
 
     @lazyproperty
     def telescope_types(self) -> List[TelescopeDescription]:
-        """list of telescope types in the array"""
+        """ list of telescope types in the array"""
         return list({tel for tel in self.tel.values()})
 
     @lazyproperty
     def camera_types(self) -> List[CameraDescription]:
-        """list of camera types in the array"""
+        """ list of camera types in the array """
         return list({tel.camera for tel in self.tel.values()})
 
     @lazyproperty
     def optics_types(self) -> List[OpticsDescription]:
-        """list of optics types in the array"""
+        """ list of optics types in the array """
         return list({tel.optics for tel in self.tel.values()})
 
     def get_tel_ids_for_type(self, tel_type):
diff --git a/ctapipe/io/datawriter.py b/ctapipe/io/datawriter.py
index ed8d811e..712e733a 100644
--- a/ctapipe/io/datawriter.py
+++ b/ctapipe/io/datawriter.py
@@ -268,7 +268,7 @@ class DataWriter(Component):
             table_name="dl1/event/subarray/trigger",
             containers=[event.index, event.trigger],
         )
-        if self._is_simulation:
+        if event.simulation is not None and event.simulation.shower is not None:
             self._writer.write(
                 table_name="simulation/event/subarray/shower",
                 containers=[event.index, event.simulation.shower],
diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py
index cdab4f59..8eb52425 100644
--- a/ctapipe/io/simteleventsource.py
+++ b/ctapipe/io/simteleventsource.py
@@ -398,6 +398,8 @@ class SimTelEventSource(EventSource):
         data.meta["origin"] = "hessio"
         data.meta["input_url"] = self.input_url
         data.meta["max_events"] = self.max_events
+        # for events without event_id, we use negative event_ids
+        pseudo_event_id = 0
 
         self._fill_array_pointing(data)
 
@@ -413,7 +415,11 @@ class SimTelEventSource(EventSource):
             data.simulation.tel.clear()
             data.trigger.tel.clear()
 
-            event_id = array_event.get("event_id", -1)
+            event_id = array_event.get("event_id", 0)
+            if event_id == 0:
+                pseudo_event_id -= 1
+                event_id = pseudo_event_id
+
             obs_id = self.file_.header["run"]
             data.count = counter
             data.index.obs_id = obs_id
@@ -426,9 +432,13 @@ class SimTelEventSource(EventSource):
             telescope_events = array_event["telescope_events"]
             tracking_positions = array_event["tracking_positions"]
 
-            true_image_sums = array_event.get("photoelectron_sums", {}).get(
-                "n_pe", np.full(self.n_telescopes_original, np.nan)
-            )
+            photoelectron_sums = array_event.get("photoelectron_sums")
+            if photoelectron_sums is not None:
+                true_image_sums = photoelectron_sums.get(
+                    "n_pe", np.full(self.n_telescopes_original, np.nan)
+                )
+            else:
+                true_image_sums = np.full(self.n_telescopes_original, np.nan)
 
             for tel_id, telescope_event in telescope_events.items():
                 adc_samples = telescope_event.get("adc_samples")
@@ -442,12 +452,13 @@ class SimTelEventSource(EventSource):
                     .get("photoelectrons", None)
                 )
 
-                data.simulation.tel[tel_id] = SimulatedCameraContainer(
-                    true_image_sum=true_image_sums[
-                        self.telescope_indices_original[tel_id]
-                    ],
-                    true_image=true_image,
-                )
+                if data.simulation is not None:
+                    data.simulation.tel[tel_id] = SimulatedCameraContainer(
+                        true_image_sum=true_image_sums[
+                            self.telescope_indices_original[tel_id]
+                        ],
+                        true_image=true_image,
+                    )
 
                 data.pointing.tel[tel_id] = self._fill_event_pointing(
                     tracking_positions[tel_id]
@@ -628,6 +639,9 @@ class SimTelEventSource(EventSource):
     def _fill_simulated_event_information(data, array_event):
         mc_event = array_event["mc_event"]
         mc_shower = array_event["mc_shower"]
+        if mc_shower is None:
+            data.simulation.shower = None
+            return
 
         data.simulation.shower = SimulatedShowerContainer(
             energy=u.Quantity(mc_shower["energy"], u.TeV),
diff --git a/ctapipe/visualization/mpl_array.py b/ctapipe/visualization/mpl_array.py
index 1c3b3e2c..5d895ec8 100644
--- a/ctapipe/visualization/mpl_array.py
+++ b/ctapipe/visualization/mpl_array.py
@@ -62,12 +62,10 @@ class ArrayDisplay:
 
         self.frame = frame
         self.subarray = subarray
-        self.axes = axes or plt.gca()
 
         # get the telescope positions. If a new frame is set, this will
         # transform to the new frame.
-        self.tel_coords = subarray.tel_coords.transform_to(frame).cartesian
-        self.unit = self.tel_coords.x.unit
+        self.tel_coords = subarray.tel_coords.transform_to(frame)
 
         # set up colors per telescope type
         tel_types = [str(tel) for tel in subarray.tels.values()]
@@ -78,10 +76,6 @@ class ArrayDisplay:
                 for tel in subarray.tel.values()
             ]
 
-            self.radii = radius
-        else:
-            self.radii = np.ones(len(tel_types)) * radius
-
         if title is None:
             title = subarray.name
 
@@ -122,22 +116,15 @@ class ArrayDisplay:
             )
         plt.legend(handles=legend_elements)
 
-        self.add_radial_grid()
-
-        # create the plot
         self.tel_colors = tel_color
         self.autoupdate = autoupdate
         self.telescopes = PatchCollection(patches, match_original=True)
         self.telescopes.set_linewidth(2.0)
 
+        self.axes = axes or plt.gca()
         self.axes.add_collection(self.telescopes)
         self.axes.set_aspect(1.0)
         self.axes.set_title(title)
-        xunit = self.tel_coords.x.unit.to_string("latex")
-        yunit = self.tel_coords.y.unit.to_string("latex")
-        xname, yname, _ = frame.get_representation_component_names().keys()
-        self.axes.set_xlabel(f"{xname} [{xunit}] $\\rightarrow$")
-        self.axes.set_ylabel(f"{yname} [{yunit}] $\\rightarrow$")
         self._labels = []
         self._quiver = None
         self.axes.autoscale_view()
@@ -149,50 +136,12 @@ class ArrayDisplay:
 
     @values.setter
     def values(self, values):
-        """set the telescope colors to display"""
+        """ set the telescope colors to display  """
         self.telescopes.set_array(np.ma.masked_invalid(values))
         self._update()
 
-    def add_radial_grid(self, spacing=100 * u.m):
-        """add some dotted rings for distance estimation. The number of rings
-        is estimated automatically from the spacing and the array footprint.
-
-        Parameters
-        ----------
-        spacing: Quantity
-            spacing between rings
-
-        """
-
-        n_circles = np.round(
-            (np.sqrt(self.subarray.footprint / np.pi) / spacing).to_value(""),
-            0,
-        )
-        circle_radii = np.arange(1, n_circles + 2, 1) * spacing.to_value(self.unit)
-        circle_patches = PatchCollection(
-            [
-                Circle(
-                    xy=(0, 0),
-                    radius=r,
-                    fill=False,
-                    fc="none",
-                    linestyle="dotted",
-                    color="gray",
-                    alpha=0.1,
-                    lw=1,
-                )
-                for r in circle_radii
-            ],
-            color="#eeeeee",
-            ls="dotted",
-            fc="none",
-            lw=3,
-        )
-
-        self.axes.add_collection(circle_patches)
-
     def set_vector_uv(self, uu, vv, c=None, **kwargs):
-        """sets the vector field U,V and color for all telescopes
+        """ sets the vector field U,V and color for all telescopes
 
         Parameters
         ----------
@@ -340,17 +289,9 @@ class ArrayDisplay:
     def add_labels(self):
         px = self.tel_coords.x.to_value("m")
         py = self.tel_coords.y.to_value("m")
-        for tel, x, y, r in zip(self.subarray.tels, px, py, self.radii):
+        for tel, x, y in zip(self.subarray.tels, px, py):
             name = str(tel)
-            lab = self.axes.text(
-                x,
-                y - r * 1.8,
-                name,
-                fontsize=8,
-                clip_on=True,
-                horizontalalignment="center",
-                verticalalignment="top",
-            )
+            lab = self.axes.text(x, y, name, fontsize=8, clip_on=True)
             self._labels.append(lab)
 
     def remove_labels(self):
@@ -359,7 +300,7 @@ class ArrayDisplay:
         self._labels = []
 
     def _update(self):
-        """signal a redraw if necessary"""
+        """ signal a redraw if necessary """
         if self.autoupdate:
             plt.draw()
 
diff --git a/environment.yml b/environment.yml
index 0db4cbba..37d72749 100644
--- a/environment.yml
+++ b/environment.yml
@@ -4,7 +4,7 @@ channels:
   - conda-forge
   - default
 dependencies:
-  - python=3.8
+  - python=3.9
   - pip
   - astropy=5
   - black
@@ -45,5 +45,5 @@ dependencies:
   - xz
   - zlib
   - zstandard
-  - eventio>=1.5.0
+  - eventio>=1.9.1
   - jinja2=3.0 # for sphinx 3.5, remove when updating to 4.x
diff --git a/setup.py b/setup.py
index 8dff48a4..73a7c2e0 100755
--- a/setup.py
+++ b/setup.py
@@ -45,7 +45,7 @@ setup(
     install_requires=[
         "astropy~=5.0",
         "bokeh~=2.0",
-        "eventio>=1.5.0,<2.0.0a0",
+        "eventio>=1.9.1,<2.0.0a0",
         "h5py",
         "iminuit>=2",
         "joblib",


Test patch:
diff --git a/ctapipe/coordinates/tests/test_coordinates.py b/ctapipe/coordinates/tests/test_coordinates.py
index ef75a58d..23206d67 100644
--- a/ctapipe/coordinates/tests/test_coordinates.py
+++ b/ctapipe/coordinates/tests/test_coordinates.py
@@ -257,7 +257,6 @@ def test_camera_focal_length_array():
 
 
 def test_ground_frame_roundtrip():
-    """test transform from sky to ground roundtrip"""
     from ctapipe.coordinates import GroundFrame, TiltedGroundFrame
 
     normal = SkyCoord(alt=70 * u.deg, az=0 * u.deg, frame=AltAz())
@@ -269,22 +268,3 @@ def test_ground_frame_roundtrip():
     assert u.isclose(coord.x, back.x, atol=1e-12 * u.m)
     assert u.isclose(coord.y, back.y, atol=1e-12 * u.m)
     assert u.isclose(coord.z, back.z, atol=1e-12 * u.m)
-
-
-def test_ground_to_eastnorth_roundtrip():
-    """Check Ground to EastingNorthing and the round-trip"""
-    from ctapipe.coordinates import GroundFrame, EastingNorthingFrame
-
-    ground = SkyCoord(
-        x=[1, 2, 3] * u.m, y=[-2, 5, 2] * u.m, z=[1, -1, 2] * u.m, frame=GroundFrame()
-    )
-    eastnorth = ground.transform_to(EastingNorthingFrame())
-    ground2 = eastnorth.transform_to(GroundFrame())
-
-    assert u.isclose(eastnorth.easting, [2, -5, -2] * u.m).all()
-    assert u.isclose(eastnorth.northing, [1, 2, 3] * u.m).all()
-    assert u.isclose(eastnorth.height, [1, -1, 2] * u.m).all()
-
-    assert u.isclose(ground.x, ground2.x).all()
-    assert u.isclose(ground.y, ground2.y).all()
-    assert u.isclose(ground.z, ground2.z).all()
\ No newline at end of file
diff --git a/ctapipe/instrument/tests/test_subarray.py b/ctapipe/instrument/tests/test_subarray.py
index c594730c..333957cf 100644
--- a/ctapipe/instrument/tests/test_subarray.py
+++ b/ctapipe/instrument/tests/test_subarray.py
@@ -16,7 +16,7 @@ from ctapipe.instrument import (
 
 
 def example_subarray(n_tels=10):
-    """generate a simple subarray for testing purposes"""
+    """ generate a simple subarray for testing purposes """
     rng = np.random.default_rng(0)
 
     pos = {}
@@ -32,7 +32,7 @@ def example_subarray(n_tels=10):
 
 
 def test_subarray_description():
-    """Test SubarrayDescription functionality"""
+    """ Test SubarrayDescription functionality """
     n_tels = 10
     sub = example_subarray(n_tels)
     sub.peek()
@@ -67,14 +67,14 @@ def test_subarray_description():
 
 
 def test_to_table(example_subarray):
-    """Check that we can generate astropy Tables from the SubarrayDescription"""
+    """ Check that we can generate astropy Tables from the SubarrayDescription """
     sub = example_subarray
     assert len(sub.to_table(kind="subarray")) == sub.num_tels
     assert len(sub.to_table(kind="optics")) == len(sub.optics_types)
 
 
 def test_tel_indexing(example_subarray):
-    """Check that we can convert between telescope_id and telescope_index"""
+    """ Check that we can convert between telescope_id and telescope_index """
     sub = example_subarray
 
     assert sub.tel_indices[1] == 0  # first tel_id is in slot 0
@@ -218,10 +218,3 @@ def test_get_tel_ids(example_subarray):
     # test invalid telescope type
     with pytest.raises(Exception):
         tel_ids = subarray.get_tel_ids(["It's a-me, Mario!"])
-
-
-def test_unknown_telescopes(example_subarray):
-    from ctapipe.instrument import UnknownTelescopeID
-
-    with pytest.raises(UnknownTelescopeID):
-        example_subarray.select_subarray([300, 201])
diff --git a/ctapipe/io/tests/test_simteleventsource.py b/ctapipe/io/tests/test_simteleventsource.py
index 77eba9e8..935a8582 100644
--- a/ctapipe/io/tests/test_simteleventsource.py
+++ b/ctapipe/io/tests/test_simteleventsource.py
@@ -20,7 +20,9 @@ from ctapipe.io import DataLevel
 gamma_test_large_path = get_dataset_path("gamma_test_large.simtel.gz")
 gamma_test_path = get_dataset_path("gamma_test.simtel.gz")
 calib_events_path = get_dataset_path("lst_prod3_calibration_and_mcphotons.simtel.zst")
-prod5b_path = get_dataset_path("gamma_20deg_0deg_run2___cta-prod5-paranal_desert-2147m-Paranal-dark_cone10-100evts.simtel.zst")
+prod5b_path = get_dataset_path(
+    "gamma_20deg_0deg_run2___cta-prod5-paranal_desert-2147m-Paranal-dark_cone10-100evts.simtel.zst"
+)
 
 
 def test_positional_input():
@@ -32,8 +34,9 @@ def test_simtel_event_source_on_gamma_test_one_event():
     assert SimTelEventSource.is_compatible(gamma_test_large_path)
 
     with SimTelEventSource(
-        input_url=gamma_test_large_path, back_seekable=True,
-        focal_length_choice='nominal',
+        input_url=gamma_test_large_path,
+        back_seekable=True,
+        focal_length_choice="nominal",
     ) as reader:
         assert not reader.is_stream
 
@@ -66,7 +69,7 @@ def test_that_event_is_not_modified_after_loop():
 def test_additional_meta_data_from_simulation_config():
     with SimTelEventSource(
         input_url=gamma_test_large_path,
-        focal_length_choice='nominal',
+        focal_length_choice="nominal",
     ) as reader:
         pass
 
@@ -108,7 +111,7 @@ def test_additional_meta_data_from_simulation_config():
 def test_properties():
     source = SimTelEventSource(
         input_url=gamma_test_large_path,
-        focal_length_choice='nominal',
+        focal_length_choice="nominal",
     )
 
     assert source.is_simulation
@@ -123,7 +126,7 @@ def test_gamma_file_prod2():
     with pytest.warns(UnknownPixelShapeWarning):
         with SimTelEventSource(
             input_url=dataset,
-            focal_length_choice='nominal',
+            focal_length_choice="nominal",
         ) as reader:
             assert reader.is_compatible(dataset)
             assert reader.is_stream  # using gzip subprocess makes it a stream
@@ -144,7 +147,7 @@ def test_max_events():
     with SimTelEventSource(
         input_url=gamma_test_large_path,
         max_events=max_events,
-        focal_length_choice='nominal',
+        focal_length_choice="nominal",
     ) as reader:
         count = 0
         for _ in reader:
@@ -156,7 +159,7 @@ def test_pointing():
     with SimTelEventSource(
         input_url=gamma_test_large_path,
         max_events=3,
-        focal_length_choice='nominal',
+        focal_length_choice="nominal",
     ) as reader:
         for e in reader:
             assert np.isclose(e.pointing.array_altitude.to_value(u.deg), 70)
@@ -174,8 +177,10 @@ def test_allowed_telescopes():
     # test that the allowed_tels mask works:
     allowed_tels = {3, 4}
     with SimTelEventSource(
-        input_url=gamma_test_large_path, allowed_tels=allowed_tels, max_events=5,
-        focal_length_choice='nominal',
+        input_url=gamma_test_large_path,
+        allowed_tels=allowed_tels,
+        max_events=5,
+        focal_length_choice="nominal",
     ) as reader:
         assert not allowed_tels.symmetric_difference(reader.subarray.tel_ids)
         for event in reader:
@@ -202,20 +207,25 @@ def test_calibration_events():
         EventType.SUBARRAY,
         EventType.SUBARRAY,
     ]
+
+    expected_ids = [-1, -2, -3, -4, -5, -6, -7, -8, 100, 200]
     with SimTelEventSource(
-        input_url=calib_events_path, skip_calibration_events=False,
-        focal_length_choice='nominal',
+        input_url=calib_events_path,
+        skip_calibration_events=False,
+        focal_length_choice="nominal",
     ) as reader:
-
-        for event, expected_type in zip_longest(reader, expected_types):
+        for event, expected_type, expected_id in zip_longest(
+            reader, expected_types, expected_ids
+        ):
             assert event.trigger.event_type is expected_type
+            assert event.index.event_id == expected_id
 
 
 def test_trigger_times():
 
     source = SimTelEventSource(
         input_url=calib_events_path,
-        focal_length_choice='nominal',
+        focal_length_choice="nominal",
     )
     t0 = Time("2020-05-06T15:30:00")
     t1 = Time("2020-05-06T15:40:00")
@@ -230,7 +240,7 @@ def test_trigger_times():
 def test_true_image():
     with SimTelEventSource(
         input_url=calib_events_path,
-        focal_length_choice='nominal',
+        focal_length_choice="nominal",
     ) as reader:
 
         for event in reader:
@@ -242,7 +252,7 @@ def test_instrument():
     """Test if same telescope types share a single instance of CameraGeometry"""
     source = SimTelEventSource(
         input_url=gamma_test_large_path,
-        focal_length_choice='nominal',
+        focal_length_choice="nominal",
     )
     subarray = source.subarray
     assert subarray.tel[1].optics.num_mirrors == 1
@@ -309,22 +319,25 @@ def test_focal_length_choice():
         SimTelEventSource(gamma_test_large_path)
 
     with pytest.raises(RuntimeError):
-        SimTelEventSource(gamma_test_large_path, focal_length_choice='effective')
+        SimTelEventSource(gamma_test_large_path, focal_length_choice="effective")
 
-    s = SimTelEventSource(gamma_test_large_path, focal_length_choice='nominal')
+    s = SimTelEventSource(gamma_test_large_path, focal_length_choice="nominal")
     assert s.subarray.tel[1].optics.equivalent_focal_length == 28 * u.m
 
     # this file does
-    s = SimTelEventSource(prod5b_path, focal_length_choice='effective')
-    assert u.isclose(s.subarray.tel[1].optics.equivalent_focal_length, 29.3 * u.m, atol=0.05 * u.m)
+    s = SimTelEventSource(prod5b_path, focal_length_choice="effective")
+    assert u.isclose(
+        s.subarray.tel[1].optics.equivalent_focal_length, 29.3 * u.m, atol=0.05 * u.m
+    )
     # check guessing of the name is not affected by focal length choice
-    assert str(s.subarray.tel[1]) == 'LST_LST_LSTCam'
+    assert str(s.subarray.tel[1]) == "LST_LST_LSTCam"
 
-    s = SimTelEventSource(prod5b_path, focal_length_choice='nominal')
-    assert u.isclose(s.subarray.tel[1].optics.equivalent_focal_length, 28.0 * u.m, atol=0.05 * u.m)
+    s = SimTelEventSource(prod5b_path, focal_length_choice="nominal")
+    assert u.isclose(
+        s.subarray.tel[1].optics.equivalent_focal_length, 28.0 * u.m, atol=0.05 * u.m
+    )
     # check guessing of the name is not affected by focal length choice
-    assert str(s.subarray.tel[1]) == 'LST_LST_LSTCam'
-
+    assert str(s.subarray.tel[1]) == "LST_LST_LSTCam"
 
 
 def test_only_config():
@@ -379,14 +392,14 @@ def test_true_image_sum():
     # this file does not contain true pe info
     with SimTelEventSource(
         gamma_test_large_path,
-        focal_length_choice='nominal',
+        focal_length_choice="nominal",
     ) as s:
         e = next(iter(s))
         assert np.all(np.isnan(sim.true_image_sum) for sim in e.simulation.tel.values())
 
     with SimTelEventSource(
         calib_events_path,
-        focal_length_choice='nominal',
+        focal_length_choice="nominal",
     ) as s:
         e = next(iter(s))
 
@@ -400,9 +413,21 @@ def test_true_image_sum():
     # check it also works with allowed_tels, since the values
     # are stored in a flat array in simtel
     with SimTelEventSource(
-        calib_events_path, allowed_tels={2, 3},
-        focal_length_choice='nominal',
+        calib_events_path,
+        allowed_tels={2, 3},
+        focal_length_choice="nominal",
     ) as s:
         e = next(iter(s))
         assert e.simulation.tel[2].true_image_sum == true_image_sums[2]
         assert e.simulation.tel[3].true_image_sum == true_image_sums[3]
+
+
+def test_extracted_calibevents():
+    with SimTelEventSource("dataset://extracted_pedestals.simtel.zst") as s:
+        i = 0
+        for e in s:
+            i = e.count
+            # these events are simulated but do not have shower information
+            assert e.simulation is not None
+            assert e.simulation.shower is None
+        assert i == 4