Repository: cta-observatory/ctapipe (python)

Issue:
Container.validate() doesn't validate sub-containers
**Describe the bug**

Containers contained in Containers do not get validated with the parent container's `.validate()` method is called.

**To Reproduce**
```py3
from ctapipe.core import Container, Field

class Sub(Container):
    y = Field(0*u.m, "Sub.y", unit=u.m)

class A(Container):
    x = Field(0*u.m, "A.x", unit=u.m)
    sub = Field(None, "sub", type=Sub)

# works:
try: 
    A(x=10*u.s, sub=Sub(y=10*u.m)).validate()
except Exception as err:
    print(err)

#> A Field 'x': the value '10.0 s' (<class 'astropy.units.quantity.Quantity'>) is invalid: : 's' (time) and 'm' (length) are not convertible

# doesn't not raise exception (the sub field doesn't get checked recursively):
try: 
    A(x=10*u.m, sub=Sub(y=10*u.s)).validate()
except Exception as err:
    print(err)
```

**Expected behavior**
Should raise exception if any sub-Container also fails validation




Pull request description:
recursively validate containers
fixes #1975 

Patch:
diff --git a/ctapipe/calib/camera/flatfield.py b/ctapipe/calib/camera/flatfield.py
index 18411057..1e3263c6 100644
--- a/ctapipe/calib/camera/flatfield.py
+++ b/ctapipe/calib/camera/flatfield.py
@@ -115,18 +115,18 @@ class FlatFieldCalculator(Component):
 
 class FlasherFlatFieldCalculator(FlatFieldCalculator):
     """Calculates flat-field parameters from flasher data
-       based on the best algorithm described by S. Fegan in MST-CAM-TN-0060 (eq. 19)
-       Pixels are defined as outliers on the base of a cut on the pixel charge median
-       over the full sample distribution and the pixel signal time inside the
-       waveform time
+      based on the best algorithm described by S. Fegan in MST-CAM-TN-0060 (eq. 19)
+      Pixels are defined as outliers on the base of a cut on the pixel charge median
+      over the full sample distribution and the pixel signal time inside the
+      waveform time
 
 
-     Parameters
-     ----------
-     charge_cut_outliers : List[2]
-         Interval of accepted charge values (fraction with respect to camera median value)
-     time_cut_outliers : List[2]
-         Interval (in waveform samples) of accepted time values
+    Parameters
+    ----------
+    charge_cut_outliers : List[2]
+        Interval of accepted charge values (fraction with respect to camera median value)
+    time_cut_outliers : List[2]
+        Interval (in waveform samples) of accepted time values
 
     """
 
@@ -140,18 +140,18 @@ class FlasherFlatFieldCalculator(FlatFieldCalculator):
 
     def __init__(self, **kwargs):
         """Calculates flat-field parameters from flasher data
-           based on the best algorithm described by S. Fegan in MST-CAM-TN-0060 (eq. 19)
-           Pixels are defined as outliers on the base of a cut on the pixel charge median
-           over the full sample distribution and the pixel signal time inside the
-           waveform time
+          based on the best algorithm described by S. Fegan in MST-CAM-TN-0060 (eq. 19)
+          Pixels are defined as outliers on the base of a cut on the pixel charge median
+          over the full sample distribution and the pixel signal time inside the
+          waveform time
 
 
-         Parameters:
-         ----------
-         charge_cut_outliers : List[2]
-             Interval of accepted charge values (fraction with respect to camera median value)
-         time_cut_outliers : List[2]
-             Interval (in waveform samples) of accepted time values
+        Parameters:
+        ----------
+        charge_cut_outliers : List[2]
+            Interval of accepted charge values (fraction with respect to camera median value)
+        time_cut_outliers : List[2]
+            Interval (in waveform samples) of accepted time values
 
         """
         super().__init__(**kwargs)
@@ -190,14 +190,14 @@ class FlasherFlatFieldCalculator(FlatFieldCalculator):
 
     def calculate_relative_gain(self, event):
         """
-         calculate the flatfield statistical values
-         and fill mon.tel[tel_id].flatfield container
+        calculate the flatfield statistical values
+        and fill mon.tel[tel_id].flatfield container
 
-         Parameters
-         ----------
-         event : general event container
+        Parameters
+        ----------
+        event : general event container
 
-         """
+        """
 
         # initialize the np array at each cycle
         waveform = event.r1.tel[self.tel_id].waveform
@@ -235,7 +235,7 @@ class FlasherFlatFieldCalculator(FlatFieldCalculator):
 
         self.collect_sample(dl1.image, pixel_mask, dl1.peak_time)
 
-        sample_age = trigger_time - self.time_start
+        sample_age = (trigger_time - self.time_start).to_value(u.s)
 
         # check if to create a calibration event
         if (
@@ -296,7 +296,7 @@ class FlasherFlatFieldCalculator(FlatFieldCalculator):
     def calculate_time_results(
         self, trace_time, masked_pixels_of_sample, time_start, trigger_time
     ):
-        """Calculate and return the time results """
+        """Calculate and return the time results"""
         masked_trace_time = np.ma.array(trace_time, mask=masked_pixels_of_sample)
 
         # median over the sample per pixel
diff --git a/ctapipe/calib/camera/pedestals.py b/ctapipe/calib/camera/pedestals.py
index 3bdcaceb..ed8c527d 100644
--- a/ctapipe/calib/camera/pedestals.py
+++ b/ctapipe/calib/camera/pedestals.py
@@ -72,8 +72,7 @@ class PedestalCalculator(Component):
         Set to None if no configuration to pass.
 
     kwargs
-
-"""
+    """
 
     tel_id = Int(0, help="id of the telescope to calculate the pedestal values").tag(
         config=True
@@ -115,7 +114,7 @@ class PedestalCalculator(Component):
 
         kwargs
 
-    """
+        """
 
         super().__init__(**kwargs)
 
@@ -143,21 +142,21 @@ class PedestalCalculator(Component):
 
 class PedestalIntegrator(PedestalCalculator):
     """Calculates pedestal parameters integrating the charge of pedestal events:
-       the pedestal value corresponds to the charge estimated with the selected
-       charge extractor
-       The pixels are set as outliers on the base of a cut on the pixel charge median
-       over the pedestal sample and the pixel charge standard deviation over
-       the pedestal sample with respect to the camera median values
+      the pedestal value corresponds to the charge estimated with the selected
+      charge extractor
+      The pixels are set as outliers on the base of a cut on the pixel charge median
+      over the pedestal sample and the pixel charge standard deviation over
+      the pedestal sample with respect to the camera median values
 
 
-     Parameters:
-     ----------
-     charge_median_cut_outliers : List[2]
-         Interval (number of std) of accepted charge values around camera median value
-     charge_std_cut_outliers : List[2]
-         Interval (number of std) of accepted charge standard deviation around camera median value
+    Parameters:
+    ----------
+    charge_median_cut_outliers : List[2]
+        Interval (number of std) of accepted charge values around camera median value
+    charge_std_cut_outliers : List[2]
+        Interval (number of std) of accepted charge standard deviation around camera median value
 
-     """
+    """
 
     charge_median_cut_outliers = List(
         [-3, 3],
@@ -170,19 +169,19 @@ class PedestalIntegrator(PedestalCalculator):
 
     def __init__(self, **kwargs):
         """Calculates pedestal parameters integrating the charge of pedestal events:
-           the pedestal value corresponds to the charge estimated with the selected
-           charge extractor
-           The pixels are set as outliers on the base of a cut on the pixel charge median
-           over the pedestal sample and the pixel charge standard deviation over
-           the pedestal sample with respect to the camera median values
-
-
-         Parameters:
-         ----------
-         charge_median_cut_outliers : List[2]
-             Interval (number of std) of accepted charge values around camera median value
-         charge_std_cut_outliers : List[2]
-             Interval (number of std) of accepted charge standard deviation around camera median value
+          the pedestal value corresponds to the charge estimated with the selected
+          charge extractor
+          The pixels are set as outliers on the base of a cut on the pixel charge median
+          over the pedestal sample and the pixel charge standard deviation over
+          the pedestal sample with respect to the camera median values
+
+
+        Parameters:
+        ----------
+        charge_median_cut_outliers : List[2]
+            Interval (number of std) of accepted charge values around camera median value
+        charge_std_cut_outliers : List[2]
+            Interval (number of std) of accepted charge standard deviation around camera median value
         """
 
         super().__init__(**kwargs)
@@ -257,7 +256,7 @@ class PedestalIntegrator(PedestalCalculator):
 
         self.collect_sample(dl1.image, pixel_mask)
 
-        sample_age = trigger_time - self.time_start
+        sample_age = (trigger_time - self.time_start).to_value(u.s)
 
         # check if to create a calibration event
         if (
diff --git a/ctapipe/core/container.py b/ctapipe/core/container.py
index cb1b536c..faa6b325 100644
--- a/ctapipe/core/container.py
+++ b/ctapipe/core/container.py
@@ -1,10 +1,11 @@
 from collections import defaultdict
 from copy import deepcopy
 from pprint import pformat
-from textwrap import wrap
+from textwrap import wrap, dedent
 import warnings
 import numpy as np
 from astropy.units import UnitConversionError, Quantity, Unit
+from inspect import isclass
 
 import logging
 
@@ -24,25 +25,25 @@ class Field:
 
     Parameters
     ----------
-    default:
-        default value of the item (this will be set when the `Container`
+    default :
+        default value of the item. This will be set when the `Container`
         is constructed, as well as when  ``Container.reset`` is called
-    description: str
+    description : str
         Help text associated with the item
-    unit: str or astropy.units.core.UnitBase
+    unit : str or astropy.units.core.UnitBase
         unit to convert to when writing output, or None for no conversion
-    ucd: str
+    ucd : str
         universal content descriptor (see Virtual Observatory standards)
-    type: type
+    type : type
         expected type of value
-    dtype: str or np.dtype
+    dtype : str or np.dtype
         expected data type of the value, None to ignore in validation.
         Means value is expected to be a numpy array or astropy quantity
-    ndim: int or None
+    ndim : int or None
         expected dimensionality of the data, for arrays, None to ignore
-    allow_none:
+    allow_none : bool
         if the value of None is given to this Field, skip validation
-    max_len:
+    max_len : int
         if type is str, max_len is the maximum number of bytes of the utf-8
         encoded string to be used.
     """
@@ -59,7 +60,6 @@ class Field:
         allow_none=True,
         max_length=None,
     ):
-
         self.default = default
         self.description = description
         self.unit = Unit(unit) if unit is not None else None
@@ -71,7 +71,34 @@ class Field:
         self.max_length = max_length
 
     def __repr__(self):
-        desc = f"{self.description}"
+        if isinstance(self.default, Container):
+            default = f"{self.default.__class__.__name__}"
+        elif isinstance(self.default, Map):
+            if isclass(self.default.default_factory):
+                cls = self.default.default_factory
+                default = f"Map({cls.__module__}.{cls.__name__})"
+            else:
+                default = f"Map({repr(self.default.default_factory)}"
+        else:
+            default = str(self.default)
+        cmps = [f"Field(default={default}"]
+        if self.unit is not None:
+            cmps.append(f", unit={self.unit}")
+        if self.dtype is not None:
+            cmps.append(f", dtype={self.dtype}")
+        if self.ndim is not None:
+            cmps.append(f", ndim={self.ndim}")
+        if self.type is not None:
+            cmps.append(f", type={self.type.__name__}")
+        if self.allow_none is False:
+            cmps.append(", allow_none=False")
+        if self.max_length is not None:
+            cmps.append(f", max_length={self.max_length}")
+        cmps.append(")")
+        return "".join(cmps)
+
+    def __str__(self):
+        desc = f"{self.description} with default {self.default}"
         if self.unit is not None:
             desc += f" [{self.unit}]"
         if self.ndim is not None:
@@ -108,6 +135,20 @@ class Field:
                 f"{errorstr} Should be an instance of {self.type}"
             )
 
+        if isinstance(value, Container):
+            # recursively check sub-containers
+            value.validate()
+            return
+
+        if isinstance(value, Map):
+            for key, map_value in value.items():
+                if isinstance(map_value, Container):
+                    try:
+                        map_value.validate()
+                    except FieldValidationError as err:
+                        raise FieldValidationError(f"[{key}]: {err} ")
+            return
+
         if self.unit is not None:
             if not isinstance(value, Quantity):
                 raise FieldValidationError(
@@ -153,6 +194,23 @@ class DeprecatedField(Field):
         self.reason = reason
 
 
+_doc_template = """{doc}
+
+Attributes
+----------
+{fields}
+meta : dict
+    dict of attached metadata
+prefix : str
+    Prefix attached to column names when saved to a table or file
+"""
+
+
+def _build_docstring(doc, fields):
+    fields = [f"{k} : {f!r}\n    {f.description}" for k, f in fields.items()]
+    return _doc_template.format(doc=dedent(doc), fields="\n".join(fields))
+
+
 class ContainerMeta(type):
     """
     The MetaClass for `Container`
@@ -179,6 +237,8 @@ class ContainerMeta(type):
         for k in field_names:
             dct["fields"][k] = dct.pop(k)
 
+        dct["__doc__"] = _build_docstring(dct.get("__doc__", ""), dct["fields"])
+
         new_cls = type.__new__(cls, name, bases, dct)
 
         # if prefix was not set as a class variable, build a default one
@@ -352,7 +412,7 @@ class Container(metaclass=ContainerMeta):
                 extra = ".*"
             if isinstance(getattr(self, name), Map):
                 extra = "[*]"
-            desc = "{:>30s}: {}".format(name + extra, repr(item))
+            desc = "{:>30s}: {}".format(name + extra, str(item))
             lines = wrap(desc, 80, subsequent_indent=" " * 32)
             text.extend(lines)
         return "\n".join(text)
diff --git a/ctapipe/image/hillas.py b/ctapipe/image/hillas.py
index 8fdd405a..fae58ccb 100644
--- a/ctapipe/image/hillas.py
+++ b/ctapipe/image/hillas.py
@@ -7,7 +7,6 @@ Hillas-style moment-based shower image parametrization.
 import astropy.units as u
 import numpy as np
 from astropy.coordinates import Angle
-from astropy.units import Quantity
 from ..containers import CameraHillasParametersContainer, HillasParametersContainer
 
 
@@ -105,12 +104,15 @@ def hillas_parameters(geom, image):
         container of hillas parametesr
     """
     unit = geom.pix_x.unit
-    pix_x = Quantity(np.asanyarray(geom.pix_x, dtype=np.float64)).value
-    pix_y = Quantity(np.asanyarray(geom.pix_y, dtype=np.float64)).value
+    pix_x = geom.pix_x.to_value(unit)
+    pix_y = geom.pix_y.to_value(unit)
     image = np.asanyarray(image, dtype=np.float64)
-    image = np.ma.filled(image, 0)
-    msg = "Image and pixel shape do not match"
-    assert pix_x.shape == pix_y.shape == image.shape, msg
+
+    if isinstance(image, np.ma.masked_array):
+        image = np.ma.filled(image, 0)
+
+    if not (pix_x.shape == pix_y.shape == image.shape):
+        raise ValueError("Image and pixel shape do not match")
 
     size = np.sum(image)
 
@@ -158,11 +160,11 @@ def hillas_parameters(geom, image):
         # calculate higher order moments along shower axes
         longitudinal = delta_x * np.cos(psi) + delta_y * np.sin(psi)
 
-        m3_long = np.average(longitudinal ** 3, weights=image)
-        skewness_long = m3_long / length ** 3
+        m3_long = np.average(longitudinal**3, weights=image)
+        skewness_long = m3_long / length**3
 
-        m4_long = np.average(longitudinal ** 4, weights=image)
-        kurtosis_long = m4_long / length ** 4
+        m4_long = np.average(longitudinal**4, weights=image)
+        kurtosis_long = m4_long / length**4
 
     # Compute of the Hillas parameters uncertainties.
     # Implementation described in [hillas_uncertainties]_ This is an internal MAGIC document
@@ -174,8 +176,8 @@ def hillas_parameters(geom, image):
     b = (1 - cos_2psi) / 2
     c = np.sin(2 * psi)
 
-    A = ((delta_x ** 2.0) - cov[0][0]) / size
-    B = ((delta_y ** 2.0) - cov[1][1]) / size
+    A = ((delta_x**2.0) - cov[0][0]) / size
+    B = ((delta_y**2.0) - cov[1][1]) / size
     C = ((delta_x * delta_y) - cov[0][1]) / size
 
     # Hillas's uncertainties
diff --git a/ctapipe/io/tableloader.py b/ctapipe/io/tableloader.py
index 6bd6c6f9..39fe12c6 100644
--- a/ctapipe/io/tableloader.py
+++ b/ctapipe/io/tableloader.py
@@ -97,14 +97,6 @@ def _empty_telescope_events_table():
     return Table(names=TELESCOPE_EVENT_KEYS, dtype=[np.int32, np.int64, np.int16])
 
 
-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.int32, np.int64])
-
-
 def _get_structure(h5file):
     """
     Check if data is stored by telescope id or by telescope type in h5file.
@@ -135,12 +127,18 @@ def _add_column_prefix(table, prefix, ignore=()):
 
 def _join_subarray_events(table1, table2):
     """Outer join two tables on the telescope subarray keys"""
-    return join_allow_empty(table1, table2, SUBARRAY_EVENT_KEYS, "outer")
+    return join_allow_empty(table1, table2, SUBARRAY_EVENT_KEYS, "left")
 
 
 def _join_telescope_events(table1, table2):
     """Outer join two tables on the telescope event keys"""
-    return join_allow_empty(table1, table2, TELESCOPE_EVENT_KEYS, "outer")
+    # we start with an empty table, but after the first non-empty, we perform
+    # left joins
+    if len(table1) == 0:
+        how = "right"
+    else:
+        how = "left"
+    return join_allow_empty(table1, table2, TELESCOPE_EVENT_KEYS, how)
 
 
 class TableLoader(Component):
@@ -188,9 +186,7 @@ class TableLoader(Component):
     load_true_parameters = traits.Bool(
         False, help="load image parameters obtained from true images"
     ).tag(config=True)
-    load_trigger = traits.Bool(True, help="load subarray trigger information").tag(
-        config=True
-    )
+
     load_instrument = traits.Bool(
         False, help="join subarray instrument information to each event"
     ).tag(config=True)
@@ -242,6 +238,10 @@ class TableLoader(Component):
     def __exit__(self, exc_type, exc_value, traceback):
         self.close()
 
+    def __len__(self):
+        """Number of subarray events in input file"""
+        return self.h5file.root[TRIGGER_TABLE].shape[0]
+
     def _read_telescope_table(self, group, tel_id, start=None, stop=None):
         if self.structure == "by_id":
             key = f"{group}/tel_{tel_id:03d}"
@@ -259,7 +259,37 @@ class TableLoader(Component):
 
         return table
 
-    def read_subarray_events(self, start=None, stop=None):
+    def _get_sort_index(self, start=None, stop=None):
+        """
+        Get plain index of increasing integers in the order in the file.
+
+        Astropy.table.join orders by the join keys, we want to keep
+        the original order. Adding an index before the first join and then
+        using that in the end to sort back to the original order solves this.
+        """
+        table = read_table(
+            self.h5file,
+            TRIGGER_TABLE,
+            start=start,
+            stop=stop,
+        )[["obs_id", "event_id"]]
+        self._add_index_if_needed(table)
+        return table
+
+    @staticmethod
+    def _sort_to_original_order(table, include_tel_id=False):
+        if include_tel_id:
+            table.sort(("__index__", "tel_id"))
+        else:
+            table.sort("__index__")
+        table.remove_column("__index__")
+
+    @staticmethod
+    def _add_index_if_needed(table):
+        if "__index__" not in table.colnames:
+            table["__index__"] = np.arange(len(table))
+
+    def read_subarray_events(self, start=None, stop=None, keep_order=True):
         """Read subarray-based event information.
 
         Returns
@@ -267,11 +297,9 @@ class TableLoader(Component):
         table: astropy.io.Table
             Table with primary index columns "obs_id" and "event_id".
         """
-        table = _empty_subarray_events_table()
-
-        if self.load_trigger:
-            trigger = read_table(self.h5file, TRIGGER_TABLE, start=start, stop=stop)
-            table = _join_subarray_events(table, trigger)
+        table = read_table(self.h5file, TRIGGER_TABLE, start=start, stop=stop)
+        if keep_order:
+            self._add_index_if_needed(table)
 
         if self.load_simulated and SHOWER_TABLE in self.h5file:
             showers = read_table(self.h5file, SHOWER_TABLE, start=start, stop=stop)
@@ -297,6 +325,9 @@ class TableLoader(Component):
                             dl2, prefix=algorithm, ignore=SUBARRAY_EVENT_KEYS
                         )
                         table = _join_subarray_events(table, dl2)
+
+        if keep_order:
+            self._sort_to_original_order(table)
         return table
 
     def read_subarray_events_chunked(self, chunk_size):
@@ -308,11 +339,9 @@ class TableLoader(Component):
         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,
+            n_total=len(self),
             chunk_size=chunk_size,
         )
 
@@ -392,35 +421,61 @@ class TableLoader(Component):
 
         if self.load_simulated and TRUE_IMPACT_GROUP in self.h5file.root:
             impacts = self._read_telescope_table(
-                TRUE_IMPACT_GROUP, tel_id, start=start, stop=stop
+                TRUE_IMPACT_GROUP,
+                tel_id,
+                start=start,
+                stop=stop,
             )
             table = _join_telescope_events(table, impacts)
 
         return table
 
     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)
+        tel_start = tel_start if tel_start is not None else [None] * len(tel_ids)
+        tel_stop = tel_stop if tel_stop is not None else [None] * len(tel_ids)
 
-        if tel_stop is None:
-            tel_stop = [None] * len(tel_ids)
+        tables = []
+        for tel_id, start, stop in zip(tel_ids, tel_start, tel_stop):
+            # no events for this telescope in chunk
+            if start is not None and stop is not None and (stop - start) == 0:
+                continue
 
-        table = vstack(
-            [
+            tables.append(
                 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
+        return vstack(tables)
 
     def _join_subarray_info(self, table, start=None, stop=None):
-        subarray_events = self.read_subarray_events(start=start, stop=stop)
+        subarray_events = self.read_subarray_events(
+            start=start,
+            stop=stop,
+            keep_order=False,
+        )
         table = join_allow_empty(
             table, subarray_events, keys=SUBARRAY_EVENT_KEYS, join_type="left"
         )
         return table
 
+    def _get_tel_start_stop(self, tel_ids, start, stop):
+        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]
+
+        return tel_start, tel_stop
+
     def read_telescope_events(self, telescopes=None, start=None, stop=None):
         """
         Read telescope-based event information.
@@ -455,26 +510,15 @@ class TableLoader(Component):
         else:
             tel_ids = self.subarray.get_tel_ids(telescopes)
 
-        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]
-
+        tel_start, tel_stop = self._get_tel_start_stop(tel_ids, start, stop)
         table = self._read_telescope_events_for_ids(tel_ids, tel_start, tel_stop)
+        table = self._join_subarray_info(table, start=start, stop=stop)
 
-        if any([self.load_trigger, self.load_simulated, self.load_dl2]):
-            table = self._join_subarray_info(table, start=start, stop=stop)
+        # sort back to order in the file
+        table = _join_subarray_events(
+            table, self._get_sort_index(start=start, stop=stop)
+        )
+        self._sort_to_original_order(table, include_tel_id=True)
 
         return table
 
@@ -491,11 +535,9 @@ class TableLoader(Component):
 
         **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,
+            n_total=len(self),
             chunk_size=chunk_size,
             **kwargs,
         )
@@ -519,7 +561,9 @@ class TableLoader(Component):
         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]:
+    def read_telescope_events_by_type(
+        self, telescopes=None, start=None, stop=None
+    ) -> Dict[str, Table]:
         """Read telescope-based event information.
 
         Parameters
@@ -535,20 +579,54 @@ 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)
 
+        tel_start, tel_stop = self._get_tel_start_stop(tel_ids, start, stop)
+        tel_start = tel_start if tel_start is not None else [None] * len(tel_ids)
+        tel_stop = tel_stop if tel_stop is not None else [None] * len(tel_ids)
+
         by_type = defaultdict(list)
+        sort_index = self._get_sort_index(start=start, stop=stop)
+
+        for tel_id, start, stop in zip(tel_ids, tel_start, tel_stop):
+            # no events for this telescope in range start/stop
+            if start is not None and stop is not None and (stop - start) == 0:
+                continue
 
-        for tel_id in tel_ids:
             key = str(self.subarray.tel[tel_id])
-            by_type[key].append(self._read_telescope_events_for_id(tel_id))
+            by_type[key].append(
+                self._read_telescope_events_for_id(tel_id, start=start, stop=stop)
+            )
 
         by_type = {k: vstack(ts) for k, ts in by_type.items()}
 
-        if any([self.load_trigger, self.load_simulated, self.load_dl2]):
-            for key, table in by_type.items():
-                by_type[key] = self._join_subarray_info(table)
+        for key in by_type.keys():
+            by_type[key] = self._join_subarray_info(
+                by_type[key], start=start, stop=stop
+            )
+            by_type[key] = _join_subarray_events(by_type[key], sort_index)
+            self._sort_to_original_order(by_type[key], include_tel_id=True)
 
         return by_type
+
+    def read_telescope_events_by_type_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`
+        """
+        return ChunkIterator(
+            self.read_telescope_events_by_type,
+            n_total=len(self),
+            chunk_size=chunk_size,
+            **kwargs,
+        )
diff --git a/examples/load_one_event.py b/examples/load_one_event.py
index facf5cca..cf17b384 100644
--- a/examples/load_one_event.py
+++ b/examples/load_one_event.py
@@ -15,9 +15,11 @@ if __name__ == "__main__":
     if len(sys.argv) >= 2:
         filename = sys.argv[1]
     else:
-        filename = get_dataset_path("gamma_test_large.simtel.gz")
+        filename = get_dataset_path(
+            "gamma_LaPalma_baseline_20Zd_180Az_prod3b_test.simtel.gz"
+        )
 
-    with EventSource(filename, max_events=1) as source:
+    with EventSource(filename, max_events=1, focal_length_choice="nominal") as source:
         calib = CameraCalibrator(subarray=source.subarray)
         process_images = ImageProcessor(subarray=source.subarray)
         process_shower = ShowerProcessor(subarray=source.subarray)


Test patch:
diff --git a/ctapipe/conftest.py b/ctapipe/conftest.py
index 4676ca16..8742faed 100644
--- a/ctapipe/conftest.py
+++ b/ctapipe/conftest.py
@@ -233,6 +233,30 @@ def dl2_shower_geometry_file_type(dl2_tmp_path, prod5_gamma_simtel_path):
         return output
 
 
+@pytest.fixture(scope="session")
+def dl2_merged_file(dl2_tmp_path, dl2_shower_geometry_file, dl2_proton_geometry_file):
+    """
+    File containing both parameters and shower geometry from a gamma simulation set.
+    """
+    from ctapipe.core import run_tool
+    from ctapipe.tools.merge import MergeTool
+
+    output = dl2_tmp_path / "merged.training.h5"
+
+    # prevent running process multiple times in case of parallel tests
+    with FileLock(output.with_suffix(output.suffix + ".lock")):
+        if output.is_file():
+            return output
+
+        argv = [
+            f"--output={output}",
+            str(dl2_proton_geometry_file),
+            str(dl2_shower_geometry_file),
+        ]
+        assert run_tool(MergeTool(), argv=argv, cwd=dl2_tmp_path) == 0
+        return output
+
+
 @pytest.fixture(scope="session")
 def dl1_file(dl1_tmp_path, prod5_gamma_simtel_path):
     """
diff --git a/ctapipe/core/tests/test_container.py b/ctapipe/core/tests/test_container.py
index 0c0e0a01..1dc1409e 100644
--- a/ctapipe/core/tests/test_container.py
+++ b/ctapipe/core/tests/test_container.py
@@ -255,7 +255,7 @@ def test_field_validation():
 
 
 def test_container_validation():
-    """ check that we can validate all fields in a container"""
+    """check that we can validate all fields in a container"""
 
     class MyContainer(Container):
         x = Field(3.2, "test", unit="m")
@@ -264,4 +264,29 @@ def test_container_validation():
     with pytest.raises(FieldValidationError):
         MyContainer().validate()  # fails since 3.2 has no units
 
+    with pytest.raises(FieldValidationError):
+        MyContainer(x=10 * u.s).validate()  # seconds is not convertable to meters
+
     MyContainer(x=6.4 * u.m).validate()  # works
+
+
+def test_recursive_validation():
+    """
+    Check both sub-containers and Maps work with recursive validation
+    """
+
+    class ChildContainer(Container):
+        x = Field(3.2 * u.m, "test", unit="m")
+
+    class ParentContainer(Container):
+        cont = Field(None, "test sub", type=ChildContainer)
+        map = Field(Map(ChildContainer), "many children")
+
+    with pytest.raises(FieldValidationError):
+        ParentContainer(cont=ChildContainer(x=1 * u.s)).validate()
+
+    with pytest.raises(FieldValidationError):
+        cont = ParentContainer(cont=ChildContainer(x=1 * u.m))
+        cont.map[0] = ChildContainer(x=1 * u.m)
+        cont.map[1] = ChildContainer(x=1 * u.s)
+        cont.validate()
diff --git a/ctapipe/io/tests/test_table_loader.py b/ctapipe/io/tests/test_table_loader.py
index b4701be2..b8a63ff8 100644
--- a/ctapipe/io/tests/test_table_loader.py
+++ b/ctapipe/io/tests/test_table_loader.py
@@ -1,6 +1,67 @@
 import pytest
 import tables
 import numpy as np
+from astropy.table import Table
+
+from ctapipe.io.astropy_helpers import read_table
+
+
+def check_equal_array_event_order(table1, table2):
+    """
+    Check that values and order of array events is consistent in two tables.
+
+    Works for both tables of subarray events (obs_id, tel_id) and
+    tables of telescope events (obs_id, tel_id, event_id) and combinations
+    of the two.
+    """
+
+    def unique_events(table):
+        return np.unique(np.array(table[["obs_id", "event_id"]]), return_index=True)
+
+    unique_events1, indicies1 = unique_events(table1)
+    unique_events2, indicies2 = unique_events(table2)
+
+    if len(unique_events1) != len(unique_events2):
+        raise ValueError("Tables have different numbers of events")
+
+    if np.any(unique_events1 != unique_events2):
+        raise ValueError("Tables have different subarray events")
+
+    # we expect the rank in which the events appear in indicies to be the
+    # same, that means argsort should produce the same result on both:
+    order1 = np.argsort(indicies1)
+    order2 = np.argsort(indicies2)
+
+    if np.any(order1 != order2):
+        raise ValueError("Tables have subarray events in different order")
+
+
+def test_check_order():
+    with pytest.raises(ValueError, match="Tables have different numbers"):
+        check_equal_array_event_order(
+            Table({"obs_id": [1, 1, 2, 2, 3], "event_id": [1, 2, 1, 2, 1]}),
+            Table({"obs_id": [1, 1, 2, 2], "event_id": [1, 2, 1, 2]}),
+        )
+
+    with pytest.raises(ValueError, match="Tables have different subarray events"):
+        check_equal_array_event_order(
+            Table({"obs_id": [1, 1, 2, 2, 3], "event_id": [1, 2, 1, 2, 1]}),
+            Table({"obs_id": [1, 1, 2, 2, 4], "event_id": [1, 2, 1, 2, 1]}),
+        )
+
+    with pytest.raises(
+        ValueError, match="Tables have subarray events in different order"
+    ):
+        check_equal_array_event_order(
+            Table({"obs_id": [1, 1, 3, 2, 2], "event_id": [1, 2, 1, 1, 2]}),
+            Table({"obs_id": [1, 1, 2, 2, 3], "event_id": [1, 2, 1, 2, 1]}),
+        )
+
+    check_equal_array_event_order(
+        Table({"obs_id": [1, 1, 3, 2, 2], "event_id": [1, 2, 1, 1, 2]}),
+        Table({"obs_id": [1, 1, 3, 2, 2], "event_id": [1, 2, 1, 1, 2]}),
+    )
+
 
 from ctapipe.instrument.subarray import SubarrayDescription
 
@@ -43,7 +104,7 @@ def test_telescope_events_for_tel_id(test_file):
 
     _, dl1_file = test_file
 
-    loader = TableLoader(dl1_file, load_dl1_parameters=True, load_trigger=True)
+    loader = TableLoader(dl1_file, load_dl1_parameters=True)
 
     with loader as table_loader:
         table = table_loader.read_telescope_events([8])
@@ -118,20 +179,36 @@ def test_true_parameters(test_file):
 
 def test_read_subarray_events(test_file_dl2):
     """Test reading subarray events"""
-
     from ctapipe.io.tableloader import TableLoader
 
     _, dl2_file = test_file_dl2
 
-    with TableLoader(
-        dl2_file, load_dl2=True, load_simulated=True, load_trigger=True
-    ) as table_loader:
+    with TableLoader(dl2_file, load_dl2=True, load_simulated=True) as table_loader:
         table = table_loader.read_subarray_events()
         assert "HillasReconstructor_alt" in table.colnames
         assert "true_energy" in table.colnames
         assert "time" in table.colnames
 
 
+def test_table_loader_keeps_original_order(dl2_merged_file):
+    """Test reading subarray events keeps order in file"""
+    from ctapipe.io.tableloader import TableLoader
+
+    # check that the order is the same as in the file itself
+    trigger = read_table(dl2_merged_file, "/dl1/event/subarray/trigger")
+    # check we actually have unsorted input
+    assert not np.all(np.diff(trigger["obs_id"]) >= 0)
+
+    with TableLoader(
+        dl2_merged_file, load_dl2=True, load_simulated=True
+    ) as table_loader:
+        events = table_loader.read_subarray_events()
+        tel_events = table_loader.read_telescope_events()
+
+    check_equal_array_event_order(events, trigger)
+    check_equal_array_event_order(events, tel_events)
+
+
 def test_read_telescope_events_type(test_file_dl2):
     """Test reading telescope events for a given telescope type"""
 
@@ -148,7 +225,6 @@ def test_read_telescope_events_type(test_file_dl2):
         load_dl2=True,
         load_simulated=True,
         load_true_images=True,
-        load_trigger=False,
         load_instrument=True,
     ) as table_loader:
 
@@ -178,7 +254,6 @@ def test_read_telescope_events_by_type(test_file_dl2):
         load_dl2=True,
         load_simulated=True,
         load_true_images=True,
-        load_trigger=False,
         load_instrument=True,
     ) as table_loader:
 
@@ -225,6 +300,8 @@ def test_chunked(dl2_shower_geometry_file):
 
     n_chunks = 2
     chunk_size = int(np.ceil(n_events / n_chunks))
+    start = 0
+    stop = chunk_size
 
     with TableLoader(
         dl2_shower_geometry_file,
@@ -233,14 +310,18 @@ def test_chunked(dl2_shower_geometry_file):
         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)
+        by_type_it = table_loader.read_telescope_events_by_type_chunked(chunk_size)
+
+        iters = (event_it, tel_event_it, by_type_it)
 
-        for chunk, (events, tel_events) in enumerate(zip(event_it, tel_event_it)):
+        for chunk, (events, tel_events, by_type) in enumerate(zip(*iters)):
             n_read += len(events)
+            start = chunk * chunk_size
+            stop = min(n_events, (chunk + 1) * chunk_size)
 
             # last chunk might be smaller
             if chunk == (n_chunks - 1):
@@ -248,13 +329,17 @@ def test_chunked(dl2_shower_geometry_file):
             else:
                 assert len(events) == chunk_size
 
-            assert len(tel_events) == np.sum(events["tels_with_trigger"])
+            # check events are in compatible order
+            check_equal_array_event_order(events, tel_events)
+            check_equal_array_event_order(trigger[start:stop], events)
+
+            # check number of telescope events is correct
+            assert len(tel_events) == np.count_nonzero(events["tels_with_trigger"])
 
-            unique_events_from_tel_table = np.unique(tel_events[["obs_id", "event_id"]])
+            n_events_by_type = 0
+            for table in by_type.values():
+                n_events_by_type += len(table)
 
-            # 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_events_by_type == len(tel_events)
 
     assert n_read == n_events