Repository: cta-observatory/ctapipe
Pull request #1935 (2022-06-14 18:39:02)

Chunked table loader
Needs #1934 


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 6c869260..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,
@@ -494,7 +507,7 @@ class SimulatedCameraContainer(Container):
     container_prefix = ""
 
     true_image_sum = Field(
-        np.nan, "Sum of the detected Cherenkov photons in the camera"
+        np.int32(-1), "Total number of detected Cherenkov photons in the camera"
     )
     true_image = Field(
         None,
@@ -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/logging.py b/ctapipe/core/logging.py
index 5ac27eee..68ffdb4e 100644
--- a/ctapipe/core/logging.py
+++ b/ctapipe/core/logging.py
@@ -68,10 +68,15 @@ def recursive_update(d1, d2, copy=False):
 
 
 def create_logging_config(
-    log_level, log_file, log_file_level, log_config: dict, quiet: bool
+    log_level,
+    log_file,
+    log_file_level,
+    log_config: dict,
+    quiet: bool,
+    module="ctapipe",
 ):
     """Update logging level for console and file according to CLI arguments."""
-    config = recursive_update(DEFAULT_LOGGING, log_config)
+    config = recursive_update(get_default_logging(module), log_config)
 
     if quiet:
         config["handlers"]["console"] = {"class": "logging.NullHandler"}
@@ -89,15 +94,13 @@ def create_logging_config(
             }
         }
         config["handlers"].update(file_handler)
-        config["loggers"]["ctapipe"]["handlers"].append("ctapipe-file")
+        config["loggers"][module]["handlers"].append("ctapipe-file")
 
         # level of logger must be at least that of all their handlers
-        config["loggers"]["ctapipe"]["level"] = get_lower_level(
-            log_level, log_file_level
-        )
+        config["loggers"][module]["level"] = get_lower_level(log_level, log_file_level)
 
     else:
-        config["loggers"]["ctapipe"]["level"] = log_level
+        config["loggers"][module]["level"] = log_level
 
     return config
 
@@ -112,32 +115,33 @@ def get_lower_level(l0, l1):
     return l0 if l0 < l1 else l1
 
 
-DEFAULT_LOGGING = {
-    "version": 1,
-    "disable_existing_loggers": False,
-    "formatters": {
-        "file": {"()": PlainFormatter, "fmt": DEFAULT_LOGGING_FORMAT},
-        "console": {"()": ColoredFormatter, "fmt": DEFAULT_LOGGING_FORMAT},
-    },
-    "handlers": {
-        "ctapipe-console": {
-            "class": "logging.StreamHandler",
-            "formatter": "console",
-            "stream": "ext://sys.stderr",
-            "level": "NOTSET",
+def get_default_logging(module="ctapipe"):
+    return {
+        "version": 1,
+        "disable_existing_loggers": False,
+        "formatters": {
+            "file": {"()": PlainFormatter, "fmt": DEFAULT_LOGGING_FORMAT},
+            "console": {"()": ColoredFormatter, "fmt": DEFAULT_LOGGING_FORMAT},
         },
-        "console": {
-            "class": "logging.StreamHandler",
-            "formatter": "console",
-            "stream": "ext://sys.stderr",
-            "level": "NOTSET",
+        "handlers": {
+            "ctapipe-console": {
+                "class": "logging.StreamHandler",
+                "formatter": "console",
+                "stream": "ext://sys.stderr",
+                "level": "NOTSET",
+            },
+            "console": {
+                "class": "logging.StreamHandler",
+                "formatter": "console",
+                "stream": "ext://sys.stderr",
+                "level": "NOTSET",
+            },
         },
-    },
-    "loggers": {
-        "ctapipe": {
-            "level": "WARN",
-            "handlers": ["ctapipe-console"],
-            "propagate": False,
-        }
-    },
-}
+        "loggers": {
+            module: {
+                "level": "WARN",
+                "handlers": ["ctapipe-console"],
+                "propagate": False,
+            }
+        },
+    }
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..fe3c3458 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):
@@ -115,23 +102,58 @@ class QualityQuery(Component):
 
     def __str__(self):
         """Print a formatted string representation of the entire table."""
-        return self.to_table().pprint_all(show_unit=True, show_dtype=True)
+        return str(self.to_table())
 
-    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/core/tool.py b/ctapipe/core/tool.py
index b7c23718..ab0fa1de 100644
--- a/ctapipe/core/tool.py
+++ b/ctapipe/core/tool.py
@@ -24,7 +24,7 @@ from traitlets.config import Application, Config, Configurable
 from .. import __version__ as version
 from . import Provenance
 from .component import Component
-from .logging import DEFAULT_LOGGING, ColoredFormatter, create_logging_config
+from .logging import ColoredFormatter, create_logging_config
 from .traits import Bool, Dict, Enum, Path
 
 __all__ = ["Tool", "ToolConfigurationError"]
@@ -133,8 +133,6 @@ class Tool(Application):
         trait=Path(
             exists=True,
             directory_ok=False,
-            allow_none=True,
-            default_value=None,
             help=(
                 "List of configuration files with parameters to load "
                 "in addition to command-line parameters. "
@@ -146,7 +144,7 @@ class Tool(Application):
         )
     ).tag(config=True)
 
-    log_config = Dict(default_value=DEFAULT_LOGGING).tag(config=True)
+    log_config = Dict(default_value={}).tag(config=True)
     log_file = Path(
         default_value=None,
         exists=None,
@@ -197,8 +195,8 @@ class Tool(Application):
 
         # tools defined in other modules should have those modules as base
         # logging name
-        module_name = self.__class__.__module__.split(".")[0]
-        self.log = logging.getLogger(f"{module_name}.{self.name}")
+        self.module_name = self.__class__.__module__.split(".")[0]
+        self.log = logging.getLogger(f"{self.module_name}.{self.name}")
         self.trait_warning_handler = CollectTraitWarningsHandler()
         self.update_logging_config()
 
@@ -233,7 +231,6 @@ class Tool(Application):
         path: Union[str, pathlib.Path]
             config file to load. [yaml, toml, json, py] formats are supported
         """
-
         path = pathlib.Path(path)
 
         if path.suffix in [".yaml", ".yml"]:
@@ -259,6 +256,7 @@ class Tool(Application):
             log_file_level=self.log_file_level,
             log_config=self.log_config,
             quiet=self.quiet,
+            module=self.module_name,
         )
 
         logging.config.dictConfig(cfg)
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/camera/geometry.py b/ctapipe/instrument/camera/geometry.py
index 2dc6c169..09a5eb8d 100644
--- a/ctapipe/instrument/camera/geometry.py
+++ b/ctapipe/instrument/camera/geometry.py
@@ -262,7 +262,7 @@ class CameraGeometry:
         scale = np.sqrt(uv_trans_x[0] ** 2 + uv_trans_y[0] ** 2) / cam.pix_x.unit
         trans_x = getattr(trans, trans_x_name)
         trans_y = getattr(trans, trans_y_name)
-        pix_area = (cam.pix_area * scale ** 2).to(trans_x.unit ** 2)
+        pix_area = (cam.pix_area * scale**2).to(trans_x.unit**2)
 
         return CameraGeometry(
             camera_name=cam.camera_name,
@@ -282,10 +282,10 @@ class CameraGeometry:
         return hash(
             (
                 self.camera_name,
-                self.pix_x[0].value,
-                self.pix_y[0].value,
+                round(self.pix_x[0].value, 3),
+                round(self.pix_y[0].value, 3),
                 self.pix_type,
-                self.pix_rotation.deg,
+                round(self.pix_rotation.deg, 3),
             )
         )
 
@@ -321,7 +321,7 @@ class CameraGeometry:
         if pix_type == PixelShape.HEXAGON:
             area = 2 * np.sqrt(3) * (dist / 2) ** 2
         elif pix_type == PixelShape.SQUARE:
-            area = dist ** 2
+            area = dist**2
         else:
             raise KeyError("unsupported pixel type")
 
@@ -582,7 +582,7 @@ class CameraGeometry:
         return CameraGeometry.from_table(table)
 
     def to_table(self):
-        """ convert this to an `astropy.table.Table` """
+        """convert this to an `astropy.table.Table`"""
         # currently the neighbor list is not supported, since
         # var-length arrays are not supported by astropy.table.Table
         t = Table(
@@ -771,18 +771,18 @@ class CameraGeometry:
             [
                 x,
                 y,
-                x ** 2,
+                x**2,
                 x * y,
-                y ** 2,
-                x ** 3,
-                x ** 2 * y,
-                x * y ** 2,
-                y ** 3,
-                x ** 4,
-                x ** 3 * y,
-                x ** 2 * y ** 2,
-                x * y ** 3,
-                y ** 4,
+                y**2,
+                x**3,
+                x**2 * y,
+                x * y**2,
+                y**3,
+                x**4,
+                x**3 * y,
+                x**2 * y**2,
+                x * y**3,
+                y**4,
             ]
         )
 
@@ -819,7 +819,7 @@ class CameraGeometry:
         self.cam_rotation = Angle(0, unit=u.deg)
 
     def info(self, printer=print):
-        """ print detailed info about this camera """
+        """print detailed info about this camera"""
         printer(f'CameraGeometry: "{self}"')
         printer("   - num-pixels: {}".format(len(self.pix_id)))
         printer(f"   - pixel-type: {self.pix_type}")
diff --git a/ctapipe/instrument/subarray.py b/ctapipe/instrument/subarray.py
index 4ac5fd77..853ddfb6 100644
--- a/ctapipe/instrument/subarray.py
+++ b/ctapipe/instrument/subarray.py
@@ -11,7 +11,7 @@ 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.table import QTable, Table, join
 from astropy.utils import lazyproperty
 
 from .. import __version__ as CTAPIPE_VERSION
@@ -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):
@@ -238,6 +238,24 @@ class SubarrayDescription:
             which table to generate (subarray or optics)
         """
 
+        if kind == 'joined':
+            table = self.to_table()
+            optics = self.to_table(kind="optics")
+            optics["optics_index"] = np.arange(len(optics))
+            optics.remove_columns(["name", "description", "type"])
+            table = join(
+                table,
+                optics,
+                keys="optics_index",
+                # conflicts for TAB_VER, TAB_TYPE, not needed here, ignore
+                metadata_conflicts="silent",
+            )
+
+            table.remove_columns(["optics_index", "camera_index"])
+            table.add_index('tel_id')
+            return table
+
+
         meta = {
             "ORIGIN": "ctapipe.instrument.SubarrayDescription",
             "SUBARRAY": self.name,
diff --git a/ctapipe/io/astropy_helpers.py b/ctapipe/io/astropy_helpers.py
index aca84323..bcd2687c 100644
--- a/ctapipe/io/astropy_helpers.py
+++ b/ctapipe/io/astropy_helpers.py
@@ -76,6 +76,8 @@ def read_table(
         # methods
         path = os.path.join("/", path)
         table = h5file.get_node(path)
+        if not isinstance(table, tables.Table):
+            raise IOError(f"Node {path} is a {table.__class__.__name__}, must be a Table")
         transforms, descriptions, meta = _parse_hdf5_attrs(table)
 
         if condition is None:
diff --git a/ctapipe/io/datawriter.py b/ctapipe/io/datawriter.py
index 712e733a..1a4a2fa5 100644
--- a/ctapipe/io/datawriter.py
+++ b/ctapipe/io/datawriter.py
@@ -6,11 +6,12 @@ Class to write DL1 (a,b) and DL2 (a) data from an event stream
 
 import pathlib
 from collections import defaultdict
-from traitlets import Instance
+from typing import List
 
 import numpy as np
 import tables
 from astropy import units as u
+from traitlets import Dict, Instance
 
 from ..containers import (
     ArrayEventContainer,
@@ -60,8 +61,14 @@ PROV = Provenance()
 
 
 def write_reference_metadata_headers(
-    obs_ids, subarray, writer, is_simulation, data_levels, contact_info
-):
+    obs_ids: List[int],
+    subarray: SubarrayDescription,
+    writer: "DataWriter",
+    is_simulation: bool,
+    data_levels,
+    contact_info: meta.Contact,
+    instrument_info: meta.Instrument,
+) -> None:
     """
     Attaches Core Provenence headers to an output HDF5 file.
     Right now this is hard-coded for use with the ctapipe-process tool
@@ -79,6 +86,10 @@ def write_reference_metadata_headers(
     data_levels: List[DataLevel]
         list of data levels that were requested/generated
         (e.g. from `DataWriter.datalevels`)
+    contact_info: meta.Contact
+        contact metadata
+    instrument_info: meta.Instrument
+        instrument metadata
     """
     activity = PROV.current_activity.provenance
     category = "Sim" if is_simulation else "Other"
@@ -101,15 +112,12 @@ def write_reference_metadata_headers(
             id_=",".join(str(x) for x in obs_ids),
         ),
         activity=meta.Activity.from_provenance(activity),
-        instrument=meta.Instrument(
-            site="Other",  # need a way to detect site...
-            class_="Subarray",
-            type_="unknown",
-            version="unknown",
-            id_=subarray.name,
-        ),
+        instrument=instrument_info,
     )
 
+    if reference.instrument.id_ == "unspecified":
+        reference.instrument.id_ = subarray.name
+
     headers = reference.to_dict()
     meta.write_to_hdf5(headers, writer.h5file)
 
@@ -133,6 +141,17 @@ class DataWriter(Component):
 
     # pylint: disable=too-many-instance-attributes
     contact_info = Instance(meta.Contact, kw={}).tag(config=True)
+    instrument_info = Instance(meta.Instrument, kw={}).tag(config=True)
+
+    context_metadata = Dict(
+        help=(
+            "Additional metadata keywords and values that describe this data. "
+            "This should be a dictionary where the keys will be appended to the "
+            "CONTEXT section of the output file's attributes. Keys can be hierarchical "
+            "by using a space between each level, e.g. `SIMULATION PRODUCTION` "
+            "would make a key PRODUCTION grouped under the key SIMULATION"
+        )
+    ).tag(config=True)
 
     output_path = Path(
         help="output filename", default_value=pathlib.Path("events.dl1.h5")
@@ -154,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)
@@ -226,6 +241,7 @@ class DataWriter(Component):
 
         self.event_source = event_source
         self.contact_info = meta.Contact(parent=self)
+        self.instrument_info = meta.Instrument(parent=self)
 
         self._at_least_one_event = False
         self._is_simulation = event_source.is_simulation
@@ -274,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)
 
@@ -284,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):
@@ -306,7 +328,11 @@ class DataWriter(Component):
                 is_simulation=self._is_simulation,
                 data_levels=self.datalevels,
                 contact_info=self.contact_info,
+                instrument_info=self.instrument_info,
             )
+
+            self._write_context_metadata_headers()
+
             self._writer.close()
             self._writer = None
 
@@ -318,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)
@@ -357,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,
         ]
@@ -388,12 +413,6 @@ class DataWriter(Component):
             transform=tr_tel_list_to_mask,
         )
 
-        # avoid some warnings about unwritable columns (which here are just
-        # sub-containers)
-        writer.exclude("dl1/event/subarray/trigger", "tel")
-        writer.exclude("dl1/monitoring/subarray/pointing", "tel")
-        writer.exclude("/dl1/event/telescope/images/.*", "parameters")
-
         # currently the trigger info is used for the event time, but we dont'
         # want the other bits of the trigger container in the pointing or other
         # montitoring containers
@@ -405,18 +424,18 @@ class DataWriter(Component):
         writer.exclude("/dl1/monitoring/telescope/pointing/.*", "trigger_pixels")
         writer.exclude("/dl1/monitoring/event/pointing/.*", "event_type")
 
-        if self.write_parameters is False:
+        if not self.write_images:
+            writer.exclude("/simulation/event/telescope/images/.*", "true_image")
+
+        if not self.write_parameters:
             writer.exclude("/dl1/event/telescope/images/.*", "image_mask")
 
         if self._is_simulation:
-            writer.exclude("/simulation/event/telescope/images/.*", "true_parameters")
             # no timing information yet for true images
             writer.exclude("/simulation/event/telescope/parameters/.*", r"peak_time_.*")
             writer.exclude("/simulation/event/telescope/parameters/.*", "timing_.*")
-            writer.exclude("/simulation/event/subarray/shower", "true_tel")
 
         # Set up transforms
-
         if self.transform_image:
             transform = FixedPointColumnTransform(
                 scale=self.image_scale,
@@ -563,11 +582,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))
 
@@ -579,11 +594,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))
 
@@ -626,46 +637,44 @@ class DataWriter(Component):
 
             table_name = self.table_name(tel_id, str(telescope))
 
-            has_sim_camera = self._is_simulation and (
-                tel_id in event.simulation.tel
-                and event.simulation.tel[tel_id].true_image is not None
-            )
-
             if self.write_parameters:
                 writer.write(
                     table_name=f"dl1/event/telescope/parameters/{table_name}",
                     containers=[tel_index, *dl1_camera.parameters.values()],
                 )
 
-                if has_sim_camera:
-                    writer.write(
-                        f"simulation/event/telescope/parameters/{table_name}",
-                        [
-                            tel_index,
-                            *event.simulation.tel[tel_id].true_parameters.values(),
-                        ],
-                    )
-
             if self.write_images:
                 if dl1_camera.image is None:
                     raise ValueError(
                         "DataWriter.write_images is True but event does not contain image"
                     )
-                # note that we always write the image, even if the image quality
-                # criteria are not met (those are only to determine if the parameters
-                # can be computed).
-                self.log.debug("WRITING IMAGES")
+
                 writer.write(
                     table_name=f"dl1/event/telescope/images/{table_name}",
                     containers=[tel_index, dl1_camera],
                 )
 
-                if self._is_simulation and has_sim_camera:
+            if self._is_simulation:
+                # always write this, so that at least the sum is included
+                writer.write(
+                    f"simulation/event/telescope/images/{table_name}",
+                    [tel_index, event.simulation.tel[tel_id]],
+                )
+
+                has_sim_image = (
+                    tel_id in event.simulation.tel
+                    and event.simulation.tel[tel_id].true_image is not None
+                )
+                if self.write_parameters and has_sim_image:
                     writer.write(
-                        f"simulation/event/telescope/images/{table_name}",
-                        [tel_index, event.simulation.tel[tel_id]],
+                        f"simulation/event/telescope/parameters/{table_name}",
+                        [
+                            tel_index,
+                            *event.simulation.tel[tel_id].true_parameters.values()
+                        ]
                     )
 
+
     def _write_dl2_telescope_events(
         self, writer: TableWriter, event: ArrayEventContainer
     ):
@@ -681,12 +690,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 = (
@@ -750,3 +754,23 @@ class DataWriter(Component):
                 )
 
         self._generate_table_indices(self._writer.h5file, "/dl1/event/subarray")
+
+    def _write_context_metadata_headers(self):
+        """write out any user-defined metadata in the context_metadata field to the
+        headers.
+
+        This will create a set of headers that start with CONTEXT <KEY> = value
+
+        Keys can be hierarchical separated by spaces
+
+        """
+
+        # first append CONTEXT to each key in the dict
+
+        context_dict = {}
+
+        for key, value in self.context_metadata.items():
+            key = " ".join(["CONTEXT", key])
+            context_dict[key] = value
+
+        meta.write_to_hdf5(context_dict, self._writer.h5file)
diff --git a/ctapipe/io/hdf5eventsource.py b/ctapipe/io/hdf5eventsource.py
index 48f53418..fa4aaeab 100644
--- a/ctapipe/io/hdf5eventsource.py
+++ b/ctapipe/io/hdf5eventsource.py
@@ -17,6 +17,8 @@ from ..containers import (
     IntensityStatisticsContainer,
     LeakageContainer,
     MorphologyContainer,
+    ParticleClassificationContainer,
+    ReconstructedEnergyContainer,
     SimulationConfigContainer,
     SimulatedShowerContainer,
     SimulatedEventContainer,
@@ -28,12 +30,15 @@ from ..containers import (
     TelEventIndexContainer,
     TelescopeTriggerContainer,
     R1CameraContainer,
+    TelescopeImpactParameterContainer,
+    ReconstructedGeometryContainer,
 )
 from .eventsource import EventSource
 from .hdf5tableio import HDF5TableReader
 from .datalevels import DataLevel
 from ..utils import IndexFinder
 
+from .tableloader import DL2_SUBARRAY_GROUP, DL2_TELESCOPE_GROUP
 
 __all__ = ["HDF5EventSource"]
 
@@ -41,6 +46,14 @@ __all__ = ["HDF5EventSource"]
 logger = logging.getLogger(__name__)
 
 
+DL2_CONTAINERS = {
+    "energy": ReconstructedEnergyContainer,
+    "geometry": ReconstructedGeometryContainer,
+    "classification": ParticleClassificationContainer,
+    "impact": TelescopeImpactParameterContainer,
+}
+
+
 COMPATIBLE_DL1_VERSIONS = [
     "v1.0.0",
     "v1.0.1",
@@ -214,9 +227,6 @@ class HDF5EventSource(EventSource):
         """
         return self._simulation_configs
 
-    def _generator(self):
-        yield from self._generate_events()
-
     def __len__(self):
         n_events = len(self.file_.root.dl1.event.subarray.trigger)
         if self.max_events is not None:
@@ -250,17 +260,10 @@ class HDF5EventSource(EventSource):
 
         return simulation_configs
 
-    def _generate_events(self):
+    def _generator(self):
         """
         Yield ArrayEventContainer to iterate through events.
         """
-        data = ArrayEventContainer()
-        # Maybe take some other metadata, but there are still some 'unknown'
-        # written out by the stage1 tool
-        data.meta["origin"] = self.file_.root._v_attrs["CTA PROCESS TYPE"]
-        data.meta["input_url"] = self.input_url
-        data.meta["max_events"] = self.max_events
-
         self.reader = HDF5TableReader(self.file_)
 
         if DataLevel.R1 in self.datalevels:
@@ -337,6 +340,49 @@ class HDF5EventSource(EventSource):
                     for tel in self.file_.root.dl1.event.telescope.parameters
                 }
 
+        dl2_readers = {}
+        if DL2_SUBARRAY_GROUP in self.file_.root:
+            dl2_group = self.file_.root[DL2_SUBARRAY_GROUP]
+
+            for kind, group in dl2_group._v_children.items():
+
+                try:
+                    container = DL2_CONTAINERS[kind]
+                except KeyError:
+                    self.log.warning("Unknown DL2 stereo group %s", kind)
+                    continue
+
+                dl2_readers[kind] = {
+                    algorithm: HDF5TableReader(self.file_).read(
+                        table._v_pathname,
+                        containers=container,
+                    )
+                    for algorithm, table in group._v_children.items()
+                }
+
+        dl2_tel_readers = {}
+        if DL2_TELESCOPE_GROUP in self.file_.root:
+            dl2_group = self.file_.root[DL2_TELESCOPE_GROUP]
+
+            for kind, group in dl2_group._v_children.items():
+                try:
+                    container = DL2_CONTAINERS[kind]
+                except KeyError:
+                    self.log.warning("Unknown DL2 telescope group %s", kind)
+                    continue
+
+                dl2_tel_readers[kind] = {}
+                for name, algorithm_group in group._v_children.items():
+                    dl2_tel_readers[kind][name] = {
+                        key: HDF5TableReader(self.file_).read(
+                            table._v_pathname,
+                            containers=container,
+                            prefixes=True,
+                        )
+                        for key, table in algorithm_group._v_children.items()
+                    }
+
+        true_impact_readers = {}
         if self.is_simulation:
             # simulated shower wide information
             mc_shower_reader = HDF5TableReader(self.file_).read(
@@ -344,7 +390,15 @@ class HDF5EventSource(EventSource):
                 SimulatedShowerContainer,
                 prefixes="true",
             )
-            data.simulation = SimulatedEventContainer()
+            if "impact" in self.file_.root.simulation.event.telescope:
+                true_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(
@@ -369,15 +423,18 @@ class HDF5EventSource(EventSource):
 
         counter = 0
         for trigger, index in events:
-            data.dl1.tel.clear()
-            if self.is_simulation:
-                data.simulation.tel.clear()
-            data.pointing.tel.clear()
-            data.trigger.tel.clear()
+            data = ArrayEventContainer(
+                trigger=trigger,
+                count=counter,
+                index=index,
+                simulation=SimulatedEventContainer() if self.is_simulation else None,
+            )
+            # Maybe take some other metadata, but there are still some 'unknown'
+            # written out by the stage1 tool
+            data.meta["origin"] = self.file_.root._v_attrs["CTA PROCESS TYPE"]
+            data.meta["input_url"] = self.input_url
+            data.meta["max_events"] = self.max_events
 
-            data.count = counter
-            data.trigger = trigger
-            data.index = index
             data.trigger.tels_with_trigger = self._full_subarray.tel_mask_to_tel_ids(
                 data.trigger.tels_with_trigger
             )
@@ -408,52 +465,53 @@ class HDF5EventSource(EventSource):
             if self.is_simulation:
                 data.simulation.shower = next(mc_shower_reader)
 
-            for tel in data.trigger.tel.keys():
-                key = f"tel_{tel:03d}"
-                if self.allowed_tels and tel not in self.allowed_tels:
+            for tel_id in data.trigger.tel.keys():
+                key = f"tel_{tel_id:03d}"
+                if self.allowed_tels and tel_id not in self.allowed_tels:
                     continue
 
+                if key in true_impact_readers:
+                    data.simulation.tel[tel_id].impact = next(true_impact_readers[key])
+
                 if DataLevel.R1 in self.datalevels:
-                    data.r1.tel[tel] = next(waveform_readers[key])
+                    data.r1.tel[tel_id] = next(waveform_readers[key])
 
                 if self.has_simulated_dl1:
-                    simulated = data.simulation.tel[tel]
+                    simulated = data.simulation.tel[tel_id]
 
                 if DataLevel.DL1_IMAGES in self.datalevels:
                     if key not in image_readers:
                         logger.debug(
-                            f"Triggered telescope {tel} is missing "
+                            f"Triggered telescope {tel_id} is missing "
                             "from the image table."
                         )
                         continue
 
-                    data.dl1.tel[tel] = next(image_readers[key])
+                    data.dl1.tel[tel_id] = next(image_readers[key])
 
                     if self.has_simulated_dl1:
                         if key not in simulated_image_iterators:
                             logger.warning(
-                                f"Triggered telescope {tel} is missing "
+                                f"Triggered telescope {tel_id} is missing "
                                 "from the simulated image table, but was present at the "
                                 "reconstructed image table."
                             )
                             continue
-                        simulated_image_row = next(
-                            simulated_image_iterators[f"tel_{tel:03d}"]
-                        )
+                        simulated_image_row = next(simulated_image_iterators[key])
                         simulated.true_image = simulated_image_row["true_image"]
 
                 if DataLevel.DL1_PARAMETERS in self.datalevels:
-                    if f"tel_{tel:03d}" not in param_readers.keys():
+                    if f"tel_{tel_id:03d}" not in param_readers.keys():
                         logger.debug(
-                            f"Triggered telescope {tel} is missing "
+                            f"Triggered telescope {tel_id} is missing "
                             "from the parameters table."
                         )
                         continue
                     # Is there a smarter way to unpack this?
                     # Best would probbaly be if we could directly read
                     # into the ImageParametersContainer
-                    params = next(param_readers[f"tel_{tel:03d}"])
-                    data.dl1.tel[tel].parameters = ImageParametersContainer(
+                    params = next(param_readers[key])
+                    data.dl1.tel[tel_id].parameters = ImageParametersContainer(
                         hillas=params[0],
                         timing=params[1],
                         leakage=params[2],
@@ -463,15 +521,15 @@ class HDF5EventSource(EventSource):
                         peak_time_statistics=params[6],
                     )
                     if self.has_simulated_dl1:
-                        if f"tel_{tel:03d}" not in param_readers.keys():
+                        if f"tel_{tel_id:03d}" not in param_readers.keys():
                             logger.debug(
-                                f"Triggered telescope {tel} is missing "
+                                f"Triggered telescope {tel_id} is missing "
                                 "from the simulated parameters table, but was "
                                 "present at the reconstructed parameters table."
                             )
                             continue
                         simulated_params = next(
-                            simulated_param_readers[f"tel_{tel:03d}"]
+                            simulated_param_readers[f"tel_{tel_id:03d}"]
                         )
                         simulated.true_parameters = ImageParametersContainer(
                             hillas=simulated_params[0],
@@ -481,6 +539,16 @@ class HDF5EventSource(EventSource):
                             intensity_statistics=simulated_params[4],
                         )
 
+                for kind, algorithms in dl2_tel_readers.items():
+                    c = getattr(data.dl2.tel[tel_id], kind)
+                    for algorithm, readers in algorithms.items():
+                        c[algorithm] = next(readers[key])
+
+            for kind, readers in dl2_readers.items():
+                c = getattr(data.dl2.stereo, kind)
+                for algorithm, reader in readers.items():
+                    c[algorithm] = next(reader)
+
             yield data
             counter += 1
 
diff --git a/ctapipe/io/hdf5tableio.py b/ctapipe/io/hdf5tableio.py
index 8af44714..4c6bb637 100644
--- a/ctapipe/io/hdf5tableio.py
+++ b/ctapipe/io/hdf5tableio.py
@@ -17,7 +17,7 @@ from .tableio import (
     EnumColumnTransform,
     QuantityColumnTransform,
 )
-from ..core import Container
+from ..core import Container, Map
 
 __all__ = ["HDF5TableWriter", "HDF5TableReader"]
 
@@ -141,12 +141,20 @@ class HDF5TableWriter(TableWriter):
         typename = ""
         shape = 1
 
+        if isinstance(value, Container):
+            self.log.debug("Ignoring sub-container: %s/%s", table_name, name)
+            return
+
+        if isinstance(value, Map):
+            self.log.debug("Ignoring map-field: %s/%s", table_name, name)
+            return
+
         if self._is_column_excluded(table_name, name):
-            self.log.debug(f"excluded column: {table_name}/{name}")
+            self.log.debug("excluded column: %s/%s", table_name, name)
             return
 
         if name in schema.columns:
-            self.log.warning(f"Found duplicated column {name}, skipping")
+            self.log.warning("Found duplicated column %s, skipping", name)
             return
 
         # apply any user-defined transforms first
diff --git a/ctapipe/io/metadata.py b/ctapipe/io/metadata.py
index e8b08520..da77c6cd 100644
--- a/ctapipe/io/metadata.py
+++ b/ctapipe/io/metadata.py
@@ -55,7 +55,7 @@ CONVERSIONS = {Time: lambda t: t.utc.iso, list: str}
 
 
 class Contact(Configurable):
-    """ Contact information """
+    """Contact information"""
 
     name = Unicode("unknown").tag(config=True)
     email = Unicode("unknown").tag(config=True)
@@ -63,7 +63,7 @@ class Contact(Configurable):
 
     @default("name")
     def default_name(self):
-        """ if no name specified, use the system's user name"""
+        """if no name specified, use the system's user name"""
         try:
             return pwd.getpwuid(os.getuid()).pw_gecos
         except RuntimeError:
@@ -90,17 +90,17 @@ class Product(HasTraits):
     # pylint: disable=no-self-use
     @default("creation_time")
     def default_time(self):
-        """ return current time by default """
+        """return current time by default"""
         return Time.now().iso
 
     @default("id_")
     def default_product_id(self):
-        """ default id is a UUID """
+        """default id is a UUID"""
         return str(uuid.uuid4())
 
 
 class Process(HasTraits):
-    """ Process (top-level workflow) information """
+    """Process (top-level workflow) information"""
 
     type_ = Enum(["Observation", "Simulation", "Other"], "Other")
     subtype = Unicode("")
@@ -108,11 +108,11 @@ class Process(HasTraits):
 
 
 class Activity(HasTraits):
-    """ Activity (tool) information """
+    """Activity (tool) information"""
 
     @classmethod
     def from_provenance(cls, activity):
-        """ construct Activity metadata from existing ActivityProvenance object"""
+        """construct Activity metadata from existing ActivityProvenance object"""
         return Activity(
             name=activity["activity_name"],
             type_="software",
@@ -132,12 +132,12 @@ class Activity(HasTraits):
     # pylint: disable=no-self-use
     @default("start_time")
     def default_time(self):
-        """ default time is now """
+        """default time is now"""
         return Time.now().iso
 
 
-class Instrument(HasTraits):
-    """ Instrumental Context """
+class Instrument(Configurable):
+    """Instrumental Context"""
 
     site = Enum(
         [
@@ -155,7 +155,7 @@ class Instrument(HasTraits):
         "Other",
         help="Which site of CTA (or external telescope) "
         "this instrument is associated with",
-    )
+    ).tag(config=True)
     class_ = Enum(
         [
             "Array",
@@ -170,15 +170,21 @@ class Instrument(HasTraits):
             "Other",
         ],
         "Other",
-    )
-    type_ = Unicode("unspecified")
-    subtype = Unicode("unspecified")
-    version = Unicode("unspecified")
-    id_ = Unicode("unspecified")
+    ).tag(config=True)
+    type_ = Unicode("unspecified").tag(config=True)
+    subtype = Unicode("unspecified").tag(config=True)
+    version = Unicode("unspecified").tag(config=True)
+    id_ = Unicode("unspecified").tag(config=True)
+
+    def __repr__(self):
+        return (
+            f"Contact({self.site=}, {self.class_=}, {self.type_=}, "
+            f"{self.subtype=}, {self.version=}, {self.id_=})"
+        )
 
 
 def _to_dict(hastraits_instance, prefix=""):
-    """ helper to convert a HasTraits to a dict with keys
+    """helper to convert a HasTraits to a dict with keys
     in the required CTA format (upper-case, space separated)
     """
     res = {}
@@ -199,8 +205,8 @@ def _to_dict(hastraits_instance, prefix=""):
 
 
 class Reference(HasTraits):
-    """ All the reference Metadata required for a CTA output file, plus a way to turn
-    it into a dict() for easy addition to the header of a file """
+    """All the reference Metadata required for a CTA output file, plus a way to turn
+    it into a dict() for easy addition to the header of a file"""
 
     contact = Instance(Contact)
     product = Instance(Product)
@@ -225,7 +231,7 @@ class Reference(HasTraits):
         return meta
 
 
-def write_to_hdf5(metadata, h5file, path='/'):
+def write_to_hdf5(metadata, h5file, path="/"):
     """
     Write metadata fields to a PyTables HDF5 file handle.
 
@@ -246,7 +252,7 @@ def write_to_hdf5(metadata, h5file, path='/'):
             node._v_attrs[key] = value  # pylint: disable=protected-access
 
 
-def read_metadata(h5file, path='/'):
+def read_metadata(h5file, path="/"):
     """
     Read metadata from an hdf5 file
 
diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py
index 8eb52425..c6e0d3e7 100644
--- a/ctapipe/io/simteleventsource.py
+++ b/ctapipe/io/simteleventsource.py
@@ -1,3 +1,4 @@
+from collections import defaultdict
 import warnings
 from gzip import GzipFile
 from io import BufferedReader
@@ -13,21 +14,27 @@ from eventio.simtel.simtelfile import SimTelFile
 from ..calib.camera.gainselection import GainSelector
 from ..containers import (
     ArrayEventContainer,
+    EventIndexContainer,
     EventType,
-    SimulatedEventContainer,
-    SimulationConfigContainer,
+    PointingContainer,
+    R0CameraContainer,
+    R1CameraContainer,
     SimulatedCameraContainer,
+    SimulatedEventContainer,
     SimulatedShowerContainer,
+    SimulationConfigContainer,
+    TelescopeImpactParameterContainer,
     TelescopePointingContainer,
     TelescopeTriggerContainer,
+    TriggerContainer,
 )
 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 +45,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
 
@@ -393,27 +401,10 @@ class SimTelEventSource(EventSource):
             warnings.warn(msg)
 
     def _generate_events(self):
-        data = ArrayEventContainer()
-        data.simulation = SimulatedEventContainer()
-        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)
-
         for counter, array_event in enumerate(self.file_):
-            # this should be done in a nicer way to not re-allocate the
-            # data each time (right now it's just deleted and garbage
-            # collected)
-            data.r0.tel.clear()
-            data.r1.tel.clear()
-            data.dl0.tel.clear()
-            data.dl1.tel.clear()
-            data.pointing.tel.clear()
-            data.simulation.tel.clear()
-            data.trigger.tel.clear()
 
             event_id = array_event.get("event_id", 0)
             if event_id == 0:
@@ -421,13 +412,23 @@ class SimTelEventSource(EventSource):
                 event_id = pseudo_event_id
 
             obs_id = self.file_.header["run"]
-            data.count = counter
-            data.index.obs_id = obs_id
-            data.index.event_id = event_id
 
-            self._fill_trigger_info(data, array_event)
-            if data.trigger.event_type == EventType.SUBARRAY:
-                self._fill_simulated_event_information(data, array_event)
+            trigger = self._fill_trigger_info(array_event)
+            if trigger.event_type == EventType.SUBARRAY:
+                shower = self._fill_simulated_event_information(array_event)
+            else:
+                shower = None
+
+            data = ArrayEventContainer(
+                simulation=SimulatedEventContainer(shower=shower),
+                pointing=self._fill_array_pointing(),
+                index=EventIndexContainer(obs_id=obs_id, event_id=event_id),
+                count=counter,
+                trigger=trigger,
+            )
+            data.meta["origin"] = "hessio"
+            data.meta["input_url"] = self.input_url
+            data.meta["max_events"] = self.max_events
 
             telescope_events = array_event["telescope_events"]
             tracking_positions = array_event["tracking_positions"]
@@ -440,6 +441,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,20 +462,31 @@ 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(
                     tracking_positions[tel_id]
                 )
 
-                r0 = data.r0.tel[tel_id]
-                r1 = data.r1.tel[tel_id]
-                r0.waveform = adc_samples
+                data.r0.tel[tel_id] = R0CameraContainer(waveform=adc_samples)
 
                 cam_mon = array_event["camera_monitorings"][tel_id]
                 pedestal = cam_mon["pedestal"] / cam_mon["n_ped_slices"]
@@ -478,7 +498,7 @@ class SimTelEventSource(EventSource):
                 mon.calibration.dc_to_pe = dc_to_pe
                 mon.calibration.pedestal_per_sample = pedestal
 
-                r1.waveform, r1.selected_gain_channel = apply_simtel_r1_calibration(
+                r1_waveform, selected_gain_channel = apply_simtel_r1_calibration(
                     adc_samples,
                     pedestal,
                     dc_to_pe,
@@ -486,13 +506,17 @@ class SimTelEventSource(EventSource):
                     self.calib_scale,
                     self.calib_shift,
                 )
+                data.r1.tel[tel_id] = R1CameraContainer(
+                    waveform=r1_waveform,
+                    selected_gain_channel=selected_gain_channel,
+                )
 
                 # get time_shift from laser calibration
                 time_calib = array_event["laser_calibrations"][tel_id]["tm_calib"]
                 pix_index = np.arange(n_pixels)
 
                 dl1_calib = data.calibration.tel[tel_id].dl1
-                dl1_calib.time_shift = time_calib[r1.selected_gain_channel, pix_index]
+                dl1_calib.time_shift = time_calib[selected_gain_channel, pix_index]
 
             yield data
 
@@ -517,36 +541,39 @@ class SimTelEventSource(EventSource):
 
         return TelescopePointingContainer(azimuth=azimuth, altitude=altitude)
 
-    def _fill_trigger_info(self, data, array_event):
+    def _fill_trigger_info(self, array_event):
         trigger = array_event["trigger_information"]
 
         if array_event["type"] == "data":
-            data.trigger.event_type = EventType.SUBARRAY
+            event_type = EventType.SUBARRAY
 
         elif array_event["type"] == "calibration":
             # if using eventio >= 1.1.1, we can use the calibration_type
-            data.trigger.event_type = SIMTEL_TO_CTA_EVENT_TYPE.get(
+            event_type = SIMTEL_TO_CTA_EVENT_TYPE.get(
                 array_event.get("calibration_type", -1), EventType.OTHER_CALIBRATION
             )
 
         else:
-            data.trigger.event_type = EventType.UNKNOWN
+            event_type = EventType.UNKNOWN
 
-        data.trigger.tels_with_trigger = trigger["triggered_telescopes"]
         if self.allowed_tels:
-            data.trigger.tels_with_trigger = np.intersect1d(
-                data.trigger.tels_with_trigger,
+            tels_with_trigger = np.intersect1d(
+                trigger["triggered_telescopes"],
                 self.subarray.tel_ids,
                 assume_unique=True,
             )
+        else:
+            tels_with_trigger = trigger["triggered_telescopes"]
+
         central_time = parse_simtel_time(trigger["gps_time"])
-        data.trigger.time = central_time
 
+        tel = defaultdict(TelescopeTriggerContainer)
         for tel_id, time in zip(
             trigger["triggered_telescopes"], trigger["trigger_times"]
         ):
             if self.allowed_tels and tel_id not in self.allowed_tels:
                 continue
+
             # telescope time is relative to central trigger in ns
             time = Time(
                 central_time.jd1,
@@ -567,21 +594,31 @@ class SimTelEventSource(EventSource):
                     n_trigger_pixels = pixel_list["pixels"]
                     trigger_pixels = pixel_list["pixel_list"]
 
-            trigger = data.trigger.tel[tel_id] = TelescopeTriggerContainer(
+            tel[tel_id] = TelescopeTriggerContainer(
                 time=time,
                 n_trigger_pixels=n_trigger_pixels,
                 trigger_pixels=trigger_pixels,
             )
+        return TriggerContainer(
+            event_type=event_type,
+            time=central_time,
+            tels_with_trigger=tels_with_trigger,
+            tel=tel,
+        )
 
-    def _fill_array_pointing(self, data):
+    def _fill_array_pointing(self):
         if self.file_.header["tracking_mode"] == 0:
             az, alt = self.file_.header["direction"]
-            data.pointing.array_altitude = u.Quantity(alt, u.rad)
-            data.pointing.array_azimuth = u.Quantity(az, u.rad)
+            return PointingContainer(
+                array_altitude=u.Quantity(alt, u.rad),
+                array_azimuth=u.Quantity(az, u.rad),
+            )
         else:
             ra, dec = self.file_.header["direction"]
-            data.pointing.array_ra = u.Quantity(ra, u.rad)
-            data.pointing.array_dec = u.Quantity(dec, u.rad)
+            return PointingContainer(
+                array_ra=u.Quantity(ra, u.rad),
+                array_dec=u.Quantity(dec, u.rad),
+            )
 
     def _parse_simulation_header(self):
         """
@@ -636,14 +673,13 @@ class SimTelEventSource(EventSource):
         return {obs_id: simulation_config}
 
     @staticmethod
-    def _fill_simulated_event_information(data, array_event):
+    def _fill_simulated_event_information(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(
+        return SimulatedShowerContainer(
             energy=u.Quantity(mc_shower["energy"], u.TeV),
             alt=Angle(mc_shower["altitude"], u.rad),
             az=Angle(mc_shower["azimuth"], u.rad),
diff --git a/ctapipe/io/tableloader.py b/ctapipe/io/tableloader.py
index 9afea2b2..ba83e5ae 100644
--- a/ctapipe/io/tableloader.py
+++ b/ctapipe/io/tableloader.py
@@ -7,11 +7,11 @@ from pathlib import Path
 import re
 from collections import defaultdict
 from typing import Dict
+from astropy.utils.decorators import lazyproperty
 
 import numpy as np
-from astropy.table import join, vstack, Table
+from astropy.table import vstack, Table
 import tables
-from traitlets.traitlets import Instance
 
 from ..core import Component, traits, Provenance
 from ..instrument import SubarrayDescription
@@ -25,6 +25,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"
@@ -36,12 +37,64 @@ SUBARRAY_EVENT_KEYS = ["obs_id", "event_id"]
 TELESCOPE_EVENT_KEYS = ["obs_id", "event_id", "tel_id"]
 
 
+class ChunkIterator:
+    """An iterator that calls a function on advancemnt
+
+    Parameters
+    ----------
+    func : Callable
+        Signature must be ``function(*args, start=start, stop=stop, **kwargs)``
+        Where start and stop are the first and last (noninclusive) indicies of
+        a chunk.
+    n_total : int
+        Total number of elements
+    chunk_size : size of one chunk, last chunk will have size <= chunk_size
+
+    *args and **kwargs will be passed to `func`
+
+    """
+
+    def __init__(
+        self,
+        func,
+        n_total,
+        chunk_size,
+        *args,
+        **kwargs,
+    ):
+        self.func = func
+        self.n_total = n_total
+        self.chunk_size = chunk_size
+        self._current_chunk = 0
+        self.n_chunks = int(np.ceil(self.n_total / self.chunk_size))
+        self.args = args
+        self.kwargs = kwargs
+
+    def __len__(self):
+        return self.n_chunks
+
+    def __iter__(self):
+        self._current_chunk = 0
+        return self
+
+    def __next__(self):
+        if self._current_chunk == self.n_chunks:
+            raise StopIteration
+
+        chunk = self._current_chunk
+        start = chunk * self.chunk_size
+        stop = min(self.n_total, (chunk + 1) * self.chunk_size)
+
+        self._current_chunk += 1
+        return self.func(*self.args, start=start, stop=stop, **self.kwargs)
+
+
 def _empty_telescope_events_table():
     """
     Create a new astropy table with correct column names and dtypes
     for telescope event based data.
     """
-    return Table(names=TELESCOPE_EVENT_KEYS, dtype=[np.int64, np.int64, np.int16])
+    return Table(names=TELESCOPE_EVENT_KEYS, dtype=[np.int32, np.int64, np.int16])
 
 
 def _empty_subarray_events_table():
@@ -49,7 +102,7 @@ def _empty_subarray_events_table():
     Create a new astropy table with correct column names and dtypes
     for subarray event based data.
     """
-    return Table(names=SUBARRAY_EVENT_KEYS, dtype=[np.int64, np.int64])
+    return Table(names=SUBARRAY_EVENT_KEYS, dtype=[np.int32, np.int64])
 
 
 def _get_structure(h5file):
@@ -116,7 +169,6 @@ class TableLoader(Component):
     input_url = traits.Path(
         directory_ok=False, exists=True, allow_none=True, default_value=None
     ).tag(config=True)
-    h5file = Instance(tables.File, default_value=None, allow_none=True).tag(config=True)
 
     load_dl1_images = traits.Bool(False, help="load extracted images").tag(config=True)
     load_dl1_parameters = traits.Bool(
@@ -143,24 +195,25 @@ class TableLoader(Component):
         False, help="join subarray instrument information to each event"
     ).tag(config=True)
 
-    def __init__(self, input_url=None, **kwargs):
+    def __init__(self, input_url=None, h5file=None, **kwargs):
+        self._should_close = False
+
         # enable using input_url as posarg
-        if input_url not in {None, traits.Undefined} and not isinstance(
-            input_url, tables.File
-        ):
+        if input_url not in {None, traits.Undefined}:
             kwargs["input_url"] = input_url
 
         super().__init__(**kwargs)
-
-        if self.h5file is None and self.input_url is None:
+        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 self.h5file is None:
+        if h5file is None:
             self.h5file = tables.open_file(self.input_url, mode="r")
             self._should_close = True
         else:
-            self.input_url = Path(self.h5file.filename)
+            if not isinstance(h5file, tables.File):
+                raise TypeError("h5file must be a tables.File")
+            self.input_url = Path(h5file.filename)
+            self.h5file = h5file
 
         self.subarray = SubarrayDescription.from_hdf(self.h5file)
 
@@ -173,20 +226,7 @@ class TableLoader(Component):
 
         self.instrument_table = None
         if self.load_instrument:
-            table = self.subarray.to_table()
-            optics = self.subarray.to_table(kind="optics")
-            optics["optics_index"] = np.arange(len(optics))
-            optics.remove_columns(["name", "description", "type"])
-            table = join(
-                table,
-                optics,
-                keys="optics_index",
-                # conflicts for TAB_VER, TAB_TYPE, not needed here, ignore
-                metadata_conflicts="silent",
-            )
-
-            table.remove_columns(["optics_index", "camera_index"])
-            self.instrument_table = table
+            self.instrument_table = self.subarray.to_table('joined')
 
     def close(self):
         """Close the underlying hdf5 file"""
@@ -202,7 +242,7 @@ class TableLoader(Component):
     def __exit__(self, exc_type, exc_value, traceback):
         self.close()
 
-    def _read_telescope_table(self, group, tel_id):
+    def _read_telescope_table(self, group, tel_id, start=None, stop=None):
         if self.structure == "by_id":
             key = f"{group}/tel_{tel_id:03d}"
             condition = None
@@ -211,13 +251,15 @@ class TableLoader(Component):
             condition = f"tel_id == {tel_id}"
 
         if key in self.h5file:
-            table = read_table(self.h5file, key, condition=condition)
+            table = read_table(
+                self.h5file, key, condition=condition, start=start, stop=stop
+            )
         else:
             table = _empty_telescope_events_table()
 
         return table
 
-    def read_subarray_events(self):
+    def read_subarray_events(self, start=None, stop=None):
         """Read subarray-based event information.
 
         Returns
@@ -228,11 +270,11 @@ class TableLoader(Component):
         table = _empty_subarray_events_table()
 
         if self.load_trigger:
-            trigger = read_table(self.h5file, TRIGGER_TABLE)
+            trigger = read_table(self.h5file, TRIGGER_TABLE, start=start, stop=stop)
             table = _join_subarray_events(table, trigger)
 
         if self.load_simulated and SHOWER_TABLE in self.h5file:
-            showers = read_table(self.h5file, SHOWER_TABLE)
+            showers = read_table(self.h5file, SHOWER_TABLE, start=start, stop=stop)
             table = _join_subarray_events(table, showers)
 
         if self.load_dl2:
@@ -242,7 +284,12 @@ class TableLoader(Component):
                     group = self.h5file.root[group_path]
 
                     for algorithm in group._v_children:
-                        dl2 = read_table(self.h5file, f"{group_path}/{algorithm}")
+                        dl2 = read_table(
+                            self.h5file,
+                            f"{group_path}/{algorithm}",
+                            start=start,
+                            stop=stop,
+                        )
 
                         # add the algorithm as prefix to distinguish multiple
                         # algorithms predicting the same quantities
@@ -252,7 +299,24 @@ class TableLoader(Component):
                         table = _join_subarray_events(table, dl2)
         return table
 
-    def _read_telescope_events_for_id(self, tel_id):
+    def read_subarray_events_chunked(self, chunk_size):
+        """
+        Iterate over chunks of subarray events.
+
+        Parameters
+        ----------
+        chunk_size: int
+            Number of subarray events to load per chunk
+        """
+        n_events = self.h5file.root[TRIGGER_TABLE].shape[0]
+
+        return ChunkIterator(
+            self.read_subarray_events,
+            n_total=n_events,
+            chunk_size=chunk_size,
+        )
+
+    def _read_telescope_events_for_id(self, tel_id, start=None, stop=None):
         """Read telescope-based event information for a single telescope.
 
         This is the most low-level function doing the actual reading.
@@ -261,35 +325,45 @@ class TableLoader(Component):
         ----------
         tel_id: int
             Telescope identification number.
+        start: int
+            First row to read
+        stop: int
+            Last row to read (non-inclusive)
 
         Returns
         -------
         table: astropy.io.Table
             Table with primary index columns "obs_id", "event_id" and "tel_id".
         """
-
         if tel_id is None:
             raise ValueError("Please, specify a telescope ID.")
 
         table = _empty_telescope_events_table()
 
         if self.load_dl1_parameters:
-            parameters = self._read_telescope_table(PARAMETERS_GROUP, tel_id)
+            parameters = self._read_telescope_table(
+                PARAMETERS_GROUP, tel_id, start=start, stop=stop
+            )
             table = _join_telescope_events(table, parameters)
 
         if self.load_dl1_images:
-            images = self._read_telescope_table(IMAGES_GROUP, tel_id)
+            images = self._read_telescope_table(
+                IMAGES_GROUP, tel_id, start=start, stop=stop
+            )
             table = _join_telescope_events(table, images)
 
         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]
 
                     for algorithm in group._v_children:
                         path = f"{group_path}/{algorithm}"
-                        dl2 = self._read_telescope_table(path, tel_id)
+                        dl2 = self._read_telescope_table(
+                            path, tel_id, start=start, stop=stop
+                        )
 
                         # add the algorithm as prefix to distinguish multiple
                         # algorithms predicting the same quantities
@@ -299,15 +373,16 @@ class TableLoader(Component):
                         table = _join_telescope_events(table, dl2)
 
         if self.load_true_images:
-            true_images = self._read_telescope_table(TRUE_IMAGES_GROUP, tel_id)
+            true_images = self._read_telescope_table(
+                TRUE_IMAGES_GROUP, tel_id, start=start, stop=stop
+            )
             table = _join_telescope_events(table, true_images)
 
         if self.load_true_parameters:
-            true_parameters = self._read_telescope_table(TRUE_PARAMETERS_GROUP, tel_id)
-
-            for col in set(true_parameters.colnames) - set(TELESCOPE_EVENT_KEYS):
-                true_parameters.rename_column(col, f"true_{col}")
-
+            true_parameters = self._read_telescope_table(
+                TRUE_PARAMETERS_GROUP, tel_id, start=start, stop=stop
+            )
+            _add_column_prefix(true_parameters, "true", ignore=TELESCOPE_EVENT_KEYS)
             table = _join_telescope_events(table, true_parameters)
 
         if self.load_instrument:
@@ -315,35 +390,57 @@ 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):
+    def _read_telescope_events_for_ids(self, tel_ids, tel_start=None, tel_stop=None):
+        if tel_start is None:
+            tel_start = [None] * len(tel_ids)
+
+        if tel_stop is None:
+            tel_stop = [None] * len(tel_ids)
 
         table = vstack(
-            [self._read_telescope_events_for_id(tel_id) for tel_id in tel_ids]
+            [
+                self._read_telescope_events_for_id(tel_id, start=start, stop=stop)
+                for tel_id, start, stop in zip(tel_ids, tel_start, tel_stop)
+            ]
         )
 
         return table
 
-    def _join_subarray_info(self, table):
-        subarray_events = self.read_subarray_events()
+    def _join_subarray_info(self, table, start=None, stop=None):
+        subarray_events = self.read_subarray_events(start=start, stop=stop)
         table = join_allow_empty(
             table, subarray_events, keys=SUBARRAY_EVENT_KEYS, join_type="left"
         )
         return table
 
-    def read_telescope_events(self, telescopes=None):
+    def read_telescope_events(self, telescopes=None, start=None, stop=None):
         """
         Read telescope-based event information.
 
         If the corresponding traitlets are True, also subarray event information
         is joined onto the table.
 
+        The start, stop parameters enable to only load parts of the file,
+        note however, to maintain integrity of subarray events, these
+        refer to the subarray indices in the file. E.g. ``start=0``, ``stop=10``
+        would load the telescope events corresponding to the first 10 subarray
+        events in the input file.
+
         Parameters
         ----------
         telescopes: Optional[List[Union[int, str, TelescopeDescription]]]
             A list containing any combination of telescope IDs and/or
             telescope descriptions. If None, all available telescopes are read.
+        start: int
+            First *subarray* event to read
+        stop: int
+            Last *subarray* event (non-inclusive)
 
         Returns
         -------
@@ -352,17 +449,74 @@ class TableLoader(Component):
         """
 
         if telescopes is None:
-            tel_ids = self.subarray.tel.keys()
+            tel_ids = tuple(self.subarray.tel.keys())
         else:
             tel_ids = self.subarray.get_tel_ids(telescopes)
 
-        table = self._read_telescope_events_for_ids(tel_ids)
+        tel_start = None
+        tel_stop = None
+        if start is not None or stop is not None:
+
+            indices = self.subarray.tel_ids_to_indices(tel_ids)
+
+            # find first/last row for each telescope
+            if start is not None:
+                tel_start = self._n_telescope_events[start][indices]
+
+            if stop is not None:
+                if stop >= len(self._n_telescope_events):
+                    tel_stop = None
+                else:
+                    tel_stop = self._n_telescope_events[stop][indices]
+
+        table = self._read_telescope_events_for_ids(tel_ids, tel_start, tel_stop)
 
         if any([self.load_trigger, self.load_simulated, self.load_dl2]):
-            table = self._join_subarray_info(table)
+            table = self._join_subarray_info(table, start=start, stop=stop)
 
         return table
 
+    def read_telescope_events_chunked(self, chunk_size, **kwargs):
+        """
+        Iterate over chunks of telescope events.
+
+        Parameters
+        ----------
+        chunk_size: int
+            Number of subarray events to load per chunk.
+            The telescope tables might be larger or smaller than chunk_size
+            depending on the selected telescopes.
+
+        **kwargs are passed to `read_telescope_events`
+        """
+        n_events = self.h5file.root[TRIGGER_TABLE].shape[0]
+
+        return ChunkIterator(
+            self.read_telescope_events,
+            n_total=n_events,
+            chunk_size=chunk_size,
+            **kwargs,
+        )
+
+    @lazyproperty
+    def _n_telescope_events(self):
+        """
+        Number of telescope events in the file for each telescope previous
+        to the nth subarray event.
+        """
+        # we need to load the trigger table until "stop" to
+        # know which telescopes participated in which events
+        table = self.h5file.root[TRIGGER_TABLE]
+        n_events = table.shape[0]
+        n_telescopes = table.coldescrs["tels_with_trigger"].shape[0]
+
+        tels_with_trigger = np.zeros((n_events + 1, n_telescopes), dtype=np.bool_)
+        table.read(field="tels_with_trigger", out=tels_with_trigger[1:])
+        tels_with_trigger = tels_with_trigger.astype(np.uint32)
+
+        np.cumsum(tels_with_trigger, out=tels_with_trigger, axis=0)
+        return tels_with_trigger
+
     def read_telescope_events_by_type(self, telescopes=None) -> Dict[str, Table]:
         """Read telescope-based event information.
 
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..8ed30b6c 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_id in event.trigger.tels_with_trigger:
+            tel_index = self.subarray.tel_indices[tel_id]
+            event.dl2.tel[tel_id].impact[k] = TelescopeImpactParameterContainer(
+                distance=impact_distances[tel_index]
+            )
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/merge.py b/ctapipe/tools/merge.py
index d458fa32..e965d14d 100644
--- a/ctapipe/tools/merge.py
+++ b/ctapipe/tools/merge.py
@@ -11,6 +11,8 @@ import tables
 import numpy as np
 from tqdm.auto import tqdm
 
+from ctapipe.utils.arrays import recarray_drop_columns
+
 from ..io import metadata as meta, HDF5EventSource, get_hdf5_datalevels
 from ..io import HDF5TableWriter
 from ..core import Provenance, Tool, traits
@@ -58,13 +60,16 @@ nodes_with_tels = {
     "/simulation/event/telescope/parameters",
     "/simulation/event/telescope/images",
 }
-image_nodes = {"/simulation/event/telescope/images", "/dl1/event/telescope/images"}
+image_nodes = {
+    "/dl1/event/telescope/images",
+}
 parameter_nodes = {
     "/simulation/event/telescope/parameters",
     "/dl1/event/telescope/parameters",
 }
 
-simulation_images = {"/simulation/event/telescope/images"}
+SIMULATED_IMAGE_GROUP = "/simulation/event/telescope/images"
+simulation_images = {SIMULATED_IMAGE_GROUP}
 
 dl2_subarray_nodes = {"/dl2/event/subarray/geometry"}
 
@@ -235,13 +240,13 @@ class MergeTool(Tool):
         # setup required nodes
         self.usable_nodes = all_nodes.copy()
 
-        if self.skip_simu_images is True:
+        if self.skip_simu_images:
             self.usable_nodes -= simulation_images
 
-        if self.skip_images is True:
+        if self.skip_images:
             self.usable_nodes -= image_nodes
 
-        if self.skip_parameters is True:
+        if self.skip_parameters:
             self.usable_nodes -= parameter_nodes
 
         # Use first file as reference to setup nodes to merge
@@ -344,7 +349,7 @@ class MergeTool(Tool):
                 if node in file:
                     file.copy_node(node, newparent=target_group)
 
-    def _merge_tel_group(self, file, input_node):
+    def _merge_tel_group(self, file, input_node, filter_columns=None):
         """Add a group that has one child table per telescope (type) to outputfile"""
         if not isinstance(input_node, tables.Group):
             raise TypeError(f"node must be a `tables.Group`, got {input_node}")
@@ -356,23 +361,35 @@ class MergeTool(Tool):
 
         for tel_name, table in input_node._v_children.items():
             if not self.allowed_tels or tel_name in self.allowed_tel_names:
-                self._merge_table(file, table)
+                self._merge_table(file, table, filter_columns=filter_columns)
 
-    def _merge_table(self, file, input_node):
+    def _merge_table(self, file, input_node, filter_columns=None):
         if not isinstance(input_node, tables.Table):
             raise TypeError(f"node must be a `tables.Table`, got {input_node}")
 
         node_path = input_node._v_pathname
         if node_path in self.output_file:
             output_table = self.output_file.get_node(node_path)
-            output_table.append(input_node[:].astype(output_table.dtype))
+            input_table = input_node[:]
+            if filter_columns is not None:
+                input_table = recarray_drop_columns(input_table, filter_columns)
+            output_table.append(input_table.astype(output_table.dtype))
         else:
             group_path, _ = os.path.split(node_path)
             if group_path not in self.output_file:
                 self._create_group(group_path)
 
             target_group = self.output_file.root[group_path]
-            file.copy_node(node_path, newparent=target_group)
+            if filter_columns is None:
+                file.copy_node(node_path, newparent=target_group)
+            else:
+                input_table = recarray_drop_columns(input_node[:], filter_columns)
+                where, name = os.path.split(node_path)
+                self.output_file.create_table(
+                    where, name, filters=input_node.filters,
+                    createparents=True,
+                    obj=input_table,
+                )
 
     def merge_tables(self, file):
         """Go over all mergeable nodes and append to outputfile"""
@@ -387,7 +404,11 @@ class MergeTool(Tool):
 
                 # nodes with child tables per telescope (type)
                 if node_path in nodes_with_tels:
-                    self._merge_tel_group(file, node)
+                    if node_path == SIMULATED_IMAGE_GROUP and self.skip_images:
+                        filter_columns = ["true_image"]
+                    else:
+                        filter_columns = None
+                    self._merge_tel_group(file, node, filter_columns=filter_columns)
 
                 # groups of tables (e.g. dl2)
                 elif isinstance(node, tables.Group):
diff --git a/ctapipe/tools/process.py b/ctapipe/tools/process.py
index b988c196..290ae736 100644
--- a/ctapipe/tools/process.py
+++ b/ctapipe/tools/process.py
@@ -14,6 +14,7 @@ from ..io import DataLevel, DataWriter, EventSource, SimTelEventSource, write_ta
 from ..io.datawriter import DATA_MODEL_VERSION
 from ..reco import ShowerProcessor
 from ..utils import EventTypeFilter
+from ..io import metadata
 
 
 COMPATIBLE_DATALEVELS = [
@@ -111,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",
@@ -134,7 +129,14 @@ class ProcessorTool(Tool):
     }
 
     classes = (
-        [CameraCalibrator, DataWriter, ImageProcessor, ShowerProcessor]
+        [
+            CameraCalibrator,
+            DataWriter,
+            ImageProcessor,
+            ShowerProcessor,
+            metadata.Instrument,
+            metadata.Contact,
+        ]
         + classes_with_traits(EventSource)
         + classes_with_traits(ImageCleaner)
         + classes_with_traits(ImageExtractor)
@@ -184,10 +186,10 @@ class ProcessorTool(Tool):
 
     @property
     def should_compute_dl2(self):
-        """ returns true if we should compute DL2 info """
+        """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):
@@ -281,7 +283,7 @@ class ProcessorTool(Tool):
 
 
 def main():
-    """ run the tool"""
+    """run the tool"""
     tool = ProcessorTool()
     tool.run()
 
diff --git a/ctapipe/utils/arrays.py b/ctapipe/utils/arrays.py
new file mode 100644
index 00000000..2f6af42c
--- /dev/null
+++ b/ctapipe/utils/arrays.py
@@ -0,0 +1,29 @@
+'''
+Some utility functions to work with numpy (rec) arrays
+'''
+import numpy as np
+from numpy.lib.recfunctions import repack_fields
+
+
+def recarray_drop_columns(array, columns):
+    '''
+    Remove columns from rec array
+    '''
+    to_use = [col for col in array.dtype.names if col not in columns]
+    subset = array.view(subset_dtype(array.dtype, to_use))
+    return repack_fields(subset)
+
+
+def subset_dtype(dtype, names):
+    '''
+    Create a new dtype only containing a subset of the columns
+    '''
+    formats = [dtype.fields[name][0] for name in names]
+    offsets = [dtype.fields[name][1] for name in names]
+    itemsize = dtype.itemsize
+    return np.dtype(dict(
+        names=names,
+        formats=formats,
+        offsets=offsets,
+        itemsize=itemsize
+    ))
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..4676ca16 100644
--- a/ctapipe/conftest.py
+++ b/ctapipe/conftest.py
@@ -133,9 +133,7 @@ def subarray_and_event_gamma_off_axis_500_gev(
 
 @pytest.fixture(scope="session")
 def prod5_gamma_simtel_path():
-    return get_dataset_path(
-        "gamma_20deg_0deg_run2___cta-prod5-paranal_desert-2147m-Paranal-dark_cone10-100evts.simtel.zst"
-    )
+    return get_dataset_path("gamma_prod5.simtel.zst")
 
 
 @pytest.fixture(scope="session")
@@ -176,8 +174,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",
-            "--max-events=20",
+            "--write-showers",
         ]
         assert run_tool(ProcessorTool(), argv=argv, cwd=dl2_tmp_path) == 0
         return output
@@ -202,7 +199,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 +225,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..85c631c0 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,76 @@ 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"
+def test_invalid_input():
+    class ExampleQualityQuery(QualityQuery):
+        """Available variables: x"""
 
-    assert len(query) == 4  # 4 events counted
+        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):
+        query(y=5)
 
-def test_bad_selector():
-    """ ensure failure if a selector function is not a function or can't be evaluated"""
 
+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", "3"),
-                ("smallish", "lambda x: x < 10"),
-            ]
-        )
-        assert s
-
-    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])
+
+
+def test_printing():
+    """Just check the query can be stringified correctly"""
+    query = QualityQuery(
+        quality_criteria=[("check", "x>3")],
+    )
+
+    assert isinstance(str(query), str)
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 ce5d0f24..a4242b0c 100644
--- a/ctapipe/io/tests/test_datawriter.py
+++ b/ctapipe/io/tests/test_datawriter.py
@@ -16,7 +16,7 @@ from ctapipe.utils import get_dataset_path
 
 
 def generate_dummy_dl2(event):
-    """ generate some dummy DL2 info and see if we can write it """
+    """generate some dummy DL2 info and see if we can write it"""
 
     algos = ["HillasReconstructor", "ImPACTReconstructor"]
 
@@ -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:
@@ -244,7 +242,9 @@ def test_metadata(tmpdir: Path):
                     "name": "Maximilian Nöthe",
                     "email": "maximilian.noethe@tu-dortmund.de",
                     "organization": "TU Dortmund",
-                }
+                },
+                "Instrument": {"site": "CTA-North", "id_": "alpha"},
+                "context_metadata": {"EXAMPLE": "test_value"},
             }
         }
     )
@@ -266,6 +266,9 @@ def test_metadata(tmpdir: Path):
             assert meta["CTA CONTACT NAME"] == "Maximilian Nöthe"
             assert meta["CTA CONTACT EMAIL"] == "maximilian.noethe@tu-dortmund.de"
             assert meta["CTA CONTACT ORGANIZATION"] == "TU Dortmund"
+            assert meta["CTA INSTRUMENT SITE"] == "CTA-North"
+            assert meta["CTA INSTRUMENT ID"] == "alpha"
+            assert meta["CONTEXT EXAMPLE"] == "test_value"
 
 
 def test_write_only_r1(r1_hdf5_file):
diff --git a/ctapipe/io/tests/test_hdf5eventsource.py b/ctapipe/io/tests/test_hdf5eventsource.py
index d3d98f83..9d792994 100644
--- a/ctapipe/io/tests/test_hdf5eventsource.py
+++ b/ctapipe/io/tests/test_hdf5eventsource.py
@@ -19,8 +19,8 @@ def test_metadata(dl1_file):
             DataLevel.DL1_IMAGES,
             DataLevel.DL1_PARAMETERS,
         }
-        assert list(source.obs_ids) == [2]
-        assert source.simulation_config[2].corsika_version == 7710
+        assert list(source.obs_ids) == [1]
+        assert source.simulation_config[1].corsika_version == 7710
 
 
 def test_subarray(dl1_file):
@@ -151,3 +151,16 @@ def test_trigger_allowed_tels(dl1_proton_file):
             assert len(e.trigger.tels_with_trigger) > 1
 
         assert i == 1
+
+
+def test_read_dl2(dl2_shower_geometry_file):
+    algorithm = "HillasReconstructor"
+
+    with HDF5EventSource(dl2_shower_geometry_file) as s:
+        e = next(iter(s))
+        assert algorithm in e.dl2.stereo.geometry
+        tel_mask = e.dl2.stereo.geometry[algorithm].tel_ids
+        tel_ids = s.subarray.tel_mask_to_tel_ids(tel_mask)
+        for tel_id in tel_ids:
+            assert algorithm in e.dl2.tel[tel_id].impact
+            assert e.dl2.tel[tel_id].impact[algorithm].distance is not None
diff --git a/ctapipe/io/tests/test_simteleventsource.py b/ctapipe/io/tests/test_simteleventsource.py
index 935a8582..49bc2d38 100644
--- a/ctapipe/io/tests/test_simteleventsource.py
+++ b/ctapipe/io/tests/test_simteleventsource.py
@@ -350,42 +350,39 @@ def test_only_config():
 
 def test_calibscale_and_calibshift(prod5_gamma_simtel_path):
 
-    telid = 25
-
     with SimTelEventSource(input_url=prod5_gamma_simtel_path, max_events=1) as source:
+        event = next(iter(source))
 
-        for event in source:
-            pass
+    # make sure we actually have data
+    assert len(event.r1.tel) > 0
 
     calib_scale = 2.0
 
     with SimTelEventSource(
         input_url=prod5_gamma_simtel_path, max_events=1, calib_scale=calib_scale
     ) as source:
+        event_scaled = next(iter(source))
 
-        for event_scaled in source:
-            pass
-
-    np.testing.assert_allclose(
-        event.r1.tel[telid].waveform[0],
-        event_scaled.r1.tel[telid].waveform[0] / calib_scale,
-        rtol=0.1,
-    )
+    for tel_id, r1 in event.r1.tel.items():
+        np.testing.assert_allclose(
+            r1.waveform[0],
+            event_scaled.r1.tel[tel_id].waveform[0] / calib_scale,
+            rtol=0.1,
+        )
 
     calib_shift = 2.0  # p.e.
 
     with SimTelEventSource(
         input_url=prod5_gamma_simtel_path, max_events=1, calib_shift=calib_shift
     ) as source:
+        event_shifted = next(iter(source))
 
-        for event_shifted in source:
-            pass
-
-    np.testing.assert_allclose(
-        event.r1.tel[telid].waveform[0],
-        event_shifted.r1.tel[telid].waveform[0] - calib_shift,
-        rtol=0.1,
-    )
+    for tel_id, r1 in event.r1.tel.items():
+        np.testing.assert_allclose(
+            r1.waveform[0],
+            event_shifted.r1.tel[tel_id].waveform[0] - calib_shift,
+            rtol=0.1,
+        )
 
 
 def test_true_image_sum():
diff --git a/ctapipe/io/tests/test_table_loader.py b/ctapipe/io/tests/test_table_loader.py
index b5b3e7bc..b4701be2 100644
--- a/ctapipe/io/tests/test_table_loader.py
+++ b/ctapipe/io/tests/test_table_loader.py
@@ -2,6 +2,8 @@ import pytest
 import tables
 import numpy as np
 
+from ctapipe.instrument.subarray import SubarrayDescription
+
 
 @pytest.fixture(params=["by_type", "by_id"])
 def test_file(request, dl1_file, dl1_by_type_file):
@@ -44,16 +46,17 @@ def test_telescope_events_for_tel_id(test_file):
     loader = TableLoader(dl1_file, load_dl1_parameters=True, load_trigger=True)
 
     with loader as table_loader:
-        table = table_loader.read_telescope_events([25])
+        table = table_loader.read_telescope_events([8])
         assert "hillas_length" in table.colnames
         assert "time" in table.colnames
         assert "event_type" in table.colnames
-        assert np.all(table["tel_id"] == 25)
+        assert np.all(table["tel_id"] == 8)
 
     with TableLoader(dl1_file, load_dl1_images=True) as table_loader:
-        table = table_loader.read_telescope_events([25])
+        table = table_loader.read_telescope_events([8])
         assert "image" in table.colnames
-        assert np.all(table["tel_id"] == 25)
+        assert np.all(table["tel_id"] == 8)
+        assert table["obs_id"].dtype == np.int32
 
     assert not table_loader.h5file.isopen
 
@@ -65,8 +68,8 @@ def test_load_instrument(test_file):
     _, dl1_file = test_file
 
     with TableLoader(dl1_file, load_instrument=True) as table_loader:
-        expected = table_loader.subarray.tel[25].optics.equivalent_focal_length
-        table = table_loader.read_telescope_events([25])
+        expected = table_loader.subarray.tel[8].optics.equivalent_focal_length
+        table = table_loader.read_telescope_events([8])
         assert "equivalent_focal_length" in table.colnames
         assert np.all(table["equivalent_focal_length"] == expected)
 
@@ -80,9 +83,11 @@ def test_load_simulated(test_file):
     with TableLoader(dl1_file, load_simulated=True) as table_loader:
         table = table_loader.read_subarray_events()
         assert "true_energy" in table.colnames
+        assert table["obs_id"].dtype == np.int32
 
-        table = table_loader.read_telescope_events([25])
+        table = table_loader.read_telescope_events([8])
         assert "true_energy" in table.colnames
+        assert "true_impact_distance" in table.colnames
 
 
 def test_true_images(test_file):
@@ -134,6 +139,8 @@ def test_read_telescope_events_type(test_file_dl2):
 
     _, dl2_file = test_file_dl2
 
+    subarray = SubarrayDescription.from_hdf(dl2_file)
+
     with TableLoader(
         dl2_file,
         load_dl1_images=False,
@@ -150,8 +157,10 @@ 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])
+        expected_ids = subarray.get_tel_ids_for_type("MST_MST_FlashCam")
+        assert set(table["tel_id"].data).issubset(expected_ids)
         assert "equivalent_focal_length" in table.colnames
+        assert "HillasReconstructor_impact_distance" in table.colnames
 
 
 def test_read_telescope_events_by_type(test_file_dl2):
@@ -160,6 +169,7 @@ def test_read_telescope_events_by_type(test_file_dl2):
     from ctapipe.io.tableloader import TableLoader
 
     _, dl2_file = test_file_dl2
+    subarray = SubarrayDescription.from_hdf(dl2_file)
 
     with TableLoader(
         dl2_file,
@@ -181,7 +191,8 @@ def test_read_telescope_events_by_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])
+            expected_ids = subarray.get_tel_ids_for_type(tel_type)
+            assert set(table["tel_id"].data).issubset(expected_ids)
             assert "equivalent_focal_length" in table.colnames
 
 
@@ -202,3 +213,48 @@ def test_h5file(test_file_dl2):
             assert 25 in loader.subarray.tel
             loader.read_subarray_events()
             loader.read_telescope_events()
+
+
+def test_chunked(dl2_shower_geometry_file):
+    """Test chunked reading"""
+    from ctapipe.io.tableloader import TableLoader, read_table
+
+    trigger = read_table(dl2_shower_geometry_file, "/dl1/event/subarray/trigger")
+    n_events = len(trigger)
+    n_read = 0
+
+    n_chunks = 2
+    chunk_size = int(np.ceil(n_events / n_chunks))
+
+    with TableLoader(
+        dl2_shower_geometry_file,
+        load_dl1_images=False,
+        load_true_images=False,
+        load_dl1_parameters=True,
+        load_dl2=True,
+        load_simulated=True,
+        load_trigger=True,
+    ) as table_loader:
+
+        tel_event_it = table_loader.read_telescope_events_chunked(chunk_size)
+        event_it = table_loader.read_subarray_events_chunked(chunk_size)
+
+        for chunk, (events, tel_events) in enumerate(zip(event_it, tel_event_it)):
+            n_read += len(events)
+
+            # last chunk might be smaller
+            if chunk == (n_chunks - 1):
+                assert len(events) == n_events % chunk_size
+            else:
+                assert len(events) == chunk_size
+
+            assert len(tel_events) == np.sum(events["tels_with_trigger"])
+
+            unique_events_from_tel_table = np.unique(tel_events[["obs_id", "event_id"]])
+
+            # test we get matching obs_ids / event_ids
+            idx = events[["obs_id", "event_id"]].argsort()
+            sorted_index = np.array(events[["obs_id", "event_id"]][idx])
+            assert (unique_events_from_tel_table == sorted_index).all()
+
+    assert n_read == n_events
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_merge.py b/ctapipe/tools/tests/test_merge.py
index 446bff21..3fc4512c 100644
--- a/ctapipe/tools/tests/test_merge.py
+++ b/ctapipe/tools/tests/test_merge.py
@@ -97,12 +97,17 @@ def test_skip_images(tmp_path, dl1_file, dl1_proton_file):
         ],
         cwd=tmp_path,
     )
+    assert ret == 0
 
     with tables.open_file(output, "r") as f:
         assert "images" not in f.root.dl1.event.telescope
+        assert "images" in f.root.simulation.event.telescope
         assert "parameters" in f.root.dl1.event.telescope
 
-    assert ret == 0
+    t = read_table(output, "/simulation/event/telescope/images/tel_001")
+    assert "true_image" not in t.colnames
+    assert "true_image_sum" in t.colnames
+
 
 
 def test_allowed_tels(tmp_path, dl1_file, dl1_proton_file):
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
     )
 
 
diff --git a/ctapipe/utils/tests/test_arrays.py b/ctapipe/utils/tests/test_arrays.py
new file mode 100644
index 00000000..2d780f21
--- /dev/null
+++ b/ctapipe/utils/tests/test_arrays.py
@@ -0,0 +1,16 @@
+"""Tests for the array utils"""
+import numpy as np
+
+
+def test_drop_columns():
+    """Test drop columns"""
+    from ctapipe.utils.arrays import recarray_drop_columns
+
+    array = np.empty(
+        1000, dtype=[("x", float), ("N", int), ("valid", bool), ("tel_id", np.uint16)]
+    )
+    assert array.dtype.itemsize == 19
+
+    dropped = recarray_drop_columns(array, columns=["N", "tel_id"])
+    assert dropped.dtype.names == ("x", "valid")
+    assert dropped.dtype.itemsize == 9