Repository: cta-observatory/ctapipe
Pull request #1943 (2022-06-22 15:13:27)

Add QualityQuery.get_table_mask, compile python code
In preparation of implemting the quality query for tables, I tried to come up with a way these quality queries are pure expressions, that would be suitable for numexpr also in the "apply to single element case".

I do this here by evaluating the expressions with the kwargs of call as locals. The expressions are compiled to bytecode in `__init__` and then evaled, which should not make a difference in performance.

I haven't measured this yet though.


Patch:
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 00b1bec5..bc377a3b 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -54,7 +54,7 @@ jobs:
           fetch-depth: 0
 
       - name: cache test data
-        uses: actions/cache@v2
+        uses: actions/cache@v3
         with:
           path: ~/.cache/ctapipe
           key: ctapipe-test-data
@@ -115,20 +115,20 @@ jobs:
 
           ctapipe-info --version
 
-      - uses: codecov/codecov-action@v1
+      - uses: codecov/codecov-action@v2
         if: matrix.codecov
 
   docs:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v3
         with:
           fetch-depth: 0
 
       - name: Set up Python
-        uses: actions/setup-python@v2
+        uses: actions/setup-python@v4
         with:
-          python-version: 3.8
+          python-version: "3.8"
 
       - name: Install doc dependencies
         run: |
diff --git a/ctapipe/containers.py b/ctapipe/containers.py
index a4a8c99f..f12f5c4a 100644
--- a/ctapipe/containers.py
+++ b/ctapipe/containers.py
@@ -92,7 +92,7 @@ class EventType(enum.Enum):
 
 
 class EventIndexContainer(Container):
-    """ index columns to include in event lists, common to all data levels"""
+    """index columns to include in event lists, common to all data levels"""
 
     container_prefix = ""  # don't want to prefix these
     obs_id = Field(0, "observation identifier")
@@ -251,7 +251,7 @@ class TimingParametersContainer(BaseTimingParametersContainer):
 
 
 class MorphologyContainer(Container):
-    """ Parameters related to pixels surviving image cleaning """
+    """Parameters related to pixels surviving image cleaning"""
 
     num_pixels = Field(-1, "Number of usable pixels")
     num_islands = Field(-1, "Number of distinct islands in the image")
@@ -287,7 +287,7 @@ class CoreParametersContainer(Container):
 
 
 class ImageParametersContainer(Container):
-    """ Collection of image parameters """
+    """Collection of image parameters"""
 
     container_prefix = "params"
     hillas = Field(
@@ -353,7 +353,7 @@ class DL1CameraContainer(Container):
 
 
 class DL1Container(Container):
-    """ DL1 Calibrated Camera Images and associated data"""
+    """DL1 Calibrated Camera Images and associated data"""
 
     tel = Field(Map(DL1CameraContainer), "map of tel_id to DL1CameraContainer")
 
@@ -466,6 +466,19 @@ class DL0Container(Container):
     tel = Field(Map(DL0CameraContainer), "map of tel_id to DL0CameraContainer")
 
 
+class TelescopeImpactParameterContainer(Container):
+    """
+    Impact Parameter computed from reconstructed shower geometry
+    """
+
+    container_prefix = "impact"
+
+    distance = Field(
+        nan * u.m, "distance of the telescope to the shower axis", unit=u.m
+    )
+    distance_uncert = Field(nan * u.m, "uncertainty in impact_parameter", unit=u.m)
+
+
 class SimulatedShowerContainer(Container):
     container_prefix = "true"
     energy = Field(nan * u.TeV, "Simulated Energy", unit=u.TeV)
@@ -475,7 +488,7 @@ class SimulatedShowerContainer(Container):
     core_y = Field(nan * u.m, "Simulated core position (y)", unit=u.m)
     h_first_int = Field(nan * u.m, "Height of first interaction", unit=u.m)
     x_max = Field(
-        nan * u.g / (u.cm ** 2), "Simulated Xmax value", unit=u.g / (u.cm ** 2)
+        nan * u.g / (u.cm**2), "Simulated Xmax value", unit=u.g / (u.cm**2)
     )
     shower_primary_id = Field(
         -1,
@@ -508,6 +521,8 @@ class SimulatedCameraContainer(Container):
         None, "Parameters derived from the true_image", type=ImageParametersContainer
     )
 
+    impact = Field(TelescopeImpactParameterContainer(), "true impact parameter")
+
 
 class SimulatedEventContainer(Container):
     shower = Field(SimulatedShowerContainer(), "True event information")
@@ -606,37 +621,33 @@ class ReconstructedGeometryContainer(Container):
         nan * u.m, "reconstructed x coordinate of the core position", unit=u.m
     )
     core_y = Field(
-        nan * u.m,
-        "reconstructed y coordinate of the core position",
-        unit=u.m
+        nan * u.m, "reconstructed y coordinate of the core position", unit=u.m
     )
     core_uncert_x = Field(
         nan * u.m,
         "reconstructed core position uncertainty along ground frame X axis",
-        unit=u.m
+        unit=u.m,
     )
     core_uncert_y = Field(
         nan * u.m,
         "reconstructed core position uncertainty along ground frame Y axis",
-        unit=u.m
+        unit=u.m,
     )
     core_tilted_x = Field(
         nan * u.m, "reconstructed x coordinate of the core position", unit=u.m
     )
     core_tilted_y = Field(
-        nan * u.m,
-        "reconstructed y coordinate of the core position",
-        unit=u.m
+        nan * u.m, "reconstructed y coordinate of the core position", unit=u.m
     )
     core_tilted_uncert_x = Field(
         nan * u.m,
         "reconstructed core position uncertainty along tilted frame X axis",
-        unit=u.m
+        unit=u.m,
     )
     core_tilted_uncert_y = Field(
         nan * u.m,
         "reconstructed core position uncertainty along tilted frame Y axis",
-        unit=u.m
+        unit=u.m,
     )
     h_max = Field(nan * u.m, "reconstructed height of the shower maximum", unit=u.m)
     h_max_uncert = Field(nan * u.m, "uncertainty of h_max", unit=u.m)
@@ -700,7 +711,7 @@ class ParticleClassificationContainer(Container):
 
 
 class ReconstructedContainer(Container):
-    """ Reconstructed shower info from multiple algorithms """
+    """Reconstructed shower info from multiple algorithms"""
 
     # Note: there is a reason why the hiererchy is
     # `event.dl2.stereo.geometry[algorithm]` and not
@@ -724,6 +735,15 @@ class ReconstructedContainer(Container):
     )
 
 
+class TelescopeReconstructedContainer(ReconstructedContainer):
+    """Telescope-wise reconstructed quantities"""
+
+    impact = Field(
+        Map(TelescopeImpactParameterContainer),
+        "map of algorithm to impact parameter info",
+    )
+
+
 class DL2Container(Container):
     """Reconstructed Shower information for a given reconstruction algorithm,
     including optionally both per-telescope mono reconstruction and per-shower
@@ -731,7 +751,7 @@ class DL2Container(Container):
     """
 
     tel = Field(
-        Map(ReconstructedContainer),
+        Map(TelescopeReconstructedContainer),
         "map of tel_id to single-telescope reconstruction (DL2a)",
     )
     stereo = Field(ReconstructedContainer(), "Stereo Shower reconstruction results")
@@ -1005,7 +1025,7 @@ class SimulatedShowerDistribution(Container):
 
 
 class ArrayEventContainer(Container):
-    """ Top-level container for all event information """
+    """Top-level container for all event information"""
 
     index = Field(EventIndexContainer(), "event indexing information")
     r0 = Field(R0Container(), "Raw Data")
diff --git a/ctapipe/coordinates/__init__.py b/ctapipe/coordinates/__init__.py
index d32680ed..8b4f127f 100644
--- a/ctapipe/coordinates/__init__.py
+++ b/ctapipe/coordinates/__init__.py
@@ -1,11 +1,15 @@
 """
 Coordinates.
 """
+
+import numpy as np
+
 from astropy.coordinates import (
     AltAz,
     FunctionTransformWithFiniteDifference,
     CIRS,
     frame_transform_graph,
+    spherical_to_cartesian,
 )
 import warnings
 from .telescope_frame import TelescopeFrame
@@ -29,6 +33,7 @@ __all__ = [
     "EastingNorthingFrame",
     "MissingFrameAttributeWarning",
     "project_to_ground",
+    "altaz_to_righthanded_cartesian",
 ]
 
 
@@ -65,3 +70,19 @@ def altaz_to_altaz(from_coo, to_frame):
         return to_frame.realize_frame(from_coo.spherical)
 
     return from_coo.transform_to(CIRS(obstime=obstime)).transform_to(to_frame)
+
+
+def altaz_to_righthanded_cartesian(alt, az):
+    """Turns an alt/az coordinate into a 3d direction in a right-handed coordinate
+    system.  This is because azimuths are in a left-handed system.
+
+    See e.g: https://github.com/astropy/astropy/issues/5300
+
+    Parameters
+    ----------
+    alt: u.Quantity
+        altitude
+    az: u.Quantity
+        azimuth
+    """
+    return np.array(spherical_to_cartesian(r=1, lat=alt, lon=-az))
diff --git a/ctapipe/core/provenance.py b/ctapipe/core/provenance.py
index b0cc65ee..ed916ecc 100644
--- a/ctapipe/core/provenance.py
+++ b/ctapipe/core/provenance.py
@@ -127,7 +127,7 @@ class Provenance(metaclass=Singleton):
         Parameters
         ----------
         config: dict
-            configuration paramters
+            configuration parameters
         """
         self.current_activity.register_config(config)
 
diff --git a/ctapipe/core/qualityquery.py b/ctapipe/core/qualityquery.py
index 8bfd3781..3fd23b31 100644
--- a/ctapipe/core/qualityquery.py
+++ b/ctapipe/core/qualityquery.py
@@ -4,8 +4,6 @@ Data Quality selection
 
 __all__ = ["QualityQuery", "QualityCriteriaError"]
 
-from collections.abc import Callable
-
 import astropy.units as u  # for use in selection functions
 import numpy as np  # for use in selection functions
 
@@ -16,11 +14,12 @@ from .traits import List
 # in selection functions (passed to eval())
 ALLOWED_GLOBALS = {"u": u, "np": np}  # astropy units  # numpy
 
+for func in ("sin", "cos", "tan", "arctan2", "log", "log10", "exp", "sqrt"):
+    ALLOWED_GLOBALS[func] = getattr(np, func)
 
-class QualityCriteriaError(TypeError):
-    """ Signal a problem with a user-defined selection criteria function"""
 
-    pass
+class QualityCriteriaError(TypeError):
+    """Signal a problem with a user-defined selection criteria function"""
 
 
 class QualityQuery(Component):
@@ -34,8 +33,8 @@ class QualityQuery(Component):
 
     quality_criteria = List(
         help=(
-            "list of tuples of ('<description', 'function string') to accept "
-            "(select) a given data value.  E.g. `[('mycut', 'lambda x: x > 3'),]. "
+            "list of tuples of ('<description', 'expression string') to accept "
+            "(select) a given data value.  E.g. `[('mycut', 'x > 3'),]. "
             "You may use `numpy` as `np` and `astropy.units` as `u`, but no other"
             " modules."
         )
@@ -45,40 +44,28 @@ class QualityQuery(Component):
         super().__init__(config=config, parent=parent, **kwargs)
 
         # add a selection to count all entries and make it the first one
-        self.quality_criteria.insert(0, ("TOTAL", "lambda x: True"))
-        self.criteria_names = []
-        self.selection_function_strings = []
-        self._selectors = []
-
-        for name, func_str in self.quality_criteria:
-            try:  # generate real functions from the selection function strings
-                self.criteria_names.append(name)
-                self.selection_function_strings.append(func_str)
-
-                func = eval(func_str, ALLOWED_GLOBALS)
-                if not isinstance(func, Callable):
-                    raise QualityCriteriaError(
-                        f"Selection criterion '{name}' cannot be evaluated because "
-                        f" '{func_str}' is not a callable function"
-                    )
-                self._selectors.append(func)
-
-            except NameError as err:
-                # catch functions that cannot be defined. Note that this cannot check
-                # that the function can run, that only happens the first time it's
-                # called.
+        self.criteria_names = [n for (n, _) in self.quality_criteria]
+        self.expressions = [e for (_, e) in self.quality_criteria]
+        self._compiled_expressions = []
+        for name, e in self.quality_criteria:
+            if "lambda" in e:
+                raise ValueError(
+                    "As of ctapipe 0.16, do not give lambda expressions"
+                    " to QualityQuery. Directly give the expression."
+                    " E.g. instead of `lambda p: p.hillas.width.value > 0`"
+                    " use `parameters.hillas.width.value > 0`"
+                )
+            try:
+                self._compiled_expressions.append(compile(e, __name__, mode="eval"))
+            except Exception:
                 raise QualityCriteriaError(
-                    f"Couldn't evaluate selection function '{name}' -> '{func_str}' "
-                    f"because: {err}"
+                    f"Error compiling expression '{e}' for check {name}"
                 )
 
-        # arrays for recording overall statistics
-        self._counts = np.zeros(len(self._selectors), dtype=np.int64)
-        self._cumulative_counts = np.zeros(len(self._selectors), dtype=np.int64)
-
-    def __len__(self):
-        """ return number of events processed"""
-        return self._counts[0]
+        # arrays for recording overall statistics, add one for total count
+        n = len(self.quality_criteria) + 1
+        self._counts = np.zeros(n, dtype=np.int64)
+        self._cumulative_counts = np.zeros(n, dtype=np.int64)
 
     def to_table(self, functions=False):
         """
@@ -101,12 +88,12 @@ class QualityQuery(Component):
         from astropy.table import Table
 
         cols = {
-            "criteria": self.criteria_names,
+            "criteria": ["TOTAL"] + self.criteria_names,
             "counts": self._counts,
             "cumulative_counts": self._cumulative_counts,
         }
         if functions:
-            cols["func"] = self.selection_function_strings
+            cols["func"] = ["True"] + self.expressions
         return Table(cols)
 
     def _repr_html_(self):
@@ -117,21 +104,56 @@ class QualityQuery(Component):
         """Print a formatted string representation of the entire table."""
         return self.to_table().pprint_all(show_unit=True, show_dtype=True)
 
-    def __call__(self, value) -> np.ndarray:
+    def __call__(self, **kwargs) -> np.ndarray:
         """
         Test that value passes all cuts
 
         Parameters
         ----------
-        value:
-            the value to pass to each selection function
+        **kwargs:
+            Are passed as locals to evaluate the given expression
 
         Returns
         -------
         np.ndarray:
             array of booleans with results of each selection criterion in order
         """
-        result = np.array(list(map(lambda f: f(value), self._selectors)))
+        # add 1 for total
+        result = np.ones(len(self.quality_criteria) + 1, dtype=bool)
+
+        _evaluate_expression(self._compiled_expressions, result, kwargs)
         self._counts += result.astype(int)
         self._cumulative_counts += result.cumprod()
         return result[1:]  # strip off TOTAL criterion, since redundant
+
+    def get_table_mask(self, table):
+        """
+        Get a boolean mask for the entries that pass the quality checks.
+
+        Parameters
+        ----------
+        table : `~astropy.table.Table`
+            Table with columns matching the expressions used in the
+            `QualityQuery.quality_criteria`.
+
+        Returns
+        -------
+        mask : np.ndarray[bool]
+            Boolean mask of valid entries.
+        """
+        n_criteria = len(self.quality_criteria) + 1
+        result = np.ones((n_criteria, len(table)), dtype=bool)
+
+        _evaluate_expression(self._compiled_expressions, result, table)
+
+        self._counts += np.count_nonzero(result, axis=1)
+        self._cumulative_counts += np.count_nonzero(np.cumprod(result, axis=0), axis=1)
+        return np.all(result, axis=0)
+
+
+def _evaluate_expression(expressions, result, locals):
+    try:
+        for i, expression in enumerate(expressions, start=1):
+            result[i] = eval(expression, ALLOWED_GLOBALS, locals)
+    except Exception:
+        raise QualityCriteriaError("Error evaluating quality query")
diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py
index bf07ea89..81effd10 100644
--- a/ctapipe/image/image_processor.py
+++ b/ctapipe/image/image_processor.py
@@ -46,7 +46,7 @@ class ImageQualityQuery(QualityQuery):
 
     quality_criteria = List(
         default_value=[
-            ("size_greater_0", "lambda image_selected: image_selected.sum() > 0")
+            ("size_greater_0", "image.sum() > 0")
         ],
         help=QualityQuery.quality_criteria.help,
     ).tag(config=True)
@@ -140,7 +140,7 @@ class ImageProcessor(TelescopeComponent):
         image_selected = image[signal_pixels]
 
         # check if image can be parameterized:
-        image_criteria = self.check_image(image_selected)
+        image_criteria = self.check_image(image=image_selected)
         self.log.debug(
             "image_criteria: %s",
             list(zip(self.check_image.criteria_names[1:], image_criteria)),
diff --git a/ctapipe/instrument/subarray.py b/ctapipe/instrument/subarray.py
index 4ac5fd77..91780f3e 100644
--- a/ctapipe/instrument/subarray.py
+++ b/ctapipe/instrument/subarray.py
@@ -137,11 +137,11 @@ class SubarrayDescription:
     def tel_coords(self):
         """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
-        pos_z = np.array([p[2].to("m").value for p in self.positions.values()]) * u.m
+        pos_x = [p[0].to_value(u.m) for p in self.positions.values()]
+        pos_y = [p[1].to_value(u.m) for p in self.positions.values()]
+        pos_z = [p[2].to_value(u.m) for p in self.positions.values()]
 
-        return SkyCoord(x=pos_x, y=pos_y, z=pos_z, frame=GroundFrame())
+        return SkyCoord(x=pos_x, y=pos_y, z=pos_z, unit=u.m, frame=GroundFrame())
 
     @lazyproperty
     def tel_ids(self):
diff --git a/ctapipe/io/datawriter.py b/ctapipe/io/datawriter.py
index fa371877..fedb8406 100644
--- a/ctapipe/io/datawriter.py
+++ b/ctapipe/io/datawriter.py
@@ -173,14 +173,10 @@ class DataWriter(Component):
         help="Store DL1 image parameters if available", default_value=True
     ).tag(config=True)
 
-    write_stereo_shower = Bool(
+    write_showers = Bool(
         help="Store DL2 stereo shower parameters if available", default_value=False
     ).tag(config=True)
 
-    write_mono_shower = Bool(
-        help="Store DL2 mono parameters if available", default_value=False
-    ).tag(config=True)
-
     compression_level = Int(
         help="compression level, 0=None, 9=maximum", default_value=5, min=0, max=9
     ).tag(config=True)
@@ -294,6 +290,14 @@ class DataWriter(Component):
                 containers=[event.index, event.simulation.shower],
             )
 
+            for tel_id, sim in event.simulation.tel.items():
+                table_name = self.table_name(tel_id, self._subarray.tel[tel_id])
+                tel_index = _get_tel_index(event, tel_id)
+                self._writer.write(
+                    f"simulation/event/telescope/impact/{table_name}",
+                    [tel_index, sim.impact],
+                )
+
         if self.write_waveforms:
             self._write_r1_telescope_events(self._writer, event)
 
@@ -304,10 +308,8 @@ class DataWriter(Component):
         self._write_dl1_telescope_events(self._writer, event)
 
         # write DL2 info if requested
-        if self.write_mono_shower:
+        if self.write_showers:
             self._write_dl2_telescope_events(self._writer, event)
-
-        if self.write_stereo_shower:
             self._write_dl2_stereo_event(self._writer, event)
 
     def finish(self):
@@ -342,7 +344,7 @@ class DataWriter(Component):
             data_levels.append(DataLevel.DL1_IMAGES)
         if self.write_parameters:
             data_levels.append(DataLevel.DL1_PARAMETERS)
-        if self.write_stereo_shower or self.write_mono_shower:
+        if self.write_showers:
             data_levels.append(DataLevel.DL2)
         if self.write_raw_waveforms:
             data_levels.append(DataLevel.R0)
@@ -381,8 +383,7 @@ class DataWriter(Component):
         writable_things = [
             self.write_parameters,
             self.write_images,
-            self.write_mono_shower,
-            self.write_stereo_shower,
+            self.write_showers,
             self.write_waveforms,
             self.write_parameters,
         ]
@@ -437,6 +438,7 @@ class DataWriter(Component):
 
         if self._is_simulation:
             writer.exclude("/simulation/event/telescope/images/.*", "true_parameters")
+            writer.exclude("/simulation/event/telescope/images/.*", "impact")
             # no timing information yet for true images
             writer.exclude("/simulation/event/telescope/parameters/.*", r"peak_time_.*")
             writer.exclude("/simulation/event/telescope/parameters/.*", "timing_.*")
@@ -590,11 +592,7 @@ class DataWriter(Component):
     ):
         for tel_id, r1_tel in event.r1.tel.items():
 
-            tel_index = TelEventIndexContainer(
-                obs_id=event.index.obs_id,
-                event_id=event.index.event_id,
-                tel_id=np.int16(tel_id),
-            )
+            tel_index = _get_tel_index(event, tel_id)
             telescope = self._subarray.tel[tel_id]
             table_name = self.table_name(tel_id, str(telescope))
 
@@ -606,11 +604,7 @@ class DataWriter(Component):
     ):
         for tel_id, r0_tel in event.r0.tel.items():
 
-            tel_index = TelEventIndexContainer(
-                obs_id=event.index.obs_id,
-                event_id=event.index.event_id,
-                tel_id=np.int16(tel_id),
-            )
+            tel_index = _get_tel_index(event, tel_id)
             telescope = self._subarray.tel[tel_id]
             table_name = self.table_name(tel_id, str(telescope))
 
@@ -686,8 +680,8 @@ class DataWriter(Component):
                         f"simulation/event/telescope/parameters/{table_name}",
                         [
                             tel_index,
-                            *event.simulation.tel[tel_id].true_parameters.values(),
-                        ],
+                            *event.simulation.tel[tel_id].true_parameters.values()
+                        ]
                     )
 
 
@@ -706,12 +700,7 @@ class DataWriter(Component):
             telescope = self._subarray.tel[tel_id]
             table_name = self.table_name(tel_id, str(telescope))
 
-            tel_index = TelEventIndexContainer(
-                obs_id=event.index.obs_id,
-                event_id=event.index.event_id,
-                tel_id=np.int16(tel_id),
-            )
-
+            tel_index = _get_tel_index(event, tel_id)
             for container_name, algorithm_map in dl2_tel.items():
                 for algorithm, container in algorithm_map.items():
                     name = (
diff --git a/ctapipe/io/hdf5eventsource.py b/ctapipe/io/hdf5eventsource.py
index 48f53418..227eff0d 100644
--- a/ctapipe/io/hdf5eventsource.py
+++ b/ctapipe/io/hdf5eventsource.py
@@ -28,6 +28,7 @@ from ..containers import (
     TelEventIndexContainer,
     TelescopeTriggerContainer,
     R1CameraContainer,
+    TelescopeImpactParameterContainer,
 )
 from .eventsource import EventSource
 from .hdf5tableio import HDF5TableReader
@@ -337,6 +338,7 @@ class HDF5EventSource(EventSource):
                     for tel in self.file_.root.dl1.event.telescope.parameters
                 }
 
+        impact_readers = {}
         if self.is_simulation:
             # simulated shower wide information
             mc_shower_reader = HDF5TableReader(self.file_).read(
@@ -345,6 +347,15 @@ class HDF5EventSource(EventSource):
                 prefixes="true",
             )
             data.simulation = SimulatedEventContainer()
+            if "impact" in self.file_.root.simulation.event.telescope:
+                impact_readers = {
+                    tel.name: self.reader.read(
+                        f"/simulation/event/telescope/impact/{tel.name}",
+                        containers=TelescopeImpactParameterContainer,
+                        prefixes=["true_impact"],
+                    )
+                    for tel in self.file_.root.simulation.event.telescope.impact
+                }
 
         # Setup iterators for the array events
         events = HDF5TableReader(self.file_).read(
@@ -413,6 +424,9 @@ class HDF5EventSource(EventSource):
                 if self.allowed_tels and tel not in self.allowed_tels:
                     continue
 
+                if key in impact_readers:
+                    data.simulation.tel[tel].impact = next(impact_readers[key])
+
                 if DataLevel.R1 in self.datalevels:
                     data.r1.tel[tel] = next(waveform_readers[key])
 
diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py
index 8eb52425..571e337a 100644
--- a/ctapipe/io/simteleventsource.py
+++ b/ctapipe/io/simteleventsource.py
@@ -14,20 +14,21 @@ from ..calib.camera.gainselection import GainSelector
 from ..containers import (
     ArrayEventContainer,
     EventType,
-    SimulatedEventContainer,
-    SimulationConfigContainer,
     SimulatedCameraContainer,
+    SimulatedEventContainer,
     SimulatedShowerContainer,
+    SimulationConfigContainer,
+    TelescopeImpactParameterContainer,
     TelescopePointingContainer,
     TelescopeTriggerContainer,
 )
 from ..coordinates import CameraFrame
 from ..core.traits import (
     Bool,
-    Float,
     CaselessStrEnum,
-    create_class_enum_trait,
+    Float,
     Undefined,
+    create_class_enum_trait,
 )
 from ..instrument import (
     CameraDescription,
@@ -38,7 +39,8 @@ from ..instrument import (
     TelescopeDescription,
 )
 from ..instrument.camera import UnknownPixelShapeWarning
-from ..instrument.guess import unknown_telescope, guess_telescope
+from ..instrument.guess import guess_telescope, unknown_telescope
+from ..reco.impact_distance import shower_impact_distance
 from .datalevels import DataLevel
 from .eventsource import EventSource
 
@@ -440,6 +442,14 @@ class SimTelEventSource(EventSource):
             else:
                 true_image_sums = np.full(self.n_telescopes_original, np.nan)
 
+            if data.simulation.shower is not None:
+                # compute impact distances of the shower to the telescopes
+                impact_distances = shower_impact_distance(
+                    shower_geom=data.simulation.shower, subarray=self.subarray
+                )
+            else:
+                impact_distances = np.full(len(self.subarray), np.nan) * u.m
+
             for tel_id, telescope_event in telescope_events.items():
                 adc_samples = telescope_event.get("adc_samples")
                 if adc_samples is None:
@@ -453,11 +463,22 @@ class SimTelEventSource(EventSource):
                 )
 
                 if data.simulation is not None:
+                    if data.simulation.shower is not None:
+                        impact_container = TelescopeImpactParameterContainer(
+                            distance=impact_distances[self.subarray.tel_index_array[tel_id]],
+                            distance_uncert=0 * u.m,
+                        )
+                    else:
+                        impact_container = TelescopeImpactParameterContainer()
+
+                    impact_container.prefix = "true_impact"
+
                     data.simulation.tel[tel_id] = SimulatedCameraContainer(
                         true_image_sum=true_image_sums[
                             self.telescope_indices_original[tel_id]
                         ],
                         true_image=true_image,
+                        impact=impact_container,
                     )
 
                 data.pointing.tel[tel_id] = self._fill_event_pointing(
diff --git a/ctapipe/io/tableloader.py b/ctapipe/io/tableloader.py
index 37dd55bd..a2fae811 100644
--- a/ctapipe/io/tableloader.py
+++ b/ctapipe/io/tableloader.py
@@ -24,6 +24,7 @@ TRIGGER_TABLE = "/dl1/event/subarray/trigger"
 SHOWER_TABLE = "/simulation/event/subarray/shower"
 TRUE_IMAGES_GROUP = "/simulation/event/telescope/images"
 TRUE_PARAMETERS_GROUP = "/simulation/event/telescope/parameters"
+TRUE_IMPACT_GROUP = "/simulation/event/telescope/impact"
 
 DL2_SUBARRAY_GROUP = "/dl2/event/subarray"
 DL2_TELESCOPE_GROUP = "/dl2/event/telescope"
@@ -152,7 +153,6 @@ class TableLoader(Component):
         if h5file is None and self.input_url is None:
             raise ValueError("Need to specify either input_url or h5file")
 
-        self._should_close = False
         if h5file is None:
             self.h5file = tables.open_file(self.input_url, mode="r")
             self._should_close = True
@@ -283,7 +283,8 @@ class TableLoader(Component):
 
         if self.load_dl2:
             if DL2_TELESCOPE_GROUP in self.h5file:
-                for group_name in self.h5file[DL2_TELESCOPE_GROUP]._v_children:
+                dl2_tel_group = self.h5file.root[DL2_TELESCOPE_GROUP]
+                for group_name in dl2_tel_group._v_children:
                     group_path = f"{DL2_TELESCOPE_GROUP}/{group_name}"
                     group = self.h5file.root[group_path]
 
@@ -315,6 +316,10 @@ class TableLoader(Component):
                 table, self.instrument_table, keys=["tel_id"], join_type="left"
             )
 
+        if self.load_simulated and TRUE_IMPACT_GROUP in self.h5file.root:
+            impacts = self._read_telescope_table(TRUE_IMPACT_GROUP, tel_id)
+            table = _join_telescope_events(table, impacts)
+
         return table
 
     def _read_telescope_events_for_ids(self, tel_ids):
diff --git a/ctapipe/reco/__init__.py b/ctapipe/reco/__init__.py
index d5c89e54..44d33307 100644
--- a/ctapipe/reco/__init__.py
+++ b/ctapipe/reco/__init__.py
@@ -9,9 +9,10 @@ from .hillas_intersection import HillasIntersection
 from .shower_processor import ShowerProcessor
 from .impact import ImPACTReconstructor
 
-
-__all__ = ["Reconstructor",
-           "ShowerProcessor",
-           "HillasReconstructor",
-           "ImPACTReconstructor",
-           "HillasIntersection"]
+__all__ = [
+    "Reconstructor",
+    "ShowerProcessor",
+    "HillasReconstructor",
+    "ImPACTReconstructor",
+    "HillasIntersection",
+]
diff --git a/ctapipe/reco/hillas_reconstructor.py b/ctapipe/reco/hillas_reconstructor.py
index a30e5bfa..f9156bb5 100644
--- a/ctapipe/reco/hillas_reconstructor.py
+++ b/ctapipe/reco/hillas_reconstructor.py
@@ -21,11 +21,11 @@ from ctapipe.coordinates import (
     TiltedGroundFrame,
     project_to_ground,
     MissingFrameAttributeWarning,
+    altaz_to_righthanded_cartesian,
 )
 from astropy.coordinates import (
     SkyCoord,
     AltAz,
-    spherical_to_cartesian,
     cartesian_to_spherical,
 )
 import warnings
@@ -138,8 +138,8 @@ class HillasPlane:
 
         # astropy's coordinates system rotates counter clockwise. Apparently we assume it to
         # be clockwise
-        self.a = np.array(spherical_to_cartesian(1, p1.alt, -p1.az)).ravel()
-        self.b = np.array(spherical_to_cartesian(1, p2.alt, -p2.az)).ravel()
+        self.a = altaz_to_righthanded_cartesian(alt=p1.alt, az=p1.az)
+        self.b = altaz_to_righthanded_cartesian(alt=p2.alt, az=p2.az)
 
         # a and c form an orthogonal basis for the great circle
         # not really necessary since the norm can be calculated
diff --git a/ctapipe/reco/impact_distance.py b/ctapipe/reco/impact_distance.py
new file mode 100644
index 00000000..98b38cde
--- /dev/null
+++ b/ctapipe/reco/impact_distance.py
@@ -0,0 +1,117 @@
+#!/usr/bin/env python3
+"""
+Functions to compute the impact distance from a simulated or reconstructed
+shower axis (Defined by the line from the impact point on the ground in the
+reconstructed sky direction) to each telescope's ground position.
+"""
+
+from typing import Union
+
+import numpy as np
+from astropy import units as u
+from astropy.coordinates import SkyCoord
+
+from ..containers import ReconstructedGeometryContainer, SimulatedShowerContainer
+from ..coordinates import GroundFrame, TiltedGroundFrame, altaz_to_righthanded_cartesian
+from ..instrument.subarray import SubarrayDescription
+
+__all__ = ["shower_impact_distance", "impact_distance"]
+
+
+def impact_distance(point: np.ndarray, direction: np.ndarray, test_points: np.ndarray):
+    """Compute impact distance from a line defined by a point and a direction vector
+    with an array of points
+
+    Parameters
+    ----------
+    point: np.ndarray
+       point defining the line
+    direction: np.ndarray
+       direction vector
+    test_points: np.ndarray
+       array of test points
+    """
+    return np.linalg.norm(
+        np.cross(test_points - point, direction), axis=1
+    ) / np.linalg.norm(direction)
+
+
+def shower_impact_distance(
+    shower_geom: Union[ReconstructedGeometryContainer, SimulatedShowerContainer],
+    subarray: SubarrayDescription,
+):
+    """computes the impact distance of the shower axis to the telescope positions
+
+    Parameters
+    ----------
+    shower_geom: ReconstructedGeometryContainer
+        reconstructed shower geometry, must contain a core and alt/az position
+    subarray: SubarrayDescription
+        SubarrayDescription from which to extract the telescope positions
+
+    Returns
+    -------
+    np.ndarray:
+       array of impact distances to each telescope (packed by tel index)
+    """
+
+    # core position in the ground frame
+    core_position = np.array(
+        [
+            shower_geom.core_x.to_value(u.m),
+            shower_geom.core_y.to_value(u.m),
+            0,
+        ]
+    )
+
+    # sky direction of the shower (defines the shower axis)
+    sky_direction = altaz_to_righthanded_cartesian(
+        alt=shower_geom.alt, az=shower_geom.az
+    )
+
+    # telescope positions in the ground frame
+    telescope_positions = subarray.tel_coords.cartesian.xyz.to_value(u.m).T
+
+    return u.Quantity(
+        impact_distance(
+            point=core_position,
+            direction=sky_direction,
+            test_points=telescope_positions,
+        ),
+        u.m,
+        copy=False,
+    )
+
+
+def shower_impact_distance_with_frames(
+    shower_geom: ReconstructedGeometryContainer, subarray: SubarrayDescription
+):
+    """an alternate and slower implementation for cross-check testing purposes.
+
+    Here we make a TiltedFrame that is aligned with the shower axis and
+    transform both the telescopes and the impact point into that frame. Then the
+    cartesian distance in the X-Y plane should be the impact distance.
+    """
+
+    tilted = TiltedGroundFrame(
+        pointing_direction=SkyCoord(
+            alt=shower_geom.alt, az=shower_geom.az, frame="altaz"
+        )
+    )
+
+    core_pos_shower_tilted = (
+        SkyCoord(
+            x=shower_geom.core_x,
+            y=shower_geom.core_y,
+            z=0 * u.m,
+            frame=GroundFrame(),
+        )
+        .transform_to(tilted)
+        .cartesian.xyz
+    )
+
+    tel_pos_shower_tilted = subarray.tel_coords.transform_to(tilted).cartesian.xyz.T
+
+    # compute cartesian distance in the tilted frame (only x,y coordinates needed):
+    delta = tel_pos_shower_tilted - core_pos_shower_tilted
+    return np.hypot(delta[:, 0], delta[:, 1])
diff --git a/ctapipe/reco/reco_algorithms.py b/ctapipe/reco/reco_algorithms.py
index bd6115c8..19012278 100644
--- a/ctapipe/reco/reco_algorithms.py
+++ b/ctapipe/reco/reco_algorithms.py
@@ -24,9 +24,9 @@ class StereoQualityQuery(QualityQuery):
 
     quality_criteria = List(
         default_value=[
-            ("> 50 phe", "lambda p: p.hillas.intensity > 50"),
-            ("Positive width", "lambda p: p.hillas.width.value > 0"),
-            ("> 3 pixels", "lambda p: p.morphology.num_pixels > 3"),
+            ("> 50 phe", "parameters.hillas.intensity > 50"),
+            ("Positive width", "parameters.hillas.width.value > 0"),
+            ("> 3 pixels", "parameters.morphology.num_pixels > 3"),
         ],
         help=QualityQuery.quality_criteria.help,
     ).tag(config=True)
@@ -64,7 +64,7 @@ class Reconstructor(Component):
         hillas_dict = {
             tel_id: dl1.parameters.hillas
             for tel_id, dl1 in event.dl1.tel.items()
-            if all(self.check_parameters(dl1.parameters))
+            if all(self.check_parameters(parameters=dl1.parameters))
         }
 
         if len(hillas_dict) < 2:
diff --git a/ctapipe/reco/shower_processor.py b/ctapipe/reco/shower_processor.py
index 48f7af59..40c1efa7 100644
--- a/ctapipe/reco/shower_processor.py
+++ b/ctapipe/reco/shower_processor.py
@@ -7,11 +7,12 @@ This processor will be able to process a shower/event in 3 steps:
 - estimation of classification (optional, currently unavailable)
 
 """
-from ctapipe.core import Component
-from ctapipe.core.traits import create_class_enum_trait
-from ctapipe.containers import ArrayEventContainer
-from ctapipe.instrument import SubarrayDescription
-from ctapipe.reco import Reconstructor
+from ..containers import ArrayEventContainer, TelescopeImpactParameterContainer
+from ..core import Component
+from ..core.traits import create_class_enum_trait
+from ..instrument import SubarrayDescription
+from . import Reconstructor
+from .impact_distance import shower_impact_distance
 
 
 class ShowerProcessor(Component):
@@ -23,8 +24,10 @@ class ShowerProcessor(Component):
 
     Input events must already contain dl1 parameters.
     """
+
     reconstructor_type = create_class_enum_trait(
-        Reconstructor, default_value="HillasReconstructor",
+        Reconstructor,
+        default_value="HillasReconstructor",
         help="The stereo geometry reconstructor to be used",
     )
 
@@ -69,3 +72,18 @@ class ShowerProcessor(Component):
         """
         k = self.reconstructor_type
         event.dl2.stereo.geometry[k] = self.reconstructor(event)
+
+        # compute and store the impact parameter for each reconstruction (for
+        # now there is only one, but in the future this should be a loop over
+        # reconstructors)
+
+        # for the stereo reconstructor:
+        impact_distances = shower_impact_distance(
+            shower_geom=event.dl2.stereo.geometry[k], subarray=self.subarray
+        )
+
+        for tel_index, impact_distance in enumerate(impact_distances):
+            tel_id = self.subarray.tel_ids[tel_index]
+            event.dl2.tel[tel_id].impact[k] = TelescopeImpactParameterContainer(
+                distance=impact_distance
+            )
diff --git a/ctapipe/resources/base_config.yaml b/ctapipe/resources/base_config.yaml
index 713bf491..80141764 100644
--- a/ctapipe/resources/base_config.yaml
+++ b/ctapipe/resources/base_config.yaml
@@ -23,7 +23,7 @@ DataWriter:
   overwrite: false # do not overwrite existing files
   write_images: false # store DL1 images
   write_parameters: false # store DL1 parameters
-  write_stereo_shower: false # store DL2 stereo geometry
+  write_showers: false # store DL2 stereo geometry
   write_raw_waveforms: false # write R0 waveforms
   write_waveforms: false # write R1 waveforms
 
@@ -87,11 +87,11 @@ ImageProcessor:
     # filter_function_string] The filter function should take a single value
     # which is the image itself, a 1D np.array.
     quality_criteria:
-      - ["enough_pixels", "lambda im: np.count_nonzero(im) > 2"]
-      - ["enough_charge", "lambda im: im.sum() > 50"]
+      - ["enough_pixels", "np.count_nonzero(image) > 2"]
+      - ["enough_charge", "image.sum() > 50"]
 
 # The ShowerProcessor performs the DL1 to DL2a (reconstructed shower geometry)
-# transition. It is run only if the parameter DataWriter.write_stereo_shower=True.
+# transition. It is run only if the parameter DataWriter.write_showers=True.
 ShowerProcessor:
   # choose between HillasReconstructor and HillasIntersection (two
   # implementations of the standard stereo line-intersection method)
@@ -105,7 +105,7 @@ ShowerProcessor:
       # image parameters, `p` (a `ctapipe.containers.ImageParametersContainer`), as
       # an argument.
       quality_criteria:
-        - [enough intensity, "lambda p: p.hillas.intensity > 50"]
-        - [Positive width, "lambda p: p.hillas.width.value > 0"]
-        - [enough pixels, "lambda p: p.morphology.num_pixels > 3"]
-        - [not clipped, "lambda p: p.leakage.intensity_width_2 < 0.8"]
+        - [enough intensity, "parameters.hillas.intensity > 50"]
+        - [Positive width, "parameters.hillas.width.value > 0"]
+        - [enough pixels, "parameters.morphology.num_pixels > 3"]
+        - [not clipped, "parameters.leakage.intensity_width_2 < 0.8"]
diff --git a/ctapipe/resources/image_modification_config.json b/ctapipe/resources/image_modification_config.json
index 5da4060e..984e7a5e 100644
--- a/ctapipe/resources/image_modification_config.json
+++ b/ctapipe/resources/image_modification_config.json
@@ -8,8 +8,7 @@
         "overwrite": false,
         "write_images": true,
         "write_parameters": true,
-        "write_stereo_shower": false,
-        "write_mono_shower": false,
+        "write_showers": false,
         "transform_image": true,
         "transform_peak_time": true
     },
@@ -36,8 +35,8 @@
         },
         "ImageQualityQuery": {
             "quality_criteria": [
-                ["enough_pixels", "lambda im: np.count_nonzero(im) > 2"],
-                ["enough_charge", "lambda im: im.sum() > 50"]
+                ["enough_pixels", "np.count_nonzero(image) > 2"],
+                ["enough_charge", "image.sum() > 50"]
             ]
         },
         "ImageModifier": {
@@ -66,9 +65,9 @@
         "HillasReconstructor": {
             "StereoQualityQuery": {
                 "quality_criteria": [
-                    ["enough intensity", "lambda p: p.hillas.intensity > 50"],
-                    ["Positive width", "lambda p: p.hillas.width.value > 0"],
-                    ["enough pixels", "lambda p: p.morphology.num_pixels > 3"]
+                    ["enough intensity", "parameters.hillas.intensity > 50"],
+                    ["Positive width", "parameters.hillas.width.value > 0"],
+                    ["enough pixels", "parameters.morphology.num_pixels > 3"]
                 ]
             }
         }
diff --git a/ctapipe/resources/stage1_config.json b/ctapipe/resources/stage1_config.json
index 2fbef2b7..591cf1ce 100644
--- a/ctapipe/resources/stage1_config.json
+++ b/ctapipe/resources/stage1_config.json
@@ -8,8 +8,7 @@
         "overwrite": false,
         "write_images": true,
         "write_parameters": true,
-        "write_stereo_shower": false,
-        "write_mono_shower": false,
+        "write_showers": false,
         "transform_image": true,
         "transform_peak_time": true
     },
@@ -35,8 +34,8 @@
         },
         "ImageQualityQuery": {
             "quality_criteria": [
-                ["enough_pixels", "lambda im: np.count_nonzero(im) > 2"],
-                ["enough_charge", "lambda im: im.sum() > 50"]
+                ["enough_pixels", "np.count_nonzero(image) > 2"],
+                ["enough_charge", "image.sum() > 50"]
             ]
         }
     },
@@ -45,9 +44,9 @@
         "HillasReconstructor": {
             "StereoQualityQuery": {
                 "quality_criteria": [
-                    ["enough intensity", "lambda p: p.hillas.intensity > 50"],
-                    ["Positive width", "lambda p: p.hillas.width.value > 0"],
-                    ["enough pixels", "lambda p: p.morphology.num_pixels > 3"]
+                    ["enough intensity", "parameters.hillas.intensity > 50"],
+                    ["Positive width", "parameters.hillas.width.value > 0"],
+                    ["enough pixels", "parameters.morphology.num_pixels > 3"]
                 ]
             }
         }
diff --git a/ctapipe/resources/stage1_config.toml b/ctapipe/resources/stage1_config.toml
index 4a2b1740..38183d4d 100644
--- a/ctapipe/resources/stage1_config.toml
+++ b/ctapipe/resources/stage1_config.toml
@@ -10,8 +10,7 @@ organization = "YOUR-ORGANIZATION"
 overwrite = false
 write_images = true
 write_parameters = true
-write_stereo_shower = false
-write_mono_shower = false
+write_showers = false
 transform_image = true
 transform_peak_time = true
 
@@ -41,15 +40,15 @@ min_picture_neighbors = [ [ "type", "*", 2,],]
 [ImageProcessor.ImageQualityQuery]
 # These specify which images should be parameterized:
 quality_criteria = [
-    [ "enough_pixels", "lambda im: np.count_nonzero(im) > 2",],
-    [ "enough_charge", "lambda im: im.sum() > 50",]
+    [ "enough_pixels", "np.count_nonzero(image) > 2",],
+    [ "enough_charge", "image.sum() > 50",]
 ]
 
 [ShowerProcessor.ShowerQualityQuery]
 # These specify criteria for telescopes that should be included in stereo
 # reconstruction:
 quality_criteria = [
-    [ "enough intensity", "lambda p: p.hillas.intensity > 50",],
-    [ "Positive width", "lambda p: p.hillas.width.value > 0",],
-    [ "enough pixels", "lambda p: p.morphology.num_pixels > 3",],
+    [ "enough intensity", "parameters.hillas.intensity > 50",],
+    [ "Positive width", "parameters.hillas.width.value > 0",],
+    [ "enough pixels", "parameters.morphology.num_pixels > 3",],
 ]
diff --git a/ctapipe/resources/stage1_config.yaml b/ctapipe/resources/stage1_config.yaml
index 77bac36c..d0413061 100644
--- a/ctapipe/resources/stage1_config.yaml
+++ b/ctapipe/resources/stage1_config.yaml
@@ -12,4 +12,4 @@
 DataWriter:
   write_images: true
   write_parameters: true
-  write_stereo_shower: false
+  write_showers: false
diff --git a/ctapipe/resources/stage2_config.json b/ctapipe/resources/stage2_config.json
index 586faa64..24890065 100644
--- a/ctapipe/resources/stage2_config.json
+++ b/ctapipe/resources/stage2_config.json
@@ -8,8 +8,7 @@
         "overwrite": false,
         "write_images": false,
         "write_parameters": false,
-        "write_stereo_shower": true,
-        "write_mono_shower": true,
+        "write_showers": true,
         "transform_image": true,
         "transform_peak_time": true
     },
@@ -35,8 +34,8 @@
         },
         "ImageQualityQuery": {
             "quality_criteria": [
-                ["enough_pixels", "lambda im: np.count_nonzero(im) > 2"],
-                ["enough_charge", "lambda im: im.sum() > 50"]
+                ["enough_pixels", "np.count_nonzero(image) > 2"],
+                ["enough_charge", "image.sum() > 50"]
             ]
         }
     },
@@ -45,9 +44,9 @@
         "HillasReconstructor": {
             "StereoQualityQuery": {
                 "quality_criteria": [
-                    ["enough intensity", "lambda p: p.hillas.intensity > 50"],
-                    ["Positive width", "lambda p: p.hillas.width.value > 0"],
-                    ["enough pixels", "lambda p: p.morphology.num_pixels > 3"]
+                    ["enough intensity", "parameters.hillas.intensity > 50"],
+                    ["Positive width", "parameters.hillas.width.value > 0"],
+                    ["enough pixels", "parameters.morphology.num_pixels > 3"]
                 ]
             }
         }
diff --git a/ctapipe/resources/stage2_config.yaml b/ctapipe/resources/stage2_config.yaml
index 00ae9725..992fbe52 100644
--- a/ctapipe/resources/stage2_config.yaml
+++ b/ctapipe/resources/stage2_config.yaml
@@ -13,4 +13,4 @@
 DataWriter:
   write_images: false
   write_parameters: false
-  write_stereo_shower: true
+  write_showers: true
diff --git a/ctapipe/resources/training_config.json b/ctapipe/resources/training_config.json
index 35120dcb..cf569d68 100644
--- a/ctapipe/resources/training_config.json
+++ b/ctapipe/resources/training_config.json
@@ -8,8 +8,7 @@
         "overwrite": false,
         "write_images": false,
         "write_parameters": true,
-        "write_stereo_shower": true,
-        "write_mono_shower": true,
+        "write_showers": true,
         "transform_image": true,
         "transform_peak_time": true
     },
@@ -35,8 +34,8 @@
         },
         "ImageQualityQuery": {
             "quality_criteria": [
-                ["enough_pixels", "lambda im: np.count_nonzero(im) > 2"],
-                ["enough_charge", "lambda im: im.sum() > 50"]
+                ["enough_pixels", "np.count_nonzero(image) > 2"],
+                ["enough_charge", "image.sum() > 50"]
             ]
         }
     },
@@ -45,9 +44,9 @@
         "HillasReconstructor": {
             "StereoQualityQuery": {
                 "quality_criteria": [
-                    ["enough intensity", "lambda p: p.hillas.intensity > 50"],
-                    ["Positive width", "lambda p: p.hillas.width.value > 0"],
-                    ["enough pixels", "lambda p: p.morphology.num_pixels > 3"]
+                    ["enough intensity", "parameters.hillas.intensity > 50"],
+                    ["Positive width", "parameters.hillas.width.value > 0"],
+                    ["enough pixels", "parameters.morphology.num_pixels > 3"]
                 ]
             }
         }
diff --git a/ctapipe/resources/training_config.yaml b/ctapipe/resources/training_config.yaml
index 32583424..55ccc3a3 100644
--- a/ctapipe/resources/training_config.yaml
+++ b/ctapipe/resources/training_config.yaml
@@ -14,4 +14,4 @@
 DataWriter:
   write_images: false
   write_parameters: true
-  write_stereo_shower: true
+  write_showers: true
diff --git a/ctapipe/tools/process.py b/ctapipe/tools/process.py
index f9b0bef1..290ae736 100644
--- a/ctapipe/tools/process.py
+++ b/ctapipe/tools/process.py
@@ -112,16 +112,10 @@ class ProcessorTool(Tool):
             "don't store DL1/Event/Telescope parameters in output",
         ),
         **flag(
-            "write-stereo-shower",
-            "DataWriter.write_stereo_shower",
-            "store DL2/Event/Subarray parameters in output",
-            "don't DL2/Event/Subarray parameters in output",
-        ),
-        **flag(
-            "write-mono-shower",
-            "DataWriter.write_mono_shower",
-            "store DL2/Event/Telescope parameters in output",
-            "don't store DL2/Event/Telescope parameters in output",
+            "write-showers",
+            "DataWriter.write_showers",
+            "store DL2/Event parameters in output",
+            "don't DL2/Event parameters in output",
         ),
         **flag(
             "write-index-tables",
@@ -195,7 +189,7 @@ class ProcessorTool(Tool):
         """returns true if we should compute DL2 info"""
         if self.force_recompute_dl2:
             return True
-        return self.write.write_stereo_shower or self.write.write_mono_shower
+        return self.write.write_showers
 
     @property
     def should_compute_dl1(self):
diff --git a/ctapipe/visualization/mpl_array.py b/ctapipe/visualization/mpl_array.py
index 1c3b3e2c..95a83182 100644
--- a/ctapipe/visualization/mpl_array.py
+++ b/ctapipe/visualization/mpl_array.py
@@ -120,7 +120,7 @@ class ArrayDisplay:
                     linewidth=0,
                 )
             )
-        plt.legend(handles=legend_elements)
+        self.axes.legend(handles=legend_elements)
 
         self.add_radial_grid()
 
diff --git a/docs/tutorials/ctapipe_overview.ipynb b/docs/tutorials/ctapipe_overview.ipynb
index a8057e81..f0a638a6 100644
--- a/docs/tutorials/ctapipe_overview.ipynb
+++ b/docs/tutorials/ctapipe_overview.ipynb
@@ -699,7 +699,7 @@
     "\n",
     "f = tempfile.NamedTemporaryFile(suffix='.hdf5')\n",
     "\n",
-    "with DataWriter(source, output_path=f.name, overwrite=True, write_stereo_shower=True) as writer:\n",
+    "with DataWriter(source, output_path=f.name, overwrite=True, write_showers=True) as writer:\n",
     "    \n",
     "    for event in source:\n",
     "        energy = event.simulation.shower.energy\n",
diff --git a/pyproject.toml b/pyproject.toml
index df02cd11..3472edfb 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [tool.black]
 line-length = 88
-target-version = ['py36', 'py37', 'py38']
+target-version = ['py37', 'py38', 'py39', 'py310']
 include = '\.pyi?$'
 exclude = '''
 /(


Test patch:
diff --git a/ctapipe/conftest.py b/ctapipe/conftest.py
index 53d8621b..e90dfa4e 100644
--- a/ctapipe/conftest.py
+++ b/ctapipe/conftest.py
@@ -176,7 +176,7 @@ def dl2_shower_geometry_file(dl2_tmp_path, prod5_gamma_simtel_path):
             f"--input={prod5_gamma_simtel_path}",
             f"--output={output}",
             "--write-images",
-            "--write-stereo-shower",
+            "--write-showers",
             "--max-events=20",
         ]
         assert run_tool(ProcessorTool(), argv=argv, cwd=dl2_tmp_path) == 0
@@ -202,7 +202,7 @@ def dl2_proton_geometry_file(dl2_tmp_path, prod5_proton_simtel_path):
             f"--input={prod5_proton_simtel_path}",
             f"--output={output}",
             "--write-images",
-            "--write-stereo-shower",
+            "--write-showers",
             "--max-events=20",
         ]
         assert run_tool(ProcessorTool(), argv=argv, cwd=dl2_tmp_path) == 0
@@ -228,7 +228,7 @@ def dl2_shower_geometry_file_type(dl2_tmp_path, prod5_gamma_simtel_path):
             f"--input={prod5_gamma_simtel_path}",
             f"--output={output}",
             "--write-images",
-            "--write-stereo-shower",
+            "--write-showers",
             "--max-events=20",
             "--DataWriter.split_datasets_by=tel_type",
         ]
diff --git a/ctapipe/coordinates/tests/test_coordinates.py b/ctapipe/coordinates/tests/test_coordinates.py
index ef75a58d..cec9630d 100644
--- a/ctapipe/coordinates/tests/test_coordinates.py
+++ b/ctapipe/coordinates/tests/test_coordinates.py
@@ -3,10 +3,20 @@ import astropy.units as u
 from astropy.coordinates import SkyCoord, EarthLocation, AltAz
 from astropy.time import Time
 from pytest import approx, raises
+from ctapipe.coordinates import altaz_to_righthanded_cartesian
 
 location = EarthLocation.of_site("Roque de los Muchachos")
 
 
+def test_altaz_to_righthanded_cartesian():
+    """
+    check the handedness of the transform
+    """
+
+    vec = altaz_to_righthanded_cartesian(alt=0 * u.deg, az=90 * u.deg)
+    assert np.allclose(vec, [0, -1, 0])
+
+
 def test_cam_to_nominal():
     from ctapipe.coordinates import CameraFrame, NominalFrame
 
diff --git a/ctapipe/core/tests/test_qualityquery.py b/ctapipe/core/tests/test_qualityquery.py
index 085ea0e4..b70f7dc3 100644
--- a/ctapipe/core/tests/test_qualityquery.py
+++ b/ctapipe/core/tests/test_qualityquery.py
@@ -1,35 +1,39 @@
 """ Tests of Selectors """
 import pytest
+import numpy as np
+from astropy.table import Table
 
 from ctapipe.core.qualityquery import QualityQuery, QualityCriteriaError
 from ctapipe.core.traits import List
 
 
 def test_selector():
-    """ test the functionality of an example Selector subclass"""
+    """test the functionality of an example Selector subclass"""
 
     class ExampleQualityQuery(QualityQuery):
+        """Available variables: x"""
+
         quality_criteria = List(
             default_value=[
-                ("high_enough", "lambda x: x > 3"),
-                ("a_value_not_too_high", "lambda x: x < 100"),
-                ("smallish", "lambda x: x < np.sqrt(100)"),
+                ("high_enough", "x > 3"),
+                ("a_value_not_too_high", "x < 100"),
+                ("smallish", "x < sqrt(100)"),
             ],
         ).tag(config=True)
 
     query = ExampleQualityQuery()
 
-    criteria1 = query(0)  # pass smallish
+    criteria1 = query(x=0)  # pass smallish
     assert len(criteria1) == 3
     assert (criteria1 == [False, True, True]).all()
 
-    criteria2 = query(20)  # pass high_enough + not_too_high
+    criteria2 = query(x=20)  # pass high_enough + not_too_high
     assert (criteria2 == [True, True, False]).all()
 
-    criteria3 = query(200)  # pass high_enough, fail not_too_high
+    criteria3 = query(x=200)  # pass high_enough, fail not_too_high
     assert (criteria3 == [True, False, False]).all()
 
-    criteria4 = query(8)  # pass all
+    criteria4 = query(x=8)  # pass all
     assert (criteria4 == True).all()
 
     tab = query.to_table()
@@ -55,44 +59,67 @@ def test_selector():
     assert tab["cumulative_counts"][2] == 2
     assert tab["cumulative_counts"][3] == 1
 
-    # check that the order is preserved
-    assert query.criteria_names[1] == "high_enough"
-    assert query.criteria_names[2] == "a_value_not_too_high"
-    assert query.criteria_names[3] == "smallish"
-
-    # check we can get back the correct function string:
-    assert query.selection_function_strings[1] == "lambda x: x > 3"
-
-    assert len(query) == 4  # 4 events counted
 
+def test_invalid_input():
+    class ExampleQualityQuery(QualityQuery):
+        """Available variables: x"""
 
-def test_bad_selector():
-    """ ensure failure if a selector function is not a function or can't be evaluated"""
+        quality_criteria = List(
+            default_value=[
+                ("high_enough", "x > 3"),
+                ("a_value_not_too_high", "x < 100"),
+                ("smallish", "x < np.sqrt(100)"),
+            ],
+        ).tag(config=True)
 
+    query = ExampleQualityQuery()
     with pytest.raises(QualityCriteriaError):
-        s = QualityQuery(
-            quality_criteria=[
-                ("high_enough", "lambda x: x > 3"),
-                ("not_good", "3"),
-                ("smallish", "lambda x: x < 10"),
-            ]
-        )
-        assert s
+        query(y=5)
+
 
+def test_bad_selector():
+    """ensure failure if a selector function is not a function or can't be evaluated"""
+
+    query = QualityQuery(
+        quality_criteria=[
+            ("high_enough", "x > 3"),
+            ("not_good", "foo"),
+            ("smallish", "x < 10"),
+        ]
+    )
     with pytest.raises(QualityCriteriaError):
-        s = QualityQuery(
-            quality_criteria=[
-                ("high_enough", "lambda x: x > 3"),
-                ("not_good", "x == 3"),
-                ("smallish", "lambda x: x < 10"),
-            ]
-        )
-        assert s
+        query(x=5)
 
     # ensure we can't run arbitrary code.
     # try to construct something that is not in the
     # ALLOWED_GLOBALS list, but which is imported in selector.py
     # and see if it works in a function
-    with pytest.raises(NameError):
-        s = QualityQuery(quality_criteria=[("dangerous", "lambda x: Component()")])
-        s(10)
+    with pytest.raises(QualityCriteriaError):
+        query = QualityQuery(quality_criteria=[("dangerous", "Component()")])
+        query(x=10)
+
+    # test we only support expressions, not statements
+    with pytest.raises(QualityCriteriaError):
+        query = QualityQuery(
+            quality_criteria=[("dangerous", "import numpy; np.array([])")]
+        )
+
+
+def test_table_mask():
+    """Test getting a mask for a whole table"""
+    query = QualityQuery(
+        quality_criteria=[
+            ("foo", "(x**2 + y**2) < 1.0"),
+            ("bar", "x < 0.5"),
+        ],
+    )
+
+    table = Table({"x": [1.0, 0.2, -0.5, 0.6, 0.7], "y": [0.0, 0.5, 1.0, 0.2, 0.1]})
+
+    mask = query.get_table_mask(table)
+    assert len(mask) == len(table)
+    assert mask.dtype == np.bool_
+    np.testing.assert_equal(mask, [False, True, False, False, False])
+    stats = query.to_table()
+    np.testing.assert_equal(stats["counts"], [5, 3, 2])
+    np.testing.assert_equal(stats["cumulative_counts"], [5, 3, 1])
diff --git a/ctapipe/io/tests/conftest.py b/ctapipe/io/tests/conftest.py
index 8c818a23..b1499606 100644
--- a/ctapipe/io/tests/conftest.py
+++ b/ctapipe/io/tests/conftest.py
@@ -21,8 +21,7 @@ def r1_hdf5_file(prod5_proton_simtel_path, r1_path):
         output_path=path,
         write_parameters=False,
         write_images=False,
-        write_stereo_shower=False,
-        write_mono_shower=False,
+        write_showers=False,
         write_raw_waveforms=False,
         write_waveforms=True,
     )
diff --git a/ctapipe/io/tests/test_datawriter.py b/ctapipe/io/tests/test_datawriter.py
index 203f7810..a4242b0c 100644
--- a/ctapipe/io/tests/test_datawriter.py
+++ b/ctapipe/io/tests/test_datawriter.py
@@ -60,8 +60,7 @@ def test_write(tmpdir: Path):
         output_path=output_path,
         write_parameters=False,
         write_images=True,
-        write_stereo_shower=True,
-        write_mono_shower=True,
+        write_showers=True,
         write_raw_waveforms=True,
         write_waveforms=True,
     ) as writer:
@@ -149,8 +148,7 @@ def test_roundtrip(tmpdir: Path):
         transform_peak_time=True,
         peak_time_dtype="int16",
         peak_time_scale=100,
-        write_stereo_shower=True,
-        write_mono_shower=True,
+        write_showers=True,
     ) as write:
         write.log.level = logging.DEBUG
         for event in source:
diff --git a/ctapipe/io/tests/test_table_loader.py b/ctapipe/io/tests/test_table_loader.py
index b5b3e7bc..a22a3814 100644
--- a/ctapipe/io/tests/test_table_loader.py
+++ b/ctapipe/io/tests/test_table_loader.py
@@ -83,6 +83,7 @@ def test_load_simulated(test_file):
 
         table = table_loader.read_telescope_events([25])
         assert "true_energy" in table.colnames
+        assert "true_impact_distance" in table.colnames
 
 
 def test_true_images(test_file):
@@ -150,8 +151,9 @@ def test_read_telescope_events_type(test_file_dl2):
         assert "HillasReconstructor_alt" in table.colnames
         assert "true_energy" in table.colnames
         assert "true_image" in table.colnames
-        assert set(table["tel_id"].data).issubset([25, 125, 130])
+        assert {25, 125, 127}.issubset(set(table["tel_id"].data))
         assert "equivalent_focal_length" in table.colnames
+        assert "HillasReconstructor_impact_distance" in table.colnames
 
 
 def test_read_telescope_events_by_type(test_file_dl2):
diff --git a/ctapipe/reco/tests/test_HillasReconstructor.py b/ctapipe/reco/tests/test_HillasReconstructor.py
index 0b594cab..b926184f 100644
--- a/ctapipe/reco/tests/test_HillasReconstructor.py
+++ b/ctapipe/reco/tests/test_HillasReconstructor.py
@@ -208,7 +208,7 @@ def test_CameraFrame_against_TelescopeFrame(filename):
     config = Config({
         "StereoQualityQuery": {
             "quality_criteria": [
-                ("valid_width", "lambda p: p.hillas.width.value > 0"),
+                ("valid_width", "parameters.hillas.width.value > 0"),
             ]
         }
     })
diff --git a/ctapipe/reco/tests/test_impact_distance.py b/ctapipe/reco/tests/test_impact_distance.py
new file mode 100644
index 00000000..960006b4
--- /dev/null
+++ b/ctapipe/reco/tests/test_impact_distance.py
@@ -0,0 +1,115 @@
+#!/usr/bin/env python3
+
+"""
+tests for impact distance computations
+"""
+
+import numpy as np
+from astropy import units as u
+from ctapipe.containers import ReconstructedGeometryContainer
+from ctapipe.instrument.subarray import SubarrayDescription
+from ctapipe.reco.impact_distance import (
+    impact_distance,
+    shower_impact_distance,
+    shower_impact_distance_with_frames,
+)
+
+
+def test_impact_distance():
+    """just test a few limit cases. the more detailed tests are done in
+    test_shower_impact_distance_off_axis"""
+
+    # test on axis (coming from  zenith onto an array)
+    point = np.array([0, 0, 0])
+    direction = np.array([0, 0, 1])  # coming perfectly down from zenith
+    test_points = np.array([[0, 0, 0], [0, 1, 0], [0, 2, 100], [0, 3, -100]])
+    impacts = impact_distance(point=point, direction=direction, test_points=test_points)
+
+    expected_impacts = np.array([0, 1, 2, 3])
+    assert np.allclose(impacts, expected_impacts)
+
+    # test at perfectly off-axis (coming in parallel to the array from x-direction)
+
+    direction = np.array([1, 0, 0])  # coming from the x-direction
+    test_points = np.array(
+        [[0, 0, 0], [0, 1, 0], [0, 2, 0], [0, 3, 0]]
+    )  # points along y axis
+    impacts = impact_distance(point=point, direction=direction, test_points=test_points)
+    expected_impacts = np.array(
+        [0, 1, 2, 3]
+    )  # expect the impacts to be equal to the y's
+    assert np.allclose(impacts, expected_impacts)
+
+    # test at perfectly off-axis (coming in parallel to the array from y-direction)
+
+    direction = np.array([0, 1, 0])  # coming from the x-direction
+    test_points = np.array(
+        [[0, 0, 0], [0, 1, 0], [0, 2, 0], [0, 3, 0]]
+    )  # points along y axis
+    impacts = impact_distance(point=point, direction=direction, test_points=test_points)
+    expected_impacts = np.array(
+        [0, 0, 0, 0]
+    )  # expect the impacts to be equal to the y's
+    assert np.allclose(impacts, expected_impacts)
+
+
+def test_shower_impact_distance():
+    """test several boundary cases using the function that takes a Subarray and
+    Container
+    """
+
+    sub = SubarrayDescription(
+        name="test",
+        tel_positions={1: [0, 0, 0] * u.m, 2: [0, 1, 0] * u.m, 3: [0, 2, 0] * u.m},
+        tel_descriptions={1: None, 2: None, 3: None},
+    )
+
+    # coming from zenith to the center of the array, the impact distance should
+    # be the cartesian distance
+    shower_geom = ReconstructedGeometryContainer(
+        core_x=0 * u.m, core_y=0 * u.m, alt=90 * u.deg, az=0 * u.deg
+    )
+    impact_distances = shower_impact_distance(shower_geom=shower_geom, subarray=sub)
+    assert np.allclose(impact_distances, [0, 1, 2] * u.m)
+
+    # alt=0  az=0 should be aligned to x-axis (north in ground frame)
+    # therefore the distances should also be just the y-offset of the telecope
+    shower_geom = ReconstructedGeometryContainer(
+        core_x=0 * u.m, core_y=0 * u.m, alt=0 * u.deg, az=0 * u.deg
+    )
+    impact_distances = shower_impact_distance(shower_geom=shower_geom, subarray=sub)
+    assert np.allclose(impact_distances, [0, 1, 2] * u.m)
+
+    # alt=0  az=90 should be aligned to y-axis (east in ground frame)
+    # therefore the distances should also be just the x-offset of the telecope
+    shower_geom = ReconstructedGeometryContainer(
+        core_x=0 * u.m, core_y=0 * u.m, alt=0 * u.deg, az=90 * u.deg
+    )
+    impact_distances = shower_impact_distance(shower_geom=shower_geom, subarray=sub)
+    assert np.allclose(impact_distances, [0, 0, 0] * u.m)
+
+
+def test_compare_3d_and_frame_impact_distance(
+    example_subarray: SubarrayDescription,
+):
+    """Test another (slower) way of computing the impact distance, using Frames
+    and compare to the implemented method.
+    """
+
+    for alt in [90, 0, 50, 60] * u.deg:
+        for az in [0, 90, 45, 270, 360] * u.deg:
+
+            shower_geom = ReconstructedGeometryContainer(
+                core_x=0 * u.m, core_y=0 * u.m, alt=alt, az=az
+            )
+            impact_distances_3d = shower_impact_distance(
+                shower_geom=shower_geom, subarray=example_subarray
+            )
+
+            impact_distances_frame = shower_impact_distance_with_frames(
+                shower_geom=shower_geom, subarray=example_subarray
+            )
+
+            assert np.allclose(
+                impact_distances_frame, impact_distances_3d
+            ), f"failed at {alt=} {az=}"
diff --git a/ctapipe/reco/tests/test_shower_processor.py b/ctapipe/reco/tests/test_shower_processor.py
index cd49fff7..fe3b87dd 100644
--- a/ctapipe/reco/tests/test_shower_processor.py
+++ b/ctapipe/reco/tests/test_shower_processor.py
@@ -45,7 +45,7 @@ def test_shower_processor_geometry(example_event, example_subarray, reconstructo
 
     # Increase some quality cuts and check that we get defaults
     config[reconstructor_type].StereoQualityQuery.quality_criteria = [
-        ("> 500 phes", "lambda p: p.hillas.intensity > 500")
+        ("> 500 phes", "parameters.hillas.intensity > 500")
     ]
 
     process_shower = ShowerProcessor(
diff --git a/ctapipe/tools/tests/test_process.py b/ctapipe/tools/tests/test_process.py
index c9f4fbb5..b041df41 100644
--- a/ctapipe/tools/tests/test_process.py
+++ b/ctapipe/tools/tests/test_process.py
@@ -63,8 +63,8 @@ def test_multiple_configs(dl1_image_file):
 
     # ensure the overwriting works (base config has this option disabled)
     assert (
-        tool.get_current_config()["ProcessorTool"]["DataWriter"]["write_stereo_shower"]
-        == True
+        tool.get_current_config()["ProcessorTool"]["DataWriter"]["write_showers"]
+        is True
     )