Repository: cta-observatory/ctapipe
Pull request #1875 (2022-04-12 08:50:15)

Fix bug in ShowerProcessor / HillasReconstructor not using QualityQuery correctly, add support for HillasIntersection in ShowerProcessor
The code was requiring at least three! telescopes, while also the debug log output implied 2 and we should not hard code a minimum of 3 telescopes for stereo reconstruction.


Patch:
diff --git a/ctapipe/image/image_processor.py b/ctapipe/image/image_processor.py
index bf07ea89..cbffd1ca 100644
--- a/ctapipe/image/image_processor.py
+++ b/ctapipe/image/image_processor.py
@@ -12,7 +12,7 @@ from ..containers import (
     PeakTimeStatisticsContainer,
 )
 from ..core import QualityQuery, TelescopeComponent
-from ..core.traits import Bool, BoolTelescopeParameter, List, create_class_enum_trait
+from ..core.traits import Bool, List, create_class_enum_trait
 from ..instrument import SubarrayDescription
 from . import (
     ImageCleaner,
@@ -61,17 +61,11 @@ class ImageProcessor(TelescopeComponent):
     image_cleaner_type = create_class_enum_trait(
         base_class=ImageCleaner, default_value="TailcutsImageCleaner"
     )
-
     use_telescope_frame = Bool(
         default_value=True,
         help="Whether to calculate parameters in the telescope or camera frame",
     ).tag(config=True)
 
-    apply_image_modifier = BoolTelescopeParameter(
-        default_value=False,
-        help="If true, apply ImageModifier to dl1 images"
-    ).tag(config=True)
-
     def __init__(
         self, subarray: SubarrayDescription, config=None, parent=None, **kwargs
     ):
@@ -202,10 +196,8 @@ class ImageProcessor(TelescopeComponent):
                 geometry = self.telescope_frame_geometries[tel_id]
             else:
                 geometry = self.subarray.tel[tel_id].camera.geometry
-
-            if self.apply_image_modifier.tel[tel_id]:
-                dl1_camera.image = self.modify(tel_id=tel_id, image=dl1_camera.image)
-
+            # compute image parameters only if requested to write them
+            dl1_camera.image = self.modify(tel_id=tel_id, image=dl1_camera.image)
             dl1_camera.image_mask = self.clean(
                 tel_id=tel_id,
                 image=dl1_camera.image,
diff --git a/ctapipe/image/modifications.py b/ctapipe/image/modifications.py
index 5a6d4120..53e4b4b8 100644
--- a/ctapipe/image/modifications.py
+++ b/ctapipe/image/modifications.py
@@ -160,36 +160,30 @@ class ImageModifier(TelescopeComponent):
         self.rng = np.random.default_rng(self.rng_seed)
 
     def __call__(self, tel_id, image, rng=None):
-        dtype = image.dtype
-
-        if self.psf_smear_factor.tel[tel_id] > 0:
-            geom = self.subarray.tel[tel_id].camera.geometry
-            image = _smear_psf_randomly(
-                image=image,
-                fraction=self.psf_smear_factor.tel[tel_id],
-                indices=geom.neighbor_matrix_sparse.indices,
-                indptr=geom.neighbor_matrix_sparse.indptr,
-                smear_probabilities=np.full(geom.max_neighbors, 1 / geom.max_neighbors),
-                seed=self.rng.integers(0, np.iinfo(np.int64).max),
-            )
-
-        if (
-            self.noise_level_dim_pixels.tel[tel_id] > 0
-            or self.noise_level_bright_pixels.tel[tel_id] > 0
-        ):
-            bright_pixels = image > self.noise_transition_charge.tel[tel_id]
-            noise = np.where(
-                bright_pixels,
-                self.noise_level_bright_pixels.tel[tel_id],
-                self.noise_level_dim_pixels.tel[tel_id],
-            )
-            image = _add_noise(
-                image,
-                noise,
-                rng=self.rng,
-                correct_bias=self.noise_correct_bias.tel[tel_id],
-            )
-
-            image[~bright_pixels] += self.noise_bias_dim_pixels.tel[tel_id]
-
-        return image.astype(dtype, copy=False)
+        geom = self.subarray.tel[tel_id].camera.geometry
+        smeared_image = _smear_psf_randomly(
+            image=image,
+            fraction=self.psf_smear_factor.tel[tel_id],
+            indices=geom.neighbor_matrix_sparse.indices,
+            indptr=geom.neighbor_matrix_sparse.indptr,
+            smear_probabilities=np.full(geom.max_neighbors, 1 / geom.max_neighbors),
+            seed=self.rng.integers(0, np.iinfo(np.int64).max),
+        )
+
+        noise = np.where(
+            image > self.noise_transition_charge.tel[tel_id],
+            self.noise_level_bright_pixels.tel[tel_id],
+            self.noise_level_dim_pixels.tel[tel_id],
+        )
+        image_with_noise = _add_noise(
+            smeared_image,
+            noise,
+            rng=self.rng,
+            correct_bias=self.noise_correct_bias.tel[tel_id],
+        )
+        bias = np.where(
+            image > self.noise_transition_charge.tel[tel_id],
+            0,
+            self.noise_bias_dim_pixels.tel[tel_id],
+        )
+        return (image_with_noise + bias).astype(np.float32)
diff --git a/ctapipe/reco/__init__.py b/ctapipe/reco/__init__.py
index 68bbd011..d5c89e54 100644
--- a/ctapipe/reco/__init__.py
+++ b/ctapipe/reco/__init__.py
@@ -1,8 +1,12 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
-from .hillas_reconstructor import HillasReconstructor
+
+# reconstructors must be imported before ShowerProcessor, so
+# they are available there
 from .reco_algorithms import Reconstructor
-from .shower_processor import ShowerProcessor
+from .hillas_reconstructor import HillasReconstructor
 from .hillas_intersection import HillasIntersection
+
+from .shower_processor import ShowerProcessor
 from .impact import ImPACTReconstructor
 
 
diff --git a/ctapipe/reco/hillas_intersection.py b/ctapipe/reco/hillas_intersection.py
index b52d1382..613b0e48 100644
--- a/ctapipe/reco/hillas_intersection.py
+++ b/ctapipe/reco/hillas_intersection.py
@@ -23,7 +23,7 @@ from ctapipe.containers import (
 )
 from ctapipe.instrument import get_atmosphere_profile_functions
 
-from astropy.coordinates import SkyCoord
+from astropy.coordinates import SkyCoord, AltAz
 from ctapipe.coordinates import (
     NominalFrame,
     CameraFrame,
@@ -39,6 +39,9 @@ from ctapipe.core import traits
 __all__ = ["HillasIntersection"]
 
 
+INVALID = ReconstructedGeometryContainer(tel_ids=[])
+
+
 class HillasIntersection(Reconstructor):
     """
     This class is a simple re-implementation of Hillas parameter based event
@@ -69,11 +72,11 @@ class HillasIntersection(Reconstructor):
         ["Konrad", "hess"], default_value="Konrad", help="Weighting Method name"
     ).tag(config=True)
 
-    def __init__(self, config=None, parent=None, **kwargs):
+    def __init__(self, subarray, **kwargs):
         """
         Weighting must be a function similar to the weight_konrad already implemented
         """
-        super().__init__(config=config, parent=parent, **kwargs)
+        super().__init__(subarray, **kwargs)
 
         # We need a conversion function from height above ground to depth of maximum
         # To do this we need the conversion table from CORSIKA
@@ -84,7 +87,39 @@ class HillasIntersection(Reconstructor):
         if self.weighting == "Konrad":
             self._weight_method = self.weight_konrad
 
-    def predict(self, hillas_dict, subarray, array_pointing, telescopes_pointings=None):
+
+    def __call__(self, event):
+        '''
+        Perform stereo reconstruction on event.
+
+        Parameters
+        ----------
+        event: `ctapipe.containers.ArrayEventContainer`
+            The event, needs to have dl1 parameters
+
+        Returns
+        -------
+        ReconstructedGeometryContainer
+        '''
+        
+        try:
+            hillas_dict = self._create_hillas_dict(event)
+        except (TooFewTelescopesException, InvalidWidthException):
+            return INVALID
+
+        # Due to tracking the pointing of the array will never be a constant
+        array_pointing = SkyCoord(
+            az=event.pointing.array_azimuth,
+            alt=event.pointing.array_altitude,
+            frame=AltAz(),
+        )
+
+        telescope_pointings = self._get_telescope_pointings(event)
+
+        return self._predict(hillas_dict, array_pointing, telescope_pointings)
+
+
+    def _predict(self, hillas_dict, array_pointing, telescopes_pointings=None):
         """
 
         Parameters
@@ -115,12 +150,12 @@ class HillasIntersection(Reconstructor):
             )
 
         # check for np.nan or 0 width's as these screw up weights
-        if any([np.isnan(hillas_dict[tel]["width"].value) for tel in hillas_dict]):
+        if any([np.isnan(h.width.value) for h in hillas_dict.values()]):
             raise InvalidWidthException(
                 "A HillasContainer contains an ellipse of width==np.nan"
             )
 
-        if any([hillas_dict[tel]["width"].value == 0 for tel in hillas_dict]):
+        if any([h.width.value == 0 for h in hillas_dict.values()]):
             raise InvalidWidthException(
                 "A HillasContainer contains an ellipse of width==0"
             )
@@ -131,11 +166,11 @@ class HillasIntersection(Reconstructor):
             }
 
         tilted_frame = TiltedGroundFrame(pointing_direction=array_pointing)
-        grd_coord = subarray.tel_coords
+        grd_coord = self.subarray.tel_coords
         tilt_coord = grd_coord.transform_to(tilted_frame)
 
         tel_ids = list(hillas_dict.keys())
-        tel_indices = subarray.tel_ids_to_indices(tel_ids)
+        tel_indices = self.subarray.tel_ids_to_indices(tel_ids)
 
         tel_x = {
             tel_id: tilt_coord.x[tel_index]
@@ -152,7 +187,7 @@ class HillasIntersection(Reconstructor):
 
         for tel_id, hillas in hillas_dict.items():
             if isinstance(hillas, CameraHillasParametersContainer):
-                focal_length = subarray.tel[tel_id].optics.equivalent_focal_length
+                focal_length = self.subarray.tel[tel_id].optics.equivalent_focal_length
                 camera_frame = CameraFrame(
                     telescope_pointing=telescopes_pointings[tel_id],
                     focal_length=focal_length,
diff --git a/ctapipe/reco/hillas_reconstructor.py b/ctapipe/reco/hillas_reconstructor.py
index 39d5cf87..316b334c 100644
--- a/ctapipe/reco/hillas_reconstructor.py
+++ b/ctapipe/reco/hillas_reconstructor.py
@@ -37,6 +37,9 @@ from astropy import units as u
 __all__ = ["HillasPlane", "HillasReconstructor"]
 
 
+INVALID = ReconstructedGeometryContainer(tel_ids=[])
+
+
 def angle(v1, v2):
     """ computes the angle between two vectors
         assuming cartesian coordinates
@@ -162,10 +165,8 @@ class HillasReconstructor(Reconstructor):
 
     """
 
-    def __init__(self, subarray, config=None, parent=None, **kwargs):
-        super().__init__(config=config, parent=parent, **kwargs)
-        self.subarray = subarray
-
+    def __init__(self, subarray, **kwargs):
+        super().__init__(subarray=subarray, **kwargs)
         self._cam_radius_m = {
             cam: cam.geometry.guess_radius() for cam in subarray.camera_types
         }
@@ -183,16 +184,17 @@ class HillasReconstructor(Reconstructor):
 
         Parameters
         ----------
-        event : container
-            `ctapipe.containers.ArrayEventContainer`
-        """
+        event: `ctapipe.containers.ArrayEventContainer`
+            The event, needs to have dl1 parameters
 
-        # Read only valid HillasContainers
-        hillas_dict = {
-            tel_id: dl1.parameters.hillas
-            for tel_id, dl1 in event.dl1.tel.items()
-            if np.isfinite(dl1.parameters.hillas.intensity)
-        }
+        Returns
+        -------
+        ReconstructedGeometryContainer
+        """
+        try:
+            hillas_dict = self._create_hillas_dict(event)
+        except (TooFewTelescopesException, InvalidWidthException):
+            return INVALID
 
         # Due to tracking the pointing of the array will never be a constant
         array_pointing = SkyCoord(
@@ -210,17 +212,12 @@ class HillasReconstructor(Reconstructor):
             for tel_id in event.dl1.tel.keys()
         }
 
-        try:
-            result = self._predict(
-                event, hillas_dict, self.subarray, array_pointing, telescope_pointings
-            )
-        except (TooFewTelescopesException, InvalidWidthException):
-            result = ReconstructedGeometryContainer()
-
-        event.dl2.stereo.geometry["HillasReconstructor"] = result
+        return self._predict(
+            event, hillas_dict, array_pointing, telescope_pointings
+        )
 
     def _predict(
-        self, event, hillas_dict, subarray, array_pointing, telescopes_pointings
+        self, event, hillas_dict, array_pointing, telescopes_pointings
     ):
         """
         The function you want to call for the reconstruction of the
@@ -243,8 +240,6 @@ class HillasReconstructor(Reconstructor):
 
         Raises
         ------
-        TooFewTelescopesException
-            if len(hillas_dict) < 2
         InvalidWidthException
             if any width is np.nan or 0
         """
@@ -256,23 +251,8 @@ class HillasReconstructor(Reconstructor):
         # This should be substituted by a DL1 QualityQuery specific to this
         # reconstructor
 
-        # stereoscopy needs at least two telescopes
-        if len(hillas_dict) < 2:
-            raise TooFewTelescopesException(
-                "need at least two telescopes, have {}".format(len(hillas_dict))
-            )
-        # check for np.nan or 0 width's as these screw up weights
-        if any([np.isnan(hillas_dict[tel]["width"].value) for tel in hillas_dict]):
-            raise InvalidWidthException(
-                "A HillasContainer contains an ellipse of width==np.nan"
-            )
-        if any([hillas_dict[tel]["width"].value == 0 for tel in hillas_dict]):
-            raise InvalidWidthException(
-                "A HillasContainer contains an ellipse of width==0"
-            )
-
         hillas_planes, psi_core_dict = self.initialize_hillas_planes(
-            hillas_dict, subarray, telescopes_pointings, array_pointing
+            hillas_dict, self.subarray, telescopes_pointings, array_pointing
         )
 
         # algebraic direction estimate
diff --git a/ctapipe/reco/impact.py b/ctapipe/reco/impact.py
index 169e526f..8cc56768 100644
--- a/ctapipe/reco/impact.py
+++ b/ctapipe/reco/impact.py
@@ -12,6 +12,7 @@ from iminuit import Minuit
 from scipy.optimize import minimize, least_squares
 from scipy.stats import norm
 
+from ctapipe.core import Component
 from ctapipe.coordinates import (
     NominalFrame,
     TiltedGroundFrame,
@@ -24,7 +25,6 @@ from ctapipe.containers import (
     ReconstructedGeometryContainer,
     ReconstructedEnergyContainer,
 )
-from ctapipe.reco.reco_algorithms import Reconstructor
 from ctapipe.utils.template_network_interpolator import (
     TemplateNetworkInterpolator,
     TimeGradientInterpolator,
@@ -63,7 +63,7 @@ def xmax_prior(energy, xmax, width=100):
     return -2 * np.log(norm.pdf(diff / width))
 
 
-class ImPACTReconstructor(Reconstructor):
+class ImPACTReconstructor(Component):
     """This class is an implementation if the impact_reco Monte Carlo
     Template based image fitting method from parsons14.  This method uses a
     comparision of the predicted image from a library of image
diff --git a/ctapipe/reco/reco_algorithms.py b/ctapipe/reco/reco_algorithms.py
index ca8e3a85..bd6115c8 100644
--- a/ctapipe/reco/reco_algorithms.py
+++ b/ctapipe/reco/reco_algorithms.py
@@ -1,5 +1,11 @@
-from ctapipe.core import Component
-from ctapipe.containers import ReconstructedGeometryContainer
+from abc import abstractmethod
+import numpy as np
+
+from ctapipe.core import Component, QualityQuery
+from ctapipe.containers import ReconstructedGeometryContainer, ArrayEventContainer
+from astropy.coordinates import SkyCoord, AltAz
+
+from ctapipe.core.traits import List
 
 __all__ = ["Reconstructor", "TooFewTelescopesException", "InvalidWidthException"]
 
@@ -12,19 +18,33 @@ class InvalidWidthException(Exception):
     pass
 
 
+class StereoQualityQuery(QualityQuery):
+    """Quality criteria for dl1 parameters checked for telescope events to enter
+    into stereo reconstruction"""
+
+    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"),
+        ],
+        help=QualityQuery.quality_criteria.help,
+    ).tag(config=True)
+
+
 class Reconstructor(Component):
     """
     This is the base class from which all direction reconstruction
     algorithms should inherit from
     """
 
-    def __init__(self, *args, **kwargs):
-        """
-        Create a new instance of ImPACTReconstructor
-        """
-        super().__init__(*args, **kwargs)
+    def __init__(self, subarray, **kwargs):
+        super().__init__(**kwargs)
+        self.subarray = subarray
+        self.check_parameters = StereoQualityQuery(parent=self)
 
-    def predict(self, tels_dict):
+    @abstractmethod
+    def __call__(self, event: ArrayEventContainer):
         """overwrite this method with your favourite direction reconstruction
         algorithm
 
@@ -39,3 +59,37 @@ class Reconstructor(Component):
 
         """
         return ReconstructedGeometryContainer()
+
+    def _create_hillas_dict(self, event):
+        hillas_dict = {
+            tel_id: dl1.parameters.hillas
+            for tel_id, dl1 in event.dl1.tel.items()
+            if all(self.check_parameters(dl1.parameters))
+        }
+
+        if len(hillas_dict) < 2:
+            raise TooFewTelescopesException()
+
+        # check for np.nan or 0 width's as these screw up weights
+        if any([np.isnan(h.width.value) for h in hillas_dict.values()]):
+            raise InvalidWidthException(
+                "A HillasContainer contains an ellipse of width=np.nan"
+            )
+
+        if any([h.width.value == 0 for h in hillas_dict.values()]):
+            raise InvalidWidthException(
+                "A HillasContainer contains an ellipse of width=0"
+            )
+
+        return hillas_dict
+
+    @staticmethod
+    def _get_telescope_pointings(event):
+        return {
+            tel_id: SkyCoord(
+                alt=event.pointing.tel[tel_id].altitude,
+                az=event.pointing.tel[tel_id].azimuth,
+                frame=AltAz(),
+            )
+            for tel_id in event.dl1.tel.keys()
+        }
diff --git a/ctapipe/reco/shower_processor.py b/ctapipe/reco/shower_processor.py
index 63ca87b0..48f7af59 100644
--- a/ctapipe/reco/shower_processor.py
+++ b/ctapipe/reco/shower_processor.py
@@ -7,40 +7,26 @@ This processor will be able to process a shower/event in 3 steps:
 - estimation of classification (optional, currently unavailable)
 
 """
-from ctapipe.core import Component, QualityQuery
-from ctapipe.core.traits import List
-from ctapipe.containers import ArrayEventContainer, ReconstructedGeometryContainer
+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 HillasReconstructor
-
-
-DEFAULT_SHOWER_PARAMETERS = ReconstructedGeometryContainer(tel_ids=[])
-
-
-class ShowerQualityQuery(QualityQuery):
-    """Configuring shower-wise data checks."""
-
-    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"),
-        ],
-        help=QualityQuery.quality_criteria.help,
-    ).tag(config=True)
+from ctapipe.reco import Reconstructor
 
 
 class ShowerProcessor(Component):
     """
-    Needs DL1_PARAMETERS as input.
-    Should be run after ImageProcessor which produces such information.
+    Run the stereo event reconstruction on the input events.
 
-    For the moment it only supports the reconstruction of the shower geometry
-    using ctapipe.reco.HillasReconstructor.
+    This is mainly needed, so that the type of reconstructor can be chosen
+    via the configuration system.
 
-    It is planned to support also energy reconstruction and particle type
-    classification.
+    Input events must already contain dl1 parameters.
     """
+    reconstructor_type = create_class_enum_trait(
+        Reconstructor, default_value="HillasReconstructor",
+        help="The stereo geometry reconstructor to be used",
+    )
 
     def __init__(
         self, subarray: SubarrayDescription, config=None, parent=None, **kwargs
@@ -63,58 +49,10 @@ class ShowerProcessor(Component):
 
         super().__init__(config=config, parent=parent, **kwargs)
         self.subarray = subarray
-        self.check_shower = ShowerQualityQuery(parent=self)
-        self.reconstructor = HillasReconstructor(self.subarray)
-
-    def reconstruct_geometry(self, event, default=DEFAULT_SHOWER_PARAMETERS):
-        """Perform shower reconstruction.
-
-        Parameters
-        ----------
-        event : container
-            A ``ctapipe`` event container
-        default: container
-            The default 'ReconstructedGeometryContainer' which is
-            filled with NaNs.
-        Returns
-        -------
-        ReconstructedGeometryContainer:
-            direction in the sky with uncertainty,
-            core position on the ground with uncertainty,
-            h_max with uncertainty,
-            is_valid boolean for successfull reconstruction,
-            average intensity of the intensities used for reconstruction,
-            measure of algorithm success (if fit),
-            list of tel_ids used if stereo, or None if Mono
-        """
-
-        # Select only images which pass the shower quality criteria
-        hillas_dict = {
-            tel_id: dl1.parameters.hillas
-            for tel_id, dl1 in event.dl1.tel.items()
-            if all(self.check_shower(dl1.parameters))
-        }
-        self.log.debug("shower_criteria:\n %s", self.check_shower)
-
-        # Reconstruct the shower only if all shower criteria are met
-        if len(hillas_dict) > 2:
-
-            self.reconstructor(event)
-
-        else:
-            self.log.debug(
-                """Less than 2 images passed the quality cuts.
-                Returning default ReconstructedGeometryContainer container"""
-            )
-            event.dl2.stereo.geometry["HillasReconstructor"] = default
-
-    def process_shower_geometry(self, event: ArrayEventContainer):
-        """Record the reconstructed shower geometry into the ArrayEventContainer."""
-
-        self.reconstruct_geometry(event)
-
-        self.log.debug(
-            "shower geometry:\n %s", event.dl2.stereo.geometry["HillasReconstructor"]
+        self.reconstructor = Reconstructor.from_name(
+            self.reconstructor_type,
+            subarray=self.subarray,
+            parent=self,
         )
 
     def __call__(self, event: ArrayEventContainer):
@@ -129,6 +67,5 @@ class ShowerProcessor(Component):
         event : ctapipe.containers.ArrayEventContainer
             Top-level container for all event information.
         """
-
-        # This is always done when calling the ShowerProcessor
-        self.process_shower_geometry(event)
+        k = self.reconstructor_type
+        event.dl2.stereo.geometry[k] = self.reconstructor(event)
diff --git a/ctapipe/tools/process.py b/ctapipe/tools/process.py
index 83bae17b..b988c196 100644
--- a/ctapipe/tools/process.py
+++ b/ctapipe/tools/process.py
@@ -230,9 +230,9 @@ class ProcessorTool(Tool):
             )
 
         if self.should_compute_dl2:
-            shower_stats = self.process_shower.check_shower.to_table(functions=True)
+            reconstructor = self.process_shower.reconstructor
             write_table(
-                shower_stats,
+                reconstructor.check_parameters.to_table(functions=True),
                 self.write.output_path,
                 "/dl2/service/image_statistics",
                 append=True,
diff --git a/docs/tutorials/lst_analysis_bootcamp_2018.ipynb b/docs/tutorials/ctapipe_overview.ipynb
similarity index 72%
rename from docs/tutorials/lst_analysis_bootcamp_2018.ipynb
rename to docs/tutorials/ctapipe_overview.ipynb
index efb49afc..0d651ab7 100644
--- a/docs/tutorials/lst_analysis_bootcamp_2018.ipynb
+++ b/docs/tutorials/ctapipe_overview.ipynb
@@ -6,7 +6,7 @@
     "pycharm": {}
    },
    "source": [
-    "# 2018 LST Bootcamp walkthrough"
+    "# Analyzing Events Using ctapipe"
    ]
   },
   {
@@ -23,7 +23,7 @@
     "<img heoght=\"300px\" src=\"https://cta-observatory.github.io/ctapipe/_images/ctapipe_logo.png\" alt=\"ctapipe\"/>\n",
     "\n",
     "\n",
-    "<p style=\"text-align: center;\">LST Analysis Bootcamp</p>\n",
+    "<p style=\"text-align: center;\">Initially presented @ LST Analysis Bootcamp</p>\n",
     "\n",
     "<p style=\"text-align: center\">Padova, 26.11.2018</p>\n",
     "\n",
@@ -103,12 +103,13 @@
     "\n",
     "* DL0 → DL3 analysis\n",
     "\n",
-    "* Currently some R0 → DL0 code to be able to analyze simtel files\n",
+    "* Currently some R0 → DL2 code to be able to analyze simtel files\n",
     "\n",
     "* ctapipe is built upon the Scientific Python Stack, core dependencies are\n",
     "  * numpy\n",
     "  * scipy\n",
-    "  * astropy"
+    "  * astropy\n",
+    "  * numba"
    ]
   },
   {
@@ -125,13 +126,11 @@
     "* ctapipe is developed as Open Source Software (Currently under MIT License) at <https://github.com/cta-observatory/ctapipe>\n",
     "\n",
     "* We use the \"Github-Workflow\": \n",
-    "  * Few people (e.g. @kosack, @mackaiver) have write access to the main repository\n",
+    "  * Few people (e.g. @kosack, @maxnoe) have write access to the main repository\n",
     "  * Contributors fork the main repository and work on branches\n",
     "  * Pull Requests are merged after Code Review and automatic execution of the test suite\n",
     "\n",
-    "* Early developement stage ⇒ backwards-incompatible API changes might and will happen \n",
-    "\n",
-    "* Many open design questions ⇒ Core Developer Meeting in the second week of December in Dortmund"
+    "* Early developement stage ⇒ backwards-incompatible API changes might and will happen \n"
    ]
   },
   {
@@ -165,16 +164,9 @@
    "source": [
     "### What's still missing?\n",
     "\n",
-    "* Easy to use IO of analysis results to standard data formats (e.g. FITS, hdf5)\n",
-    "* Easy to use \"analysis builder\"\n",
-    "* A \"Standard Analysis\"\n",
     "* Good integration with machine learning techniques\n",
     "* IRF calculation \n",
-    "* Defining APIs for IO, instrument description access etc.\n",
-    "* Most code only tested on HESSIO simulations\n",
-    "* Documentation, e.g. formal definitions of coordinate frames \n",
-    " \n",
-    " "
+    "* Documentation, e.g. formal definitions of coordinate frames  "
    ]
   },
   {
@@ -634,7 +626,17 @@
     "pycharm": {}
    },
    "source": [
-    "### Putting it all together / Stereo reconstruction"
+    "### Putting it all together / Stereo reconstruction\n",
+    "\n",
+    "\n",
+    "All these steps are now unified in several components configurable through the config system, mainly:\n",
+    "\n",
+    "* CameraCalibrator for DL0 → DL1 (Images)\n",
+    "* ImageProcessor for DL1 (Images) → DL1 (Parameters)\n",
+    "* ShowerProcessor for stereo reconstruction of the shower geometry\n",
+    "* DataWriter for writing data into HDF5\n",
+    "\n",
+    "A command line tool doing these steps and writing out data in HDF5 format is available as `ctapipe-process`"
    ]
   },
   {
@@ -655,137 +657,72 @@
     "\n",
     "from ctapipe.calib import CameraCalibrator\n",
     "\n",
-    "from ctapipe.image import tailcuts_clean, number_of_islands\n",
-    "from ctapipe.image import hillas_parameters, leakage_parameters, concentration_parameters\n",
-    "from ctapipe.image import timing_parameters\n",
-    "\n",
-    "from ctapipe.reco import HillasReconstructor\n",
+    "from ctapipe.image import ImageProcessor\n",
+    "from ctapipe.reco import ShowerProcessor\n",
     "\n",
-    "from ctapipe.io import HDF5TableWriter\n",
+    "from ctapipe.io import DataWriter\n",
     "\n",
     "from copy import deepcopy\n",
     "import tempfile\n",
     "\n",
-    "# unoptimized cleaning levels, copied from \n",
-    "# https://github.com/tudo-astroparticlephysics/cta_preprocessing\n",
-    "cleaning_level = {\n",
-    "    'ASTRICam': (5, 7, 2),  # (5, 10)?\n",
-    "    'LSTCam': (3.5, 7.5, 2),  # ?? (3, 6) for Abelardo...\n",
-    "    'FlashCam': (4, 8, 2),  # there is some scaling missing?\n",
-    "}\n",
-    "\n",
+    "from traitlets.config import Config\n",
+    "\n",
+    "\n",
+    "image_processor_config = Config({\n",
+    "    \"ImageProcessor\": {\n",
+    "        \"image_cleaner_type\": \"TailcutsImageCleaner\",\n",
+    "        \"TailcutsImageCleaner\": {\n",
+    "            \"picture_threshold_pe\": [\n",
+    "                (\"type\", \"LST_LST_LSTCam\", 7.5),\n",
+    "                (\"type\", \"MST_MST_FlashCam\", 8),\n",
+    "                (\"type\", \"SST_ASTRI_ASTRICam\", 7),\n",
+    "            ],\n",
+    "            \"boundary_threshold_pe\": [\n",
+    "                (\"type\", \"LST_LST_LSTCam\", 5),\n",
+    "                (\"type\", \"MST_MST_FlashCam\", 4),\n",
+    "                (\"type\", \"SST_ASTRI_ASTRICam\", 4),\n",
+    "            ]\n",
+    "            \n",
+    "        }\n",
+    "    }\n",
+    "})\n",
     "\n",
     "input_url = get_dataset_path('gamma_test_large.simtel.gz')\n",
     "source = EventSource(input_url)\n",
     "\n",
     "calibrator = CameraCalibrator(subarray=source.subarray)\n",
-    "\n",
+    "image_processor = ImageProcessor(subarray=source.subarray, config=image_processor_config)\n",
+    "shower_processor = ShowerProcessor(subarray=source.subarray)\n",
     "horizon_frame = AltAz()\n",
     "\n",
     "f = tempfile.NamedTemporaryFile(suffix='.hdf5')\n",
     "\n",
-    "with HDF5TableWriter(f.name, mode='w', group_name='events') as writer:\n",
+    "with DataWriter(source, output_path=f.name, overwrite=True, write_stereo_shower=True) as writer:\n",
     "    \n",
     "    for event in source:\n",
-    "        print('Id: {}, E = {:1.3f}, Telescopes: {}'.format(event.count, event.simulation.shower.energy, len(event.r0.tel)))\n",
+    "        energy = event.simulation.shower.energy\n",
+    "        n_telescopes_r1 = len(event.r1.tel)\n",
+    "        event_id = event.index.event_id\n",
+    "        print(f'Id: {event_id}, E = {energy:1.3f}, Telescopes (R1): {n_telescopes_r1}')\n",
     "        \n",
     "        calibrator(event)\n",
-    "\n",
-    "        # mapping of telescope_id to parameters for stereo reconstruction\n",
-    "        hillas_containers = {}\n",
-    "        telescope_pointings = {}\n",
-    "        time_gradients = {}\n",
-    "\n",
-    "        for telescope_id, dl1 in event.dl1.tel.items():\n",
-    "            \n",
-    "            # Initialize an empty container for DL1b data\n",
-    "            event.dl1.tel[telescope_id].parameters = ImageParametersContainer()\n",
-    "            \n",
-    "            geometry = source.subarray.tels[telescope_id].camera.geometry\n",
-    "            image = dl1.image\n",
-    "            peak_time = dl1.peak_time\n",
-    "\n",
-    "            boundary, picture, min_neighbors = cleaning_level[geometry.camera_name]\n",
-    "\n",
-    "            clean = tailcuts_clean(\n",
-    "                geometry, \n",
-    "                image,\n",
-    "                boundary_thresh=boundary,\n",
-    "                picture_thresh=picture,\n",
-    "                min_number_picture_neighbors=min_neighbors\n",
-    "            )\n",
-    "\n",
-    "            # require more than five pixels after cleaning in each telescope\n",
-    "            if clean.sum() < 5:\n",
-    "                continue\n",
-    "\n",
-    "            hillas_c = hillas_parameters(geometry[clean], image[clean])\n",
-    "            event.dl1.tel[telescope_id].parameters.hillas = hillas_c\n",
-    "            \n",
-    "            leakage_c = leakage_parameters(geometry, image, clean)\n",
-    "            n_islands, island_ids = number_of_islands(geometry, clean)\n",
-    "\n",
-    "            # remove events with high leakage\n",
-    "            if leakage_c.intensity_width_2 > 0.2:\n",
-    "                continue\n",
-    "\n",
-    "            timing_c = timing_parameters(geometry, image, peak_time, hillas_c, clean)\n",
-    "\n",
-    "            hillas_containers[telescope_id] = hillas_c\n",
-    "\n",
-    "            # ssts have no timing in prod3b, so we'll use the skewness\n",
-    "            time_gradients[telescope_id] = timing_c.slope.value if geometry.camera_name != 'ASTRICam' else hillas_c.skewness\n",
-    "\n",
-    "            # this makes sure, that we get an arrow in the array plow for each telescope\n",
-    "            # might have the wrong direction though\n",
-    "            if abs(time_gradients[telescope_id]) < 0.2:\n",
-    "                time_gradients[telescope_id] = 1.0\n",
-    "\n",
-    "            telescope_pointings[telescope_id] = SkyCoord(\n",
-    "                alt=event.pointing.tel[telescope_id].altitude,\n",
-    "                az=event.pointing.tel[telescope_id].azimuth,\n",
-    "                frame=horizon_frame\n",
-    "            )  \n",
-    "            \n",
-    "        # the array pointing is needed for the creation of the TiltedFrame to perform the \n",
-    "        # impact point reconstruction\n",
-    "        array_pointing = SkyCoord(\n",
-    "            az=event.pointing.array_azimuth,\n",
-    "            alt=event.pointing.array_altitude,\n",
-    "            frame=horizon_frame\n",
-    "        )\n",
-    "\n",
-    "        if len(hillas_containers) > 2:\n",
-    "            reco = HillasReconstructor(source.subarray)\n",
-    "            stereo = reco._predict(\n",
-    "                event,\n",
-    "                hillas_containers,\n",
-    "                source.subarray,\n",
-    "                array_pointing,\n",
-    "                telescope_pointings\n",
-    "            )\n",
-    "\n",
-    "            writer.write('reconstructed', stereo)\n",
-    "            writer.write('true', event.simulation.shower)\n",
-    "    \n",
+    "        image_processor(event)\n",
+    "        shower_processor(event)\n",
+    "        \n",
+    "        stereo = event.dl2.stereo.geometry[\"HillasReconstructor\"]\n",
+    "        if stereo.is_valid:\n",
     "            print('  Alt: {:.2f}°'.format(stereo.alt.deg))\n",
     "            print('  Az: {:.2f}°'.format(stereo.az.deg))\n",
     "            print('  Hmax: {:.0f}'.format(stereo.h_max))\n",
     "            print('  CoreX: {:.1f}'.format(stereo.core_x))\n",
     "            print('  CoreY: {:.1f}'.format(stereo.core_y))\n",
+    "            print('  Multiplicity: {:d}'.format(len(stereo.tel_ids)))\n",
     "        \n",
     "        # save a nice event for plotting later\n",
     "        if event.count == 3:\n",
-    "            \n",
-    "            core_dict = {\n",
-    "                tel_id: dl1.parameters.core.psi for tel_id, dl1 in event.dl1.tel.items()\n",
-    "            }\n",
-    "            plotting_core = core_dict\n",
     "            plotting_event = deepcopy(event)\n",
-    "            plotting_hillas = hillas_containers\n",
-    "            plotting_timing = time_gradients\n",
-    "            plotting_stereo = stereo\n",
-    "    "
+    "            \n",
+    "        writer(event)"
    ]
   },
   {
@@ -799,17 +736,27 @@
     "from astropy.coordinates.angle_utilities import angular_separation\n",
     "import pandas as pd\n",
     "\n",
+    "from ctapipe.io import TableLoader\n",
     "\n",
-    "df_rec = pd.read_hdf(f.name, key='events/reconstructed')\n",
-    "df_true = pd.read_hdf(f.name, key='events/true')\n",
-    "\n",
+    "loader = TableLoader(f.name, load_dl2_geometry=True, load_simulated=True)\n",
     "\n",
+    "events = loader.read_subarray_events()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "pycharm": {}
+   },
+   "outputs": [],
+   "source": [
     "theta = angular_separation(\n",
-    "    df_rec.az.values * u.deg, df_rec.alt.values * u.deg,\n",
-    "    df_true.az.values * u.deg, df_true.alt.values * u.deg,\n",
+    "    events[\"HillasReconstructor_az\"].quantity, events[\"HillasReconstructor_alt\"].quantity,\n",
+    "    events[\"true_az\"].quantity, events[\"true_alt\"].quantity\n",
     ")\n",
     "\n",
-    "plt.hist(theta.to(u.deg).value**2, bins=25, range=[0, 0.3])\n",
+    "plt.hist(theta.to_value(u.deg)**2, bins=25, range=[0, 0.3])\n",
     "plt.xlabel(r'$\\theta² / deg²$')\n",
     "None"
    ]
@@ -836,6 +783,16 @@
     "\n",
     "angle_offset = plotting_event.pointing.array_azimuth\n",
     "\n",
+    "plotting_hillas = {\n",
+    "    tel_id: dl1.parameters.hillas\n",
+    "    for tel_id, dl1 in plotting_event.dl1.tel.items()\n",
+    "}\n",
+    "\n",
+    "plotting_core = {\n",
+    "    tel_id: dl1.parameters.core.psi\n",
+    "    for tel_id, dl1 in plotting_event.dl1.tel.items()\n",
+    "}\n",
+    "\n",
     "\n",
     "disp = ArrayDisplay(source.subarray)\n",
     "\n",
@@ -846,7 +803,8 @@
     "    s=200, c='k', marker='x', label='True Impact',\n",
     ")\n",
     "plt.scatter(\n",
-    "    plotting_stereo.core_x, plotting_stereo.core_y,\n",
+    "    plotting_event.dl2.stereo.geometry[\"HillasReconstructor\"].core_x,\n",
+    "    plotting_event.dl2.stereo.geometry[\"HillasReconstructor\"].core_y,\n",
     "    s=200, c='r', marker='x', label='Estimated Impact',\n",
     ")\n",
     "\n",
@@ -861,87 +819,18 @@
     "pycharm": {}
    },
    "source": [
-    "## LST Mono with output\n",
-    "\n",
-    "\n",
-    "* Let's use the `HDF5TableWriter` to save the dl1 Hillas parameter data to an hdf5 file\n",
-    "* This is not ideal yet and one of the major points to be discussed in two weeks"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "pycharm": {}
-   },
-   "outputs": [],
-   "source": [
-    "from ctapipe.io import HDF5TableWriter\n",
-    "from ctapipe.core.container import Container, Field\n",
-    "\n",
-    "input_url = get_dataset_path('gamma_test_large.simtel.gz')\n",
-    "\n",
-    "source = EventSource(\n",
-    "    input_url,\n",
-    "    allowed_tels=[1, 2, 3, 4], # only use the first LST\n",
-    ")\n",
-    "\n",
-    "calibrator = CameraCalibrator(subarray=source.subarray)\n",
-    "\n",
-    "class EventInfo(Container):\n",
-    "    event_id = Field('event_id')\n",
-    "    obs_id =  Field('obs_id')\n",
-    "    telescope_id = Field('telescope_id')\n",
-    "    \n",
-    "\n",
-    "with HDF5TableWriter(filename='hillas.h5', group_name='dl1', mode='w') as writer:\n",
-    "\n",
-    "    for event in source:\n",
-    "        print('Id: {}, E = {:1.3f}, Telescopes: {}'.format(event.count, event.simulation.shower.energy, len(event.r0.tel)))\n",
-    "    \n",
-    "        calibrator(event)\n",
-    "    \n",
-    "        for telescope_id, dl1 in event.dl1.tel.items():      \n",
-    "\n",
-    "            geometry = source.subarray.tels[telescope_id].camera.geometry\n",
-    "            image = dl1.image\n",
-    "            peak_time = dl1.peak_time\n",
-    "\n",
-    "            boundary, picture, min_neighbors = cleaning_level[geometry.camera_name]\n",
-    "\n",
-    "            clean = tailcuts_clean(\n",
-    "                geometry, \n",
-    "                image,\n",
-    "                boundary_thresh=boundary,\n",
-    "                picture_thresh=picture,\n",
-    "                min_number_picture_neighbors=min_neighbors\n",
-    "            )\n",
-    "            \n",
-    "            if clean.sum() < 5:\n",
-    "                continue\n",
-    "            \n",
-    "            event_info = EventInfo(event_id=event.index.event_id, obs_id=event.index.obs_id, telescope_id=telescope_id)\n",
-    "            hillas_c = hillas_parameters(geometry[clean], image[clean])\n",
-    "            leakage_c = leakage_parameters(geometry, image, clean)\n",
-    "            timing_c = timing_parameters(geometry, image, peak_time, hillas_c, clean)\n",
-    "\n",
-    "            writer.write('events', [event_info, event.simulation.shower, hillas_c, leakage_c, timing_c])\n",
-    "    \n"
+    "### Reading the LST dl1 data\n"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
-    "pycharm": {}
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
-    "import pandas as pd\n",
+    "loader = TableLoader(f.name, load_simulated=True, load_dl1_parameters=True)\n",
     "\n",
-    "df = pd.read_hdf('hillas.h5', key='dl1/events')\n",
-    "df.set_index(['obs_id', 'event_id', 'telescope_id'], inplace=True)\n",
-    "df.head()"
+    "dl1_table = loader.read_telescope_events([\"LST_LST_LSTCam\"])"
    ]
   },
   {
@@ -952,7 +841,10 @@
    },
    "outputs": [],
    "source": [
-    "plt.scatter(np.log10(df.energy), np.log10(df.intensity))\n",
+    "plt.scatter(\n",
+    "    np.log10(dl1_table[\"true_energy\"].quantity / u.TeV),\n",
+    "    np.log10(dl1_table[\"hillas_intensity\"]),\n",
+    ")\n",
     "plt.xlabel('log10(E / TeV)')\n",
     "plt.ylabel('log10(intensity)')\n",
     "None"
@@ -968,12 +860,12 @@
     "\n",
     "* Many of you might have heard: \"Python is slow\".\n",
     "* That's trueish.\n",
-    "* All python objects are classes living on the heap, event integers.\n",
+    "* All python objects are classes living on the heap, even integers.\n",
     "* Looping over lots of \"primitives\" is quite slow compared to other languages.\n",
     "\n",
     "⇒ Vectorize as much as possible using numpy  \n",
     "⇒ Use existing interfaces to fast C / C++ / Fortran code  \n",
-    "⇒ Optimize using cython or numba  \n",
+    "⇒ Optimize using numba  \n",
     "\n",
     "**But: \"Premature Optimization is the root of all evil\" — Donald Knuth**\n",
     "\n",
@@ -1225,11 +1117,27 @@
    "source": [
     "**A lot less code, and a factor 3 speed improvement**"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Finally, current ctapipe implementation is using numba:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%timeit number_of_islands(geometry, clean)"
+   ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -1243,7 +1151,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.7"
+   "version": "3.8.13"
   }
  },
  "nbformat": 4,
diff --git a/docs/tutorials/index.rst b/docs/tutorials/index.rst
index 475ff066..bfa4bb78 100644
--- a/docs/tutorials/index.rst
+++ b/docs/tutorials/index.rst
@@ -7,9 +7,8 @@ Tutorials
 .. toctree::
 
    ctapipe_handson
+   ctapipe_overview
    coordinates_example
    raw_data_exploration
    calibrated_data_exploration
    theta_square
-   lst_analysis_bootcamp_2018
-   
diff --git a/docs/tutorials/theta_square.ipynb b/docs/tutorials/theta_square.ipynb
index 74a1fd5e..b141c586 100644
--- a/docs/tutorials/theta_square.ipynb
+++ b/docs/tutorials/theta_square.ipynb
@@ -41,14 +41,13 @@
     "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
     "\n",
-    "from ctapipe.containers import ImageParametersContainer\n",
     "from ctapipe.io import EventSource\n",
-    "from ctapipe.visualization import CameraDisplay\n",
-    "from ctapipe.instrument import CameraGeometry\n",
+    "\n",
     "from ctapipe.calib import CameraCalibrator\n",
-    "from ctapipe.reco import HillasReconstructor\n",
-    "from ctapipe.image import hillas_parameters, tailcuts_clean\n",
-    "from ctapipe.utils import datasets"
+    "from ctapipe.image import ImageProcessor\n",
+    "from ctapipe.reco import ShowerProcessor\n",
+    "\n",
+    "from tqdm.auto import tqdm"
    ]
   },
   {
@@ -60,7 +59,7 @@
     }
    },
    "source": [
-    "Get source events in MC dataset. Here we stop at 10 events, just to make this example run fast, but for real use, one would need more statistics. "
+    "Get source events in MC dataset. "
    ]
   },
   {
@@ -74,12 +73,13 @@
    },
    "outputs": [],
    "source": [
-    "filename = datasets.get_dataset_path(\"gamma_test_large.simtel.gz\")\n",
-    "source = EventSource(filename, allowed_tels={1, 2, 3, 4}, max_events=10)\n",
+    "source = EventSource(\"dataset://gamma_test_large.simtel.gz\", allowed_tels={1, 2, 3, 4})\n",
+    "\n",
     "subarray = source.subarray\n",
-    "reco = HillasReconstructor(subarray=subarray)\n",
-    "calib = CameraCalibrator(subarray=source.subarray)\n",
-    "horizon_frame = AltAz()"
+    "\n",
+    "calib = CameraCalibrator(subarray=subarray)\n",
+    "image_processor = ImageProcessor(subarray=subarray)\n",
+    "shower_processor = ShowerProcessor(subarray=subarray)"
    ]
   },
   {
@@ -94,78 +94,20 @@
    "outputs": [],
    "source": [
     "off_angles = []\n",
-    "for event in source:\n",
+    "\n",
+    "for event in tqdm(source):\n",
     "\n",
     "    # calibrating the event\n",
     "    calib(event)\n",
-    "    \n",
-    "    hillas_params = {}\n",
-    "    \n",
-    "    # pointing direction of the telescopes\n",
-    "    telescope_pointings = {}\n",
-    "    \n",
-    "    subarray = source.subarray\n",
-    "    \n",
-    "    # get hillas params for each event in different telescopes\n",
-    "    for tel_id in event.dl0.tel.keys():\n",
-    "        \n",
-    "        event.dl1.tel[tel_id].parameters = ImageParametersContainer()\n",
-    "\n",
-    "        # telescope pointing direction\n",
-    "        telescope_pointings[tel_id] = SkyCoord(\n",
-    "            alt=event.pointing.tel[tel_id].altitude,\n",
-    "            az=event.pointing.tel[tel_id].azimuth,\n",
-    "            frame=horizon_frame\n",
-    "        )  \n",
-    "\n",
-    "        # Camera Geometry required for hillas parametrization\n",
-    "        camgeom = subarray.tel[tel_id].camera.geometry\n",
-    "\n",
-    "        # note the [0] is for channel 0 which is high-gain channel\n",
-    "        image = event.dl1.tel[tel_id].image\n",
-    "\n",
-    "        # Cleaning  of the image\n",
-    "        cleaned_image = image\n",
-    "        # create a clean mask of pixels above the threshold\n",
-    "        cleanmask = tailcuts_clean(\n",
-    "            camgeom, image, picture_thresh=10, boundary_thresh=5\n",
-    "        )\n",
-    "        # set all rejected pixels to zero\n",
-    "        cleaned_image[~cleanmask] = 0\n",
+    "    image_processor(event)\n",
+    "    shower_processor(event)\n",
     "\n",
-    "        # Calulate hillas parameters\n",
-    "        # It fails for empty pixels\n",
-    "        try:\n",
-    "            hillas = hillas_parameters(camgeom, cleaned_image)\n",
-    "        except Exception as e:\n",
-    "            print(e)\n",
-    "            pass\n",
-    "        \n",
-    "        if hillas.width.value > 0:\n",
-    "            hillas_params[tel_id] = hillas\n",
-    "\n",
-    "    if len(hillas_params) < 2:\n",
-    "        continue\n",
-    "\n",
-    "    array_pointing = SkyCoord(\n",
-    "        az=event.pointing.array_azimuth,\n",
-    "        alt=event.pointing.array_altitude,\n",
-    "        frame=horizon_frame\n",
-    "    )\n",
-    "    \n",
-    "    reco = HillasReconstructor(subarray=subarray)\n",
-    "    reco_result = reco._predict(\n",
-    "        event,\n",
-    "        hillas_params,\n",
-    "        source.subarray,\n",
-    "        array_pointing,\n",
-    "        telescope_pointings,\n",
-    "    )\n",
+    "    reco_result = event.dl2.stereo.geometry['HillasReconstructor']\n",
     "\n",
     "    # get angular offset between reconstructed shower direction and MC\n",
     "    # generated shower direction\n",
-    "    mcshower = event.simulation.shower\n",
-    "    off_angle = angular_separation(mcshower.az, mcshower.alt, reco_result.az, reco_result.alt)\n",
+    "    true_shower = event.simulation.shower\n",
+    "    off_angle = angular_separation(true_shower.az, true_shower.alt, reco_result.az, reco_result.alt)\n",
     "\n",
     "    # Appending all estimated off angles\n",
     "    off_angles.append(off_angle.to(u.deg).value)"
@@ -221,23 +163,16 @@
    },
    "outputs": [],
    "source": [
-    "plt.hist(thetasquare, bins=10, range=[0,0.4])\n",
+    "plt.hist(thetasquare, bins=10, range=[0, 0.4])\n",
     "plt.xlabel(r'$\\theta^2$ (deg)')\n",
     "plt.ylabel(\"# of events\")\n",
     "plt.show()"
    ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "again, this plot is not beautiful since we have such low stats"
-   ]
   }
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -251,7 +186,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.7"
+   "version": "3.8.13"
   },
   "toc": {
    "nav_menu": {


Test patch:
diff --git a/ctapipe/image/tests/test_modifications.py b/ctapipe/image/tests/test_modifications.py
index 27304942..f6a09673 100644
--- a/ctapipe/image/tests/test_modifications.py
+++ b/ctapipe/image/tests/test_modifications.py
@@ -70,16 +70,3 @@ def test_smear_image():
             # drawn from the poissonian (especially with signal value 1)
             if fraction > 0:
                 assert not ((image > 0) == (smeared > 0)).all()
-
-
-def test_defaults_no_change(example_subarray):
-    """Test that the default settings do not change the input image"""
-    rng = np.random.default_rng(0)
-
-    modifier = modifications.ImageModifier(example_subarray)
-    tel_id = 1
-    n_pixels = example_subarray.tel[tel_id].camera.geometry.n_pixels
-    image = rng.normal(50, 15, size=n_pixels).astype(np.float32)
-
-    new_image = modifier(tel_id=tel_id, image=image)
-    assert np.all(image == new_image)
diff --git a/ctapipe/reco/tests/test_HillasReconstructor.py b/ctapipe/reco/tests/test_HillasReconstructor.py
index 418d6461..7ee4f1ad 100644
--- a/ctapipe/reco/tests/test_HillasReconstructor.py
+++ b/ctapipe/reco/tests/test_HillasReconstructor.py
@@ -3,14 +3,13 @@ import numpy as np
 from astropy import units as u
 import pytest
 
-from ctapipe.containers import ImageParametersContainer, HillasParametersContainer
+from traitlets.config import Config
+from ctapipe.containers import HillasParametersContainer
+from ctapipe.image.image_processor import ImageProcessor
 from ctapipe.instrument import SubarrayDescription, TelescopeDescription
-from ctapipe.image.cleaning import tailcuts_clean
-from ctapipe.image.hillas import hillas_parameters, HillasParameterizationError
 from ctapipe.io import SimTelEventSource
 from ctapipe.reco.hillas_reconstructor import HillasReconstructor, HillasPlane
 from ctapipe.utils import get_dataset_path
-from ctapipe.coordinates import TelescopeFrame
 from astropy.coordinates import SkyCoord, AltAz
 from ctapipe.calib import CameraCalibrator
 
@@ -117,57 +116,46 @@ def test_invalid_events(subarray_and_event_gamma_off_axis_500_gev):
     # 4-LST bright event already calibrated
     # we'll clean it and parametrize it again in the TelescopeFrame
     subarray, event = subarray_and_event_gamma_off_axis_500_gev
-
-    tel_azimuth = {}
-    tel_altitude = {}
-
-    #source = EventSource(filename, max_events=1)
-    #subarray = source.subarray
     calib = CameraCalibrator(subarray)
-    fit = HillasReconstructor(subarray)
+    image_processor = ImageProcessor(subarray)
 
-    #for event in source:
+    # perform no quality cuts, so we can see if our additional checks on valid
+    # input work
+    config = Config({
+        "StereoQualityQuery": {
+            "quality_criteria": [],
+        }
+    })
+    hillas_reconstructor = HillasReconstructor(subarray, config=config)
 
     calib(event)
+    image_processor(event)
 
-    hillas_dict = {}
-    for tel_id, dl1 in event.dl1.tel.items():
-
-        geom = subarray.tel[tel_id].camera.geometry
-        tel_azimuth[tel_id] = event.pointing.tel[tel_id].azimuth
-        tel_altitude[tel_id] = event.pointing.tel[tel_id].altitude
-
-        mask = tailcuts_clean(
-            geom, dl1.image, picture_thresh=10.0, boundary_thresh=5.0
-        )
-
-        dl1.parameters = ImageParametersContainer()
-
-        try:
-            moments = hillas_parameters(geom[mask], dl1.image[mask])
-            hillas_dict[tel_id] = moments
-            dl1.parameters.hillas = moments
-        except HillasParameterizationError:
-            dl1.parameters.hillas = HillasParametersContainer()
-            continue
+    result = hillas_reconstructor(event)
+    assert result.is_valid
 
     # copy event container to modify it
-    event_copy = deepcopy(event)
+    invalid_event = deepcopy(event)
+
     # overwrite all image parameters but the last one with dummy ones
-    for tel_id in list(event_copy.dl1.tel.keys())[:-1]:
-        event_copy.dl1.tel[tel_id].parameters.hillas = HillasParametersContainer()
-    fit(event_copy)
-    assert event_copy.dl2.stereo.geometry["HillasReconstructor"].is_valid is False
+    for tel_id in list(invalid_event.dl1.tel.keys())[:-1]:
+        invalid_event.dl1.tel[tel_id].parameters.hillas = HillasParametersContainer()
 
+    result = hillas_reconstructor(invalid_event)
+    assert not result.is_valid
+
+    tel_id = list(invalid_event.dl1.tel.keys())[-1]
     # Now use the original event, but overwrite the last width to 0
-    event.dl1.tel[tel_id].parameters.hillas.width = 0 * u.m
-    fit(event)
-    assert event.dl2.stereo.geometry["HillasReconstructor"].is_valid is False
+    invalid_event = deepcopy(event)
+    invalid_event.dl1.tel[tel_id].parameters.hillas.width = 0 * u.m
+    result = hillas_reconstructor(invalid_event)
+    assert not result.is_valid
 
     # Now use the original event, but overwrite the last width to NaN
-    event.dl1.tel[tel_id].parameters.hillas.width = np.nan * u.m
-    fit(event)
-    assert event.dl2.stereo.geometry["HillasReconstructor"].is_valid is False
+    invalid_event = deepcopy(event)
+    invalid_event.dl1.tel[tel_id].parameters.hillas.width = np.nan * u.m
+    result = hillas_reconstructor(invalid_event)
+    assert not result.is_valid
 
 
 def test_reconstruction_against_simulation(subarray_and_event_gamma_off_axis_500_gev):
@@ -180,58 +168,14 @@ def test_reconstruction_against_simulation(subarray_and_event_gamma_off_axis_500
     subarray, event = subarray_and_event_gamma_off_axis_500_gev
 
     # define reconstructor
+    calib = CameraCalibrator(subarray)
+    image_processor = ImageProcessor(subarray)
     reconstructor = HillasReconstructor(subarray)
 
-    hillas_dict = {}
-    telescope_pointings = {}
-
-    for tel_id, dl1 in event.dl1.tel.items():
-
-        telescope_pointings[tel_id] = SkyCoord(
-            alt=event.pointing.tel[tel_id].altitude,
-            az=event.pointing.tel[tel_id].azimuth,
-            frame=AltAz(),
-        )
-
-        geom_CameraFrame = subarray.tel[tel_id].camera.geometry
-
-        # this could be done also out of this loop,
-        # but in case of real data each telescope would have a
-        # different telescope_pointing
-        geom_TelescopeFrame = geom_CameraFrame.transform_to(
-            TelescopeFrame(telescope_pointing=telescope_pointings[tel_id])
-        )
-
-        mask = tailcuts_clean(
-            geom_TelescopeFrame,
-            dl1.image,
-            picture_thresh=5.0,
-            boundary_thresh=2.5,
-            keep_isolated_pixels=False,
-            min_number_picture_neighbors=2,
-        )
-
-        try:
-            hillas_dict[tel_id] = hillas_parameters(
-                geom_TelescopeFrame[mask], dl1.image[mask]
-            )
-
-            # the original event is created from a
-            # pytest fixture with "session" scope, so it's always the same
-            # and if we used the same event we would overwrite the image
-            # parameters for the next tests, thus causing their failure
-            test_event = deepcopy(event)
-            test_event.dl1.tel[tel_id].parameters = ImageParametersContainer()
-            test_event.dl1.tel[tel_id].parameters.hillas = hillas_dict[tel_id]
-
-        except HillasParameterizationError as e:
-            print(e)
-            continue
-
     # Get shower geometry
-    reconstructor(event)
-    # get the result from the correct DL2 container
-    result = event.dl2.stereo.geometry["HillasReconstructor"]
+    calib(event)
+    image_processor(event)
+    result = reconstructor(event)
 
     # get the reconstructed coordinates in the sky
     reco_coord = SkyCoord(alt=result.alt, az=result.az, frame=AltAz())
@@ -244,9 +188,13 @@ def test_reconstruction_against_simulation(subarray_and_event_gamma_off_axis_500
     assert reco_coord.separation(true_coord) < 0.1 * u.deg
 
 
-@pytest.mark.parametrize("filename", 
-                         ["gamma_divergent_LaPalma_baseline_20Zd_180Az_prod3_test.simtel.gz",
-                         "gamma_LaPalma_baseline_20Zd_180Az_prod3b_test.simtel.gz"])
+@pytest.mark.parametrize(
+    "filename",
+    [
+        "gamma_divergent_LaPalma_baseline_20Zd_180Az_prod3_test.simtel.gz",
+        "gamma_LaPalma_baseline_20Zd_180Az_prod3b_test.simtel.gz"
+    ]
+)
 def test_CameraFrame_against_TelescopeFrame(filename):
 
     input_file = get_dataset_path(
@@ -255,88 +203,54 @@ def test_CameraFrame_against_TelescopeFrame(filename):
 
     source = SimTelEventSource(input_file, max_events=10)
 
+    # too few events survive for this test with the defautl quality criteria,
+    # use less restrictive ones
+    config = Config({
+        "StereoQualityQuery": {
+            "quality_criteria": [
+                ("valid_width", "lambda p: p.hillas.width.value > 0"),
+            ]
+        }
+    })
+
     calib = CameraCalibrator(subarray=source.subarray)
-    reconstructor = HillasReconstructor(source.subarray)
+    reconstructor = HillasReconstructor(source.subarray, config=config)
+    image_processor_camera_frame = ImageProcessor(
+        source.subarray, use_telescope_frame=False
+    )
+    image_processor_telescope_frame = ImageProcessor(
+        source.subarray, use_telescope_frame=True
+    )
 
     reconstructed_events = 0
 
-    for event in source:
+    for event_telescope_frame in source:
 
-        calib(event)
+        calib(event_telescope_frame)
         # make a copy of the calibrated event for the camera frame case
         # later we clean and paramretrize the 2 events in the same way
         # but in 2 different frames to check they return compatible results
-        event_camera_frame = deepcopy(event)
-
-        telescope_pointings = {}
-        hillas_dict_camera_frame = {}
-        hillas_dict_telescope_frame = {}
+        event_camera_frame = deepcopy(event_telescope_frame)
 
-        for tel_id, dl1 in event.dl1.tel.items():
+        image_processor_telescope_frame(event_telescope_frame)
+        image_processor_camera_frame(event_camera_frame)
 
-            event_camera_frame.dl1.tel[tel_id].parameters = ImageParametersContainer()
-            event.dl1.tel[tel_id].parameters = ImageParametersContainer()
+        result_camera_frame = reconstructor(event_camera_frame)
+        result_telescope_frame = reconstructor(event_telescope_frame)
 
-            # this is needed only here to transform the camera geometries
-            telescope_pointings[tel_id] = SkyCoord(
-                alt=event.pointing.tel[tel_id].altitude,
-                az=event.pointing.tel[tel_id].azimuth,
-                frame=AltAz(),
-            )
-
-            geom_camera_frame = source.subarray.tel[tel_id].camera.geometry
-
-            # this could be done also out of this loop,
-            # but in case of real data each telescope would have a
-            # different telescope_pointing
-            geom_telescope_frame = geom_camera_frame.transform_to(
-                TelescopeFrame(telescope_pointing=telescope_pointings[tel_id])
-            )
+        assert result_camera_frame.is_valid == result_telescope_frame.is_valid
 
-            mask = tailcuts_clean(
-                geom_telescope_frame, dl1.image, picture_thresh=10.0, boundary_thresh=5.0
-            )
-
-            try:
-                moments_camera_frame = hillas_parameters(
-                    geom_camera_frame[mask], dl1.image[mask]
-                )
-                moments_telescope_frame = hillas_parameters(
-                    geom_telescope_frame[mask], dl1.image[mask]
-                )
-
-                if (moments_camera_frame.width.value > 0) and (moments_telescope_frame.width.value > 0):
-                    event_camera_frame.dl1.tel[
-                        tel_id
-                    ].parameters.hillas = moments_camera_frame
-                    dl1.parameters.hillas = moments_telescope_frame
-
-                    hillas_dict_camera_frame[tel_id] = moments_camera_frame
-                    hillas_dict_telescope_frame[tel_id] = moments_telescope_frame
-                else:
-                    continue
-
-            except HillasParameterizationError as e:
-                print(e)
-                continue
-
-        if (len(hillas_dict_camera_frame) > 2) and (len(hillas_dict_telescope_frame) > 2):
-            reconstructor(event_camera_frame)
-            reconstructor(event)
+        if result_telescope_frame.is_valid:
             reconstructed_events += 1
-        else:  # this event was not good enough to be tested on
-            continue
-
-        # Compare old approach with new approach
-        result_camera_frame = event_camera_frame.dl2.stereo.geometry["HillasReconstructor"]
-        result_telescope_frame = event.dl2.stereo.geometry["HillasReconstructor"]
 
-        assert result_camera_frame.is_valid
-        assert result_telescope_frame.is_valid
+            for field, cam in result_camera_frame.items():
+                tel = getattr(result_telescope_frame, field)
 
-        for field in event.dl2.stereo.geometry["HillasReconstructor"].as_dict():
-            C = np.asarray(result_camera_frame.as_dict()[field])
-            T = np.asarray(result_telescope_frame.as_dict()[field])
-            assert (np.isclose(C, T, rtol=1e-03, atol=1e-03, equal_nan=True)).all()
+                if hasattr(cam, "unit"):
+                    assert u.isclose(cam, tel, rtol=1e-3, atol=1e-3 * tel.unit, equal_nan=True)
+                elif isinstance(cam, list):
+                    assert cam == tel
+                else:
+                    assert np.isclose(cam, tel, rtol=1e-3, atol=1e-3, equal_nan=True)
 
-    assert reconstructed_events > 0 # check that we reconstruct at least 1 event
\ No newline at end of file
+    assert reconstructed_events > 0 # check that we reconstruct at least 1 event
diff --git a/ctapipe/reco/tests/test_hillas_intersection.py b/ctapipe/reco/tests/test_hillas_intersection.py
index e05af2d3..f10394a3 100644
--- a/ctapipe/reco/tests/test_hillas_intersection.py
+++ b/ctapipe/reco/tests/test_hillas_intersection.py
@@ -3,24 +3,15 @@ import astropy.units as u
 from numpy.testing import assert_allclose
 import numpy as np
 from astropy.coordinates import SkyCoord
-from ctapipe.coordinates import NominalFrame, AltAz, CameraFrame, TelescopeFrame
+from ctapipe.coordinates import NominalFrame, AltAz, CameraFrame
 from ctapipe.containers import HillasParametersContainer
 
-from ctapipe.io import EventSource
-
-from ctapipe.utils import get_dataset_path
-
-from ctapipe.image.cleaning import tailcuts_clean
-from ctapipe.image.hillas import hillas_parameters, HillasParameterizationError
-from ctapipe.calib import CameraCalibrator
-
 
 def test_intersect():
     """
     Simple test to check the intersection of lines. Try to intersect positions at (0,1) and (1,0)
     with angles perpendicular and test they cross at (0,0)
     """
-    hill = HillasIntersection()
     x1 = 0
     y1 = 1
     theta1 = 90 * u.deg
@@ -29,7 +20,7 @@ def test_intersect():
     y2 = 0
     theta2 = 0 * u.deg
 
-    sx, sy = hill.intersect_lines(x1, y1, theta1, x2, y2, theta2)
+    sx, sy = HillasIntersection.intersect_lines(x1, y1, theta1, x2, y2, theta2)
 
     assert_allclose(sx, 0, atol=1e-6)
     assert_allclose(sy, 0, atol=1e-6)
@@ -40,7 +31,6 @@ def test_parallel():
     Simple test to check the intersection of lines. Try to intersect positions at (0,0) and (0,1)
     with angles parallel and check the behaviour
     """
-    hill = HillasIntersection()
     x1 = 0
     y1 = 0
     theta1 = 0 * u.deg
@@ -49,29 +39,29 @@ def test_parallel():
     y2 = 0
     theta2 = 0 * u.deg
 
-    sx, sy = hill.intersect_lines(x1, y1, theta1, x2, y2, theta2)
+    sx, sy = HillasIntersection.intersect_lines(x1, y1, theta1, x2, y2, theta2)
     assert_allclose(sx, np.nan, atol=1e-6)
     assert_allclose(sy, np.nan, atol=1e-6)
 
 
-def test_intersection_xmax_reco():
+def test_intersection_xmax_reco(example_subarray):
     """
     Test the reconstruction of xmax with two LSTs that are pointing at zenith = 0.
     The telescopes are places along the x and y axis at the same distance from the center.
     The impact point is hard-coded to be happening in the center of this cartesian system.
     """
-    hill_inter = HillasIntersection()
+    hill_inter = HillasIntersection(subarray=example_subarray)
 
-    horizon_frame = AltAz()
+    altaz = AltAz()
     zen_pointing = 10 * u.deg
 
     array_direction = SkyCoord(
-        alt=90 * u.deg - zen_pointing, az=0 * u.deg, frame=horizon_frame
+        alt=90 * u.deg - zen_pointing, az=0 * u.deg, frame=altaz
     )
     nom_frame = NominalFrame(origin=array_direction)
 
     source_sky_pos_reco = SkyCoord(
-        alt=90 * u.deg - zen_pointing, az=0 * u.deg, frame=horizon_frame
+        alt=90 * u.deg - zen_pointing, az=0 * u.deg, frame=altaz
     )
 
     nom_pos_reco = source_sky_pos_reco.transform_to(nom_frame)
@@ -106,13 +96,13 @@ def test_intersection_xmax_reco():
     print(x_max)
 
 
-def test_intersection_reco_impact_point_tilted():
+def test_intersection_reco_impact_point_tilted(example_subarray):
     """
     Function to test the reconstruction of the impact point in the tilted frame.
     This is done using a squared configuration, of which the impact point occupies a vertex,
     ad the three telescopes the other three vertices.
     """
-    hill_inter = HillasIntersection()
+    hill_inter = HillasIntersection(example_subarray)
 
     delta = 100 * u.m
     tel_x_dict = {1: delta, 2: -delta, 3: -delta}
@@ -132,11 +122,11 @@ def test_intersection_reco_impact_point_tilted():
     np.testing.assert_allclose(reco_konrad[1], -delta.to_value(u.m), atol=1e-8)
 
 
-def test_intersection_weighting_spoiled_parameters():
+def test_intersection_weighting_spoiled_parameters(example_subarray):
     """
     Test that the weighting scheme is useful especially when a telescope is 90 deg with respect to the other two
     """
-    hill_inter = HillasIntersection()
+    hill_inter = HillasIntersection(example_subarray)
 
     delta = 100 * u.m
     tel_x_dict = {1: delta, 2: -delta, 3: -delta}
@@ -157,13 +147,13 @@ def test_intersection_weighting_spoiled_parameters():
     np.testing.assert_allclose(reco_konrad_spoiled[1], -delta.to_value(u.m), atol=1e-1)
 
 
-def test_intersection_nominal_reconstruction():
+def test_intersection_nominal_reconstruction(example_subarray):
     """
     Testing the reconstruction of the position in the nominal frame with a three-telescopes system.
     This is done using a squared configuration, of which the impact point occupies a vertex,
     ad the three telescopes the other three vertices.
     """
-    hill_inter = HillasIntersection()
+    hill_inter = HillasIntersection(example_subarray)
 
     delta = 1.0 * u.m
     horizon_frame = AltAz()
@@ -227,108 +217,16 @@ def test_intersection_nominal_reconstruction():
     )
 
 
-def test_reconstruction():
-    """
-    a test of the complete fit procedure on one event including:
-    • tailcut cleaning
-    • hillas parametrisation
-    • direction fit
-    • position fit
-
-    in the end, proper units in the output are asserted """
-
-    filename = get_dataset_path("gamma_test_large.simtel.gz")
-
-    fit = HillasIntersection()
-
-    source = EventSource(filename, max_events=10)
-    calib = CameraCalibrator(source.subarray)
-
-    horizon_frame = AltAz()
-
-    reconstructed_events = 0
-
-    for event in source:
-        calib(event)
-
-        mc = event.simulation.shower
-        array_pointing = SkyCoord(az=mc.az, alt=mc.alt, frame=horizon_frame)
-
-        hillas_dict = {}
-        telescope_pointings = {}
-        telescope_frame = TelescopeFrame()
-        for tel_id, dl1 in event.dl1.tel.items():
-
-            geom = source.subarray.tel[tel_id].camera.geometry
-            geom_tel_frame = geom.transform_to(telescope_frame)
-
-            telescope_pointings[tel_id] = SkyCoord(
-                alt=event.pointing.tel[tel_id].altitude,
-                az=event.pointing.tel[tel_id].azimuth,
-                frame=horizon_frame,
-            )
-
-            mask = tailcuts_clean(
-                geom, dl1.image, picture_thresh=10.0, boundary_thresh=5.0
-            )
-
-            try:
-                moments = hillas_parameters(geom_tel_frame[mask], dl1.image[mask])
-                hillas_dict[tel_id] = moments
-            except HillasParameterizationError as e:
-                print(e)
-                continue
-
-        if len(hillas_dict) < 2:
-            continue
-        else:
-            reconstructed_events += 1
-
-        # divergent mode put to on even though the file has parallel pointing.
-        fit_result = fit.predict(
-            hillas_dict, source.subarray, array_pointing, telescope_pointings
-        )
-
-        print(fit_result)
-        print(event.simulation.shower.core_x, event.simulation.shower.core_y)
-        fit_result.alt.to(u.deg)
-        fit_result.az.to(u.deg)
-        fit_result.core_x.to(u.m)
-        assert fit_result.is_valid
-
-    assert reconstructed_events > 0
-
-
 def test_reconstruction_works(subarray_and_event_gamma_off_axis_500_gev):
     subarray, event = subarray_and_event_gamma_off_axis_500_gev
-    reconstructor = HillasIntersection()
-
-    array_pointing = SkyCoord(
-        az=event.pointing.array_azimuth,
-        alt=event.pointing.array_altitude,
-        frame=AltAz(),
-    )
-
-    hillas_dict = {
-        tel_id: dl1.parameters.hillas
-        for tel_id, dl1 in event.dl1.tel.items()
-        if dl1.parameters.hillas.width.value > 0
-    }
+    reconstructor = HillasIntersection(subarray)
 
-    telescope_pointings = {
-        tel_id: SkyCoord(alt=pointing.altitude, az=pointing.azimuth, frame=AltAz())
-        for tel_id, pointing in event.pointing.tel.items()
-        if tel_id in hillas_dict
-    }
     true_coord = SkyCoord(
         alt=event.simulation.shower.alt, az=event.simulation.shower.az, frame=AltAz()
     )
 
-    result = reconstructor.predict(
-        hillas_dict, subarray, array_pointing, telescope_pointings
-    )
+    result = reconstructor(event)
     reco_coord = SkyCoord(alt=result.alt, az=result.az, frame=AltAz())
-
     assert reco_coord.separation(true_coord) < 0.1 * u.deg
 
 
@@ -338,26 +236,13 @@ def test_selected_subarray(subarray_and_event_gamma_off_axis_500_gev):
 
     # remove telescopes 2 and 3 to see that HillasIntersection can work
     # with arbirary telescope ids
-    subarray = subarray.select_subarray([1, 4])
-
-    reconstructor = HillasIntersection()
-    array_pointing = SkyCoord(
-        az=event.pointing.array_azimuth,
-        alt=event.pointing.array_altitude,
-        frame=AltAz(),
-    )
+    allowed_tels = {1, 4}
+    for tel_id in subarray.tel.keys():
+        if tel_id not in allowed_tels:
+            event.dl1.tel.pop(tel_id, None)
 
-    # again, only use telescopes 1 and 4
-    hillas_dict = {
-        tel_id: dl1.parameters.hillas
-        for tel_id, dl1 in event.dl1.tel.items()
-        if dl1.parameters.hillas.width.value > 0 and tel_id in {1, 4}
-    }
-
-    telescope_pointings = {
-        tel_id: SkyCoord(alt=pointing.altitude, az=pointing.azimuth, frame=AltAz())
-        for tel_id, pointing in event.pointing.tel.items()
-        if tel_id in hillas_dict
-    }
+    subarray = subarray.select_subarray(allowed_tels)
 
-    reconstructor.predict(hillas_dict, subarray, array_pointing, telescope_pointings)
+    reconstructor = HillasIntersection(subarray)
+    result = reconstructor(event)
+    assert result.is_valid
diff --git a/ctapipe/reco/tests/test_reconstruction_methods.py b/ctapipe/reco/tests/test_reconstruction_methods.py
index d9063513..1989e396 100644
--- a/ctapipe/reco/tests/test_reconstruction_methods.py
+++ b/ctapipe/reco/tests/test_reconstruction_methods.py
@@ -1,16 +1,10 @@
 from astropy import units as u
-import numpy as np
 
-from ctapipe.containers import ImageParametersContainer, HillasParametersContainer
-from ctapipe.image.cleaning import tailcuts_clean
-from ctapipe.image.hillas import hillas_parameters, HillasParameterizationError
+from ctapipe.image import ImageProcessor
 from ctapipe.io import EventSource
 from ctapipe.reco import HillasReconstructor, HillasIntersection
 
-from ctapipe.reco.reco_algorithms import InvalidWidthException
-
 from ctapipe.utils import get_dataset_path
-from astropy.coordinates import SkyCoord, AltAz
 from ctapipe.calib import CameraCalibrator
 
 import pytest
@@ -39,77 +33,18 @@ def test_reconstructors(reconstructors):
     source = EventSource(filename, max_events=10)
     subarray = source.subarray
     calib = CameraCalibrator(source.subarray)
-    horizon_frame = AltAz()
+    image_processor = ImageProcessor(source.subarray)
 
     for event in source:
         calib(event)
-        sim_shower = event.simulation.shower
-        array_pointing = SkyCoord(
-            az=sim_shower.az, alt=sim_shower.alt, frame=horizon_frame
-        )
-
-        telescope_pointings = {}
-
-        for tel_id, dl1 in event.dl1.tel.items():
-
-            geom = source.subarray.tel[tel_id].camera.geometry
-
-            telescope_pointings[tel_id] = SkyCoord(
-                alt=event.pointing.tel[tel_id].altitude,
-                az=event.pointing.tel[tel_id].azimuth,
-                frame=horizon_frame,
-            )
-            mask = tailcuts_clean(
-                geom, dl1.image, picture_thresh=10.0, boundary_thresh=5.0
-            )
-
-            dl1.parameters = ImageParametersContainer()
-
-            try:
-                moments = hillas_parameters(geom[mask], dl1.image[mask])
-            except HillasParameterizationError:
-                dl1.parameters.hillas = HillasParametersContainer()
-                continue
-
-            # Make sure we provide only good images for the test
-            if np.isnan(moments.width.value) or (moments.width.value == 0):
-                dl1.parameters.hillas = HillasParametersContainer()
-            else:
-                dl1.parameters.hillas = moments
-
-        hillas_dict = {
-            tel_id: dl1.parameters.hillas
-            for tel_id, dl1 in event.dl1.tel.items()
-            if np.isfinite(dl1.parameters.hillas.intensity)
-        }
-        if len(hillas_dict) < 2:
-            continue
-
-        for count, reco_method in enumerate(reconstructors):
-            if reco_method is HillasReconstructor:
-                reconstructor = HillasReconstructor(subarray)
-
-                reconstructor(event)
-
-                event.dl2.stereo.geometry["HillasReconstructor"].alt.to(u.deg)
-                event.dl2.stereo.geometry["HillasReconstructor"].az.to(u.deg)
-                event.dl2.stereo.geometry["HillasReconstructor"].core_x.to(u.m)
-                assert event.dl2.stereo.geometry["HillasReconstructor"].is_valid
+        image_processor(event)
 
-            else:
-                reconstructor = reco_method()
+        for ReconstructorType in reconstructors:
+            reconstructor = ReconstructorType(subarray)
 
-                try:
-                    reconstructor_out = reconstructor.predict(
-                        hillas_dict,
-                        source.subarray,
-                        array_pointing,
-                        telescope_pointings,
-                    )
-                except InvalidWidthException:
-                    continue
+            reconstructor(event)
 
-                reconstructor_out.alt.to(u.deg)
-                reconstructor_out.az.to(u.deg)
-                reconstructor_out.core_x.to(u.m)
-                assert reconstructor_out.is_valid
+            name = ReconstructorType.__name__
+            assert event.dl2.stereo.geometry[name].alt.unit.is_equivalent(u.deg)
+            assert event.dl2.stereo.geometry[name].az.unit.is_equivalent(u.deg)
+            assert event.dl2.stereo.geometry[name].core_x.unit.is_equivalent(u.m)
diff --git a/ctapipe/reco/tests/test_shower_processor.py b/ctapipe/reco/tests/test_shower_processor.py
index f8444d3c..cd49fff7 100644
--- a/ctapipe/reco/tests/test_shower_processor.py
+++ b/ctapipe/reco/tests/test_shower_processor.py
@@ -2,6 +2,7 @@
 Tests for ShowerProcessor functionalities.
 """
 from numpy import isfinite
+import pytest
 
 from traitlets.config.loader import Config
 
@@ -10,7 +11,8 @@ from ctapipe.image import ImageProcessor
 from ctapipe.reco import ShowerProcessor
 
 
-def test_shower_processor_geometry(example_event, example_subarray):
+@pytest.mark.parametrize("reconstructor_type", ["HillasReconstructor", "HillasIntersection"])
+def test_shower_processor_geometry(example_event, example_subarray, reconstructor_type):
     """Ensure we get shower geometry when we input an event with parametrized images."""
 
     calibrate = CameraCalibrator(subarray=example_subarray)
@@ -21,15 +23,17 @@ def test_shower_processor_geometry(example_event, example_subarray):
         subarray=example_subarray, image_cleaner_type="MARSImageCleaner"
     )
 
-    process_shower = ShowerProcessor(subarray=example_subarray)
+    process_shower = ShowerProcessor(
+        subarray=example_subarray,
+        reconstructor_type=reconstructor_type
+    )
 
     calibrate(example_event)
     process_images(example_event)
 
     process_shower(example_event)
-    print(process_shower.check_shower.to_table())
 
-    DL2a = example_event.dl2.stereo.geometry["HillasReconstructor"]
+    DL2a = example_event.dl2.stereo.geometry[reconstructor_type]
     print(DL2a)
     assert isfinite(DL2a.alt)
     assert isfinite(DL2a.az)
@@ -40,16 +44,19 @@ def test_shower_processor_geometry(example_event, example_subarray):
     assert isfinite(DL2a.average_intensity)
 
     # Increase some quality cuts and check that we get defaults
-    config.ShowerQualityQuery.quality_criteria = [
+    config[reconstructor_type].StereoQualityQuery.quality_criteria = [
         ("> 500 phes", "lambda p: p.hillas.intensity > 500")
     ]
 
-    process_shower = ShowerProcessor(config=config, subarray=example_subarray)
+    process_shower = ShowerProcessor(
+        config=config,
+        subarray=example_subarray,
+        reconstructor_type=reconstructor_type,
+    )
 
     process_shower(example_event)
-    print(process_shower.check_shower.to_table())
 
-    DL2a = example_event.dl2.stereo.geometry["HillasReconstructor"]
+    DL2a = example_event.dl2.stereo.geometry[reconstructor_type]
     print(DL2a)
     assert not isfinite(DL2a.alt)
     assert not isfinite(DL2a.az)
diff --git a/ctapipe/tools/tests/resources/image_modification_config.json b/ctapipe/tools/tests/resources/image_modification_config.json
index 7456d055..a6b657f3 100644
--- a/ctapipe/tools/tests/resources/image_modification_config.json
+++ b/ctapipe/tools/tests/resources/image_modification_config.json
@@ -18,7 +18,6 @@
     },
     "ImageProcessor": {
         "image_cleaner_type": "TailcutsImageCleaner",
-        "apply_image_modifier": [["type", "*", true]],
         "TailcutsImageCleaner": {
             "picture_threshold_pe": [
                 ["type", "*", 10.0],
@@ -62,12 +61,15 @@
         }
     },
     "ShowerProcessor": {
-        "ShowerQualityQuery": {
-            "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"]
-            ]
+        "reconstructor_type": "HillasReconstructor",
+        "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"]
+                ]
+            }
         }
     }
 }
diff --git a/ctapipe/tools/tests/resources/stage1_config.json b/ctapipe/tools/tests/resources/stage1_config.json
index 6ee328a8..2fbef2b7 100644
--- a/ctapipe/tools/tests/resources/stage1_config.json
+++ b/ctapipe/tools/tests/resources/stage1_config.json
@@ -41,12 +41,15 @@
         }
     },
     "ShowerProcessor": {
-        "ShowerQualityQuery": {
-            "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"]
-            ]
+        "reconstructor_type": "HillasReconstructor",
+        "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"]
+                ]
+            }
         }
     }
 }
diff --git a/ctapipe/tools/tests/resources/stage2_config.json b/ctapipe/tools/tests/resources/stage2_config.json
index adece5cf..586faa64 100644
--- a/ctapipe/tools/tests/resources/stage2_config.json
+++ b/ctapipe/tools/tests/resources/stage2_config.json
@@ -41,13 +41,15 @@
         }
     },
     "ShowerProcessor": {
-        "ShowerQualityQuery": {
-            "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 < 2"]
-            ]
+        "reconstructor_type": "HillasReconstructor",
+        "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"]
+                ]
+            }
         }
     }
 }
diff --git a/ctapipe/tools/tests/resources/training_config.json b/ctapipe/tools/tests/resources/training_config.json
index 0c3decb4..35120dcb 100644
--- a/ctapipe/tools/tests/resources/training_config.json
+++ b/ctapipe/tools/tests/resources/training_config.json
@@ -41,12 +41,15 @@
         }
     },
     "ShowerProcessor": {
-        "ShowerQualityQuery": {
-            "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"]
-            ]
+        "reconstructor_type": "HillasReconstructor",
+        "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"]
+                ]
+            }
         }
     }
 }