Repository: cta-observatory/ctapipe
Pull request #2168 (2022-12-15 11:54:39)

Make muon analysis part of the processing tool

RuneDominik: Use this instead of #1618 
RuneDominik: As of now, remaining tasks for this PR are (alongside the comments above ofc):

- [ ] Add/Restructure tests
- [ ] Remove old muon analysis tool
- [ ] Add docstrings/documentation and adapt config files  

Patch:
diff --git a/ctapipe/containers.py b/ctapipe/containers.py
index ebf56d83..bcccaef7 100644
--- a/ctapipe/containers.py
+++ b/ctapipe/containers.py
@@ -961,10 +961,14 @@ class EventCalibrationContainer(Container):
 
 
 class MuonRingContainer(Container):
-    """Container for the result of a ring fit, center_x, center_y"""
+    """Container for the result of a ring fit in telescope frame"""
 
-    center_x = Field(nan * u.deg, "center (x) of the fitted muon ring", unit=u.deg)
-    center_y = Field(nan * u.deg, "center (y) of the fitted muon ring", unit=u.deg)
+    center_fov_lon = Field(
+        nan * u.deg, "center (fov_lon) of the fitted muon ring", unit=u.deg
+    )
+    center_fov_lat = Field(
+        nan * u.deg, "center (fov_lat) of the fitted muon ring", unit=u.deg
+    )
     radius = Field(nan * u.deg, "radius of the fitted muon ring", unit=u.deg)
     center_phi = Field(
         nan * u.deg, "Angle of ring center within camera plane", unit=u.deg
@@ -992,7 +996,31 @@ class MuonParametersContainer(Container):
     )
     intensity_ratio = Field(nan, "Intensity ratio of pixels in the ring to all pixels")
     mean_squared_error = Field(
-        nan, "MSE of the deviation of all pixels after cleaning from the ring fit"
+        nan * u.deg**2,
+        "MSE of the deviation of all pixels after cleaning from the ring fit",
+    )
+
+
+class MuonTelescopeContainer(Container):
+    """
+    Container for muon analysis
+    """
+
+    ring = Field(default_factory=MuonRingContainer, description="muon ring fit")
+    parameters = Field(
+        default_factory=MuonParametersContainer, description="muon parameters"
+    )
+    efficiency = Field(
+        default_factory=MuonEfficiencyContainer, description="muon efficiency"
+    )
+
+
+class MuonContainer(Container):
+    """Root container for muon parameters"""
+
+    tel = Field(
+        default_factory=partial(Map, MuonTelescopeContainer),
+        description="map of tel_id to MuonTelescopeContainer",
     )
 
 
@@ -1232,6 +1260,9 @@ class ArrayEventContainer(Container):
         default_factory=MonitoringContainer,
         description="container for event-wise monitoring data (MON)",
     )
+    muon = Field(
+        default_factory=MuonContainer, description="Container for muon analysis results"
+    )
 
 
 class SchedulingBlockContainer(Container):
diff --git a/ctapipe/image/muon/__init__.py b/ctapipe/image/muon/__init__.py
index e0f7058e..f7025eb1 100644
--- a/ctapipe/image/muon/__init__.py
+++ b/ctapipe/image/muon/__init__.py
@@ -1,14 +1,13 @@
-from .fitting import kundu_chaudhuri_circle_fit
 from .features import (
-    mean_squared_error,
     intensity_ratio_inside_ring,
+    mean_squared_error,
     ring_completeness,
     ring_containment,
 )
-
-from .ring_fitter import MuonRingFitter
+from .fitting import kundu_chaudhuri_circle_fit
 from .intensity_fitter import MuonIntensityFitter
-
+from .processor import MuonProcessor
+from .ring_fitter import MuonRingFitter
 
 __all__ = [
     "MuonIntensityFitter",
@@ -18,4 +17,5 @@ __all__ = [
     "intensity_ratio_inside_ring",
     "ring_completeness",
     "ring_containment",
+    "MuonProcessor",
 ]
diff --git a/ctapipe/image/muon/processor.py b/ctapipe/image/muon/processor.py
new file mode 100644
index 00000000..30b382ae
--- /dev/null
+++ b/ctapipe/image/muon/processor.py
@@ -0,0 +1,276 @@
+"""
+High level muon analysis  (MuonProcessor Component)
+"""
+import numpy as np
+
+from ctapipe.containers import (
+    ArrayEventContainer,
+    MuonParametersContainer,
+    MuonTelescopeContainer,
+)
+from ctapipe.coordinates import TelescopeFrame
+from ctapipe.core import QualityQuery, TelescopeComponent
+from ctapipe.core.traits import FloatTelescopeParameter, List, Tuple, Unicode
+
+from .features import (
+    intensity_ratio_inside_ring,
+    mean_squared_error,
+    ring_completeness,
+    ring_containment,
+)
+from .intensity_fitter import MuonIntensityFitter
+from .ring_fitter import MuonRingFitter
+
+INVALID = MuonTelescopeContainer()
+INVALID_PARAMETERS = MuonParametersContainer()
+
+__all__ = ["MuonProcessor"]
+
+
+class ImageParameterQuery(QualityQuery):
+    """
+    For configuring quality checks performed on image parameters
+    before rings are fitted.
+    """
+
+    quality_criteria = List(
+        Tuple(Unicode(), Unicode()),
+        default_value=[
+            ("min_pixels", "dl1_params.morphology.n_pixels > 100"),
+            ("min_intensity", "dl1_params.hillas.intensity > 500"),
+        ],
+        help=(
+            "Quality cuts as list of tuples of ('description', 'expression string'), "
+            "e.g. ``[('min_pixels', 'dl1_params.morphology.n_pixels > 100'),]``, "
+            "to select muon images for analysis. "
+            "You may use ``numpy`` as ``np`` and ``astropy.units`` as ``u``, "
+            "but no other modules. "
+            "DL1 image parameters can be accessed by prefixing the wanted parameter "
+            "with ``dl1_params``. "
+        ),
+    ).tag(config=True)
+
+
+class RingQuery(QualityQuery):
+    """
+    For configuring quality checks performed on the extracted rings before
+    computing efficiency parameters.
+    """
+
+    quality_criteria = List(
+        Tuple(Unicode(), Unicode()),
+        default_value=[
+            ("radius_not_nan", "np.isfinite(ring.radius.value)"),
+            ("min_pixels", "np.count_nonzero(mask) > 50"),
+            ("ring_containment", "parameters.containment > 0.5"),
+        ],
+        help=(
+            "Quality cuts as list of tuples of ('description', 'expression string'), "
+            "e.g. ``[('radius_not_nan', 'np.isfinite(ring.radius.value)'),]``, "
+            "to select fitted muons for further intensity fitting. "
+            "You may use ``numpy`` as ``np`` and ``astropy.units`` as ``u``, but "
+            "no other modules. "
+            "Ring parameters and geometry can be accessed by prefixing the wanted "
+            "parameter with ``parameters`` or ``ring``, the ring mask can be accessed "
+            "as ``mask``."
+        ),
+    ).tag(config=True)
+
+
+class MuonProcessor(TelescopeComponent):
+    """
+    Takes cleaned images and extracts muon rings. Should be run after ImageProcessor.
+    """
+
+    completeness_threshold = FloatTelescopeParameter(
+        default_value=30.0, help="Threshold for calculating the ``ring_completeness``"
+    ).tag(config=True)
+
+    ratio_width = FloatTelescopeParameter(
+        default_value=1.5,
+        help="Ring width for intensity ratio computation in units of pixel diameter",
+    ).tag(config=True)
+
+    pedestal = FloatTelescopeParameter(
+        help="Pedestal noise rms", default_value=1.1
+    ).tag(config=True)
+
+    def __init__(self, subarray, **kwargs):
+        super().__init__(subarray, **kwargs)
+        self.dl1_query = ImageParameterQuery(parent=self)
+        self.ring_query = RingQuery(parent=self)
+
+        frame = TelescopeFrame()
+        self.geometries = {
+            tel_id: tel.camera.geometry.transform_to(frame)
+            for tel_id, tel in self.subarray.tel.items()
+        }
+        self.fov_radius = {
+            tel_id: geometry.guess_radius()
+            for tel_id, geometry in self.geometries.items()
+        }
+
+        self.ring_fitter = MuonRingFitter(parent=self)
+
+        self.intensity_fitter = MuonIntensityFitter(subarray=subarray, parent=self)
+
+    def __call__(self, event: ArrayEventContainer):
+        for tel_id in event.dl1.tel:
+            self._process_telescope_event(event, tel_id)
+
+    def _process_telescope_event(self, event, tel_id):
+        """
+        Extract and process a ring from a single image.
+
+        Parameters
+        ----------
+        event: ArrayEventContainer
+            Collection of all event information
+        tel_id: int
+            Telescope ID of the instrument that has measured the image
+        """
+        event_index = event.index
+        event_id = event_index.event_id
+
+        if self.subarray.tel[tel_id].optics.n_mirrors != 1:
+            self.log.warning(
+                f"Skipping non-single mirror telescope {tel_id},"
+                " muon-fitting for dual-mirror telescope is currently"
+                " not supported. Exclude dual mirror telescopes via setting"
+                " 'EventSource.allowed_tels'."
+            )
+            event.muon.tel[tel_id] = INVALID
+            return
+
+        self.log.debug(f"Processing event {event_id}, telescope {tel_id}")
+        dl1 = event.dl1.tel[tel_id]
+        image = dl1.image
+        mask = dl1.image_mask
+        if mask is None:
+            mask = image > 0
+
+        checks = self.dl1_query(dl1_params=dl1.parameters)
+
+        if not all(checks):
+            event.muon.tel[tel_id] = INVALID
+            return
+
+        geometry = self.geometries[tel_id]
+        fov_lon = geometry.pix_x
+        fov_lat = geometry.pix_y
+
+        # iterative ring fit.
+        # First use cleaning pixels, then only pixels close to the ring
+        # three iterations seems to be enough for most rings
+        for _ in range(3):
+            ring = self.ring_fitter(fov_lon, fov_lat, image, mask)
+            dist = np.sqrt(
+                (fov_lon - ring.center_fov_lon) ** 2
+                + (fov_lat - ring.center_fov_lat) ** 2
+            )
+            mask = np.abs(dist - ring.radius) / ring.radius < 0.4
+
+        parameters = self._calculate_muon_parameters(
+            tel_id, image, dl1.image_mask, ring
+        )
+
+        checks = self.ring_query(parameters=parameters, ring=ring, mask=mask)
+        if not all(checks):
+            event.muon.tel[tel_id] = MuonTelescopeContainer(
+                parameters=parameters, ring=ring
+            )
+            return
+
+        efficiency = self.intensity_fitter(
+            tel_id,
+            ring.center_fov_lon,
+            ring.center_fov_lat,
+            ring.radius,
+            image,
+            mask=mask,
+            pedestal=np.full(mask.shape, self.pedestal.tel[tel_id]),
+        )
+
+        self.log.debug(
+            f"Muon fit: r={ring.radius:.2f}"
+            f", width={efficiency.width:.4f}"
+            f", efficiency={efficiency.optical_efficiency:.2%}"
+        )
+
+        event.muon.tel[tel_id] = MuonTelescopeContainer(
+            ring=ring, efficiency=efficiency, parameters=parameters
+        )
+
+    def _calculate_muon_parameters(
+        self, tel_id, image, clean_mask, ring
+    ) -> MuonParametersContainer:
+        """
+        Calculate features from identified muon rings.
+
+        Parameters
+        ----------
+        tel_id: int
+            Telescope ID of the instrument that has measured the image
+        image: np.ndarray
+            Image to process
+        clean_mask: np.ndarray[bool]
+            DL1 Image cleaning mask
+        ring: MuonRingContainer
+            Collection of the fitted rings parameters
+
+        Returns
+        -------
+        MuonParametersContainer:
+            Collection of the fitted rings containment in the camera,
+            completeness, intensity ratio and the pixels MSE around
+            the fitted ring.
+        """
+        if np.isnan(ring.radius.value):
+            return INVALID_PARAMETERS
+
+        geometry = self.geometries[tel_id]
+        fov_radius = self.fov_radius[tel_id]
+        fov_lon = geometry.pix_x
+        fov_lat = geometry.pix_y
+
+        # add ring containment, not filled in fit
+        containment = ring_containment(
+            ring.radius, ring.center_fov_lon, ring.center_fov_lat, fov_radius
+        )
+
+        completeness = ring_completeness(
+            fov_lon,
+            fov_lat,
+            image,
+            ring.radius,
+            ring.center_fov_lon,
+            ring.center_fov_lat,
+            threshold=self.completeness_threshold.tel[tel_id],
+        )
+
+        pixel_width = geometry.pixel_width[clean_mask]
+        intensity_ratio = intensity_ratio_inside_ring(
+            fov_lon[clean_mask],
+            fov_lat[clean_mask],
+            image[clean_mask],
+            ring.radius,
+            ring.center_fov_lon,
+            ring.center_fov_lat,
+            width=self.ratio_width.tel[tel_id] * pixel_width,
+        )
+
+        mse = mean_squared_error(
+            fov_lon[clean_mask],
+            fov_lat[clean_mask],
+            image[clean_mask],
+            ring.radius,
+            ring.center_fov_lon,
+            ring.center_fov_lat,
+        )
+
+        return MuonParametersContainer(
+            containment=containment,
+            completeness=completeness,
+            intensity_ratio=intensity_ratio,
+            mean_squared_error=mse,
+        )
diff --git a/ctapipe/image/muon/ring_fitter.py b/ctapipe/image/muon/ring_fitter.py
index 5483a865..8f8d49ee 100644
--- a/ctapipe/image/muon/ring_fitter.py
+++ b/ctapipe/image/muon/ring_fitter.py
@@ -1,9 +1,10 @@
 import numpy as np
-from ctapipe.core import Component
-from ctapipe.containers import MuonRingContainer
-from .fitting import kundu_chaudhuri_circle_fit, taubin_circle_fit
 import traitlets as traits
 
+from ctapipe.containers import MuonRingContainer
+from ctapipe.core import Component
+
+from .fitting import kundu_chaudhuri_circle_fit, taubin_circle_fit
 
 # the fit methods do not expose the same interface, so we
 # force the same interface onto them, here.
@@ -11,14 +12,14 @@ import traitlets as traits
 # exposed to the user via the string traitlet `fit_method`
 
 
-def kundu_chaudhuri(x, y, weights, mask):
-    """kundu_chaudhuri_circle_fit with x, y, weights, mask interface"""
-    return kundu_chaudhuri_circle_fit(x[mask], y[mask], weights[mask])
+def kundu_chaudhuri(fov_lon, fov_lat, weights, mask):
+    """kundu_chaudhuri_circle_fit with fov_lon, fov_lat, weights, mask interface"""
+    return kundu_chaudhuri_circle_fit(fov_lon[mask], fov_lat[mask], weights[mask])
 
 
-def taubin(x, y, weights, mask):
-    """taubin_circle_fit with x, y, weights, mask interface"""
-    return taubin_circle_fit(x, y, mask)
+def taubin(fov_lon, fov_lat, weights, mask):
+    """taubin_circle_fit with fov_lon, fov_lat, weights, mask interface"""
+    return taubin_circle_fit(fov_lon, fov_lat, mask)
 
 
 FIT_METHOD_BY_NAME = {m.__name__: m for m in [kundu_chaudhuri, taubin]}
@@ -34,17 +35,19 @@ class MuonRingFitter(Component):
         default_value=list(FIT_METHOD_BY_NAME.keys())[0],
     ).tag(config=True)
 
-    def __call__(self, x, y, img, mask):
+    def __call__(self, fov_lon, fov_lat, img, mask):
         """allows any fit to be called in form of
-            MuonRingFitter(fit_method = "name of the fit")
+        MuonRingFitter(fit_method = "name of the fit")
         """
         fit_function = FIT_METHOD_BY_NAME[self.fit_method]
-        radius, center_x, center_y = fit_function(x, y, img, mask)
+        radius, center_fov_lon, center_fov_lat = fit_function(
+            fov_lon, fov_lat, img, mask
+        )
 
         return MuonRingContainer(
-            center_x=center_x,
-            center_y=center_y,
+            center_fov_lon=center_fov_lon,
+            center_fov_lat=center_fov_lat,
             radius=radius,
-            center_phi=np.arctan2(center_y, center_x),
-            center_distance=np.sqrt(center_x ** 2 + center_y ** 2),
+            center_phi=np.arctan2(center_fov_lat, center_fov_lon),
+            center_distance=np.sqrt(center_fov_lon**2 + center_fov_lat**2),
         )
diff --git a/ctapipe/io/datawriter.py b/ctapipe/io/datawriter.py
index 8ddf8cae..d03d7c4b 100644
--- a/ctapipe/io/datawriter.py
+++ b/ctapipe/io/datawriter.py
@@ -194,6 +194,10 @@ class DataWriter(Component):
         help="Store DL2 stereo shower parameters if available", default_value=False
     ).tag(config=True)
 
+    write_muon_parameters = Bool(
+        help="Store muon parameters if available", default_value=False
+    ).tag(config=True)
+
     compression_level = Int(
         help="compression level, 0=None, 9=maximum", default_value=5, min=0, max=9
     ).tag(config=True)
@@ -326,6 +330,9 @@ class DataWriter(Component):
             self._write_dl2_telescope_events(self._writer, event)
             self._write_dl2_stereo_event(self._writer, event)
 
+        if self.write_muon_parameters:
+            self._write_muon_telescope_events(self._writer, event)
+
     def finish(self):
         """called after all events are done"""
         self.log.info("Finishing output")
@@ -398,6 +405,7 @@ class DataWriter(Component):
             self.write_images,
             self.write_showers,
             self.write_waveforms,
+            self.write_muon_parameters,
         ]
         if not any(writable_things):
             raise ToolConfigurationError(
@@ -703,6 +711,18 @@ class DataWriter(Component):
                         ],
                     )
 
+    def _write_muon_telescope_events(
+        self, writer: TableWriter, event: ArrayEventContainer
+    ):
+
+        for tel_id, muon in event.muon.tel.items():
+            table_name = self.table_name(tel_id)
+            tel_index = _get_tel_index(event, tel_id)
+            writer.write(
+                f"dl1/event/telescope/muon/{table_name}",
+                [tel_index, muon.ring, muon.parameters, muon.efficiency],
+            )
+
     def _write_dl2_telescope_events(
         self, writer: TableWriter, event: ArrayEventContainer
     ):
diff --git a/ctapipe/io/hdf5eventsource.py b/ctapipe/io/hdf5eventsource.py
index c44bdaf8..9c82488a 100644
--- a/ctapipe/io/hdf5eventsource.py
+++ b/ctapipe/io/hdf5eventsource.py
@@ -23,6 +23,10 @@ from ..containers import (
     IntensityStatisticsContainer,
     LeakageContainer,
     MorphologyContainer,
+    MuonEfficiencyContainer,
+    MuonParametersContainer,
+    MuonRingContainer,
+    MuonTelescopeContainer,
     ObservationBlockContainer,
     ParticleClassificationContainer,
     PeakTimeStatisticsContainer,
@@ -248,10 +252,10 @@ class HDF5EventSource(EventSource):
                 return False
 
             # we can now read both R1 and DL1
+            has_muons = "/dl1/event/telescope/muon" in f.root
             datalevels = set(metadata["CTA PRODUCT DATA LEVELS"].split(","))
-            if not datalevels.intersection(
-                ("R1", "DL1_IMAGES", "DL1_PARAMETERS", "DL2")
-            ):
+            datalevels = datalevels & {"R1", "DL1_IMAGES", "DL1_PARAMETERS", "DL2"}
+            if not datalevels and not has_muons:
                 return False
 
         return True
@@ -273,6 +277,13 @@ class HDF5EventSource(EventSource):
                 return True
         return False
 
+    @lazyproperty
+    def has_muon_parameters(self):
+        """
+        True for files that contain muon parameters
+        """
+        return "/dl1/event/telescope/muon" in self.file_.root
+
     @property
     def subarray(self):
         return self._subarray
@@ -461,6 +472,19 @@ class HDF5EventSource(EventSource):
                     for table in self.file_.root.dl1.event.telescope.parameters
                 }
 
+        if self.has_muon_parameters:
+            muon_readers = {
+                table.name: self.reader.read(
+                    f"/dl1/event/telescope/muon/{table.name}",
+                    containers=[
+                        MuonRingContainer,
+                        MuonParametersContainer,
+                        MuonEfficiencyContainer,
+                    ],
+                )
+                for table in self.file_.root.dl1.event.telescope.muon
+            }
+
         dl2_readers = {}
         if DL2_SUBARRAY_GROUP in self.file_.root:
             dl2_group = self.file_.root[DL2_SUBARRAY_GROUP]
@@ -595,6 +619,7 @@ class HDF5EventSource(EventSource):
 
             for tel_id in data.trigger.tel.keys():
                 key = f"tel_{tel_id:03d}"
+
                 if self.allowed_tels and tel_id not in self.allowed_tels:
                     continue
 
@@ -625,11 +650,12 @@ class HDF5EventSource(EventSource):
                                 "reconstructed image table."
                             )
                             continue
+
                         simulated_image_row = next(simulated_image_iterators[key])
                         simulated.true_image = simulated_image_row["true_image"]
 
                 if DataLevel.DL1_PARAMETERS in self.datalevels:
-                    if f"tel_{tel_id:03d}" not in param_readers:
+                    if key not in param_readers:
                         logger.debug(
                             f"Triggered telescope {tel_id} is missing "
                             "from the parameters table."
@@ -649,16 +675,14 @@ class HDF5EventSource(EventSource):
                         peak_time_statistics=params[6],
                     )
                     if self.has_simulated_dl1:
-                        if f"tel_{tel_id:03d}" not in simulated_param_readers:
+                        if key not in simulated_param_readers:
                             logger.debug(
                                 f"Triggered telescope {tel_id} is missing "
                                 "from the simulated parameters table, but was "
                                 "present at the reconstructed parameters table."
                             )
                             continue
-                        simulated_params = next(
-                            simulated_param_readers[f"tel_{tel_id:03d}"]
-                        )
+                        simulated_params = next(simulated_param_readers[key])
                         simulated.true_parameters = ImageParametersContainer(
                             hillas=simulated_params[0],
                             leakage=simulated_params[1],
@@ -667,6 +691,14 @@ class HDF5EventSource(EventSource):
                             intensity_statistics=simulated_params[4],
                         )
 
+                if self.has_muon_parameters:
+                    ring, parameters, efficiency = next(muon_readers[key])
+                    data.muon.tel[tel_id] = MuonTelescopeContainer(
+                        ring=ring,
+                        parameters=parameters,
+                        efficiency=efficiency,
+                    )
+
                 for kind, algorithms in dl2_tel_readers.items():
                     c = getattr(data.dl2.tel[tel_id], kind)
                     for algorithm, readers in algorithms.items():
@@ -724,11 +756,13 @@ class HDF5EventSource(EventSource):
         # Same comments as to _fill_array_pointing apply
         pointing_group = self.file_.root.dl1.monitoring.telescope.pointing
         for tel_id in data.trigger.tel.keys():
+            key = f"tel_{tel_id:03d}"
+
             if self.allowed_tels and tel_id not in self.allowed_tels:
                 continue
 
-            tel_pointing_table = pointing_group[f"tel_{tel_id:03d}"]
-            closest_time_index = tel_pointing_finder[f"tel_{tel_id:03d}"].closest(
+            tel_pointing_table = pointing_group[key]
+            closest_time_index = tel_pointing_finder[key].closest(
                 data.trigger.tel[tel_id].time.mjd
             )
 
diff --git a/ctapipe/io/tableloader.py b/ctapipe/io/tableloader.py
index f769d6be..6e4b34f1 100644
--- a/ctapipe/io/tableloader.py
+++ b/ctapipe/io/tableloader.py
@@ -22,6 +22,7 @@ __all__ = ["TableLoader"]
 
 PARAMETERS_GROUP = "/dl1/event/telescope/parameters"
 IMAGES_GROUP = "/dl1/event/telescope/images"
+MUON_GROUP = "/dl1/event/telescope/muon"
 TRIGGER_TABLE = "/dl1/event/subarray/trigger"
 SHOWER_TABLE = "/simulation/event/subarray/shower"
 TRUE_IMAGES_GROUP = "/simulation/event/telescope/images"
@@ -177,6 +178,9 @@ class TableLoader(Component):
     load_dl1_parameters = traits.Bool(
         True, help="load reconstructed image parameters"
     ).tag(config=True)
+    load_dl1_muons = traits.Bool(False, help="load muon ring parameters").tag(
+        config=True
+    )
 
     load_dl2 = traits.Bool(False, help="load available dl2 stereo parameters").tag(
         config=True
@@ -435,6 +439,12 @@ class TableLoader(Component):
             )
             table = _merge_telescope_tables(table, parameters)
 
+        if self.load_dl1_muons:
+            muon_parameters = self._read_telescope_table(
+                MUON_GROUP, tel_id, start=tel_start, stop=tel_stop
+            )
+            table = _merge_telescope_tables(table, muon_parameters)
+
         if self.load_dl1_images:
             images = self._read_telescope_table(
                 IMAGES_GROUP, tel_id, start=tel_start, stop=tel_stop
diff --git a/ctapipe/resources/base_config.yaml b/ctapipe/resources/base_config.yaml
index 08fed8ab..539bfb0f 100644
--- a/ctapipe/resources/base_config.yaml
+++ b/ctapipe/resources/base_config.yaml
@@ -26,6 +26,7 @@ DataWriter:
   write_showers: false # store DL2 stereo geometry
   write_raw_waveforms: false # write R0 waveforms
   write_waveforms: false # write R1 waveforms
+  write_muon_parameters: false # analyze and write muon parameters
 
 # The CameraCalibrator takes data from R1 or DL0 to DL1a level, applying finer
 # calibration and turning waveforms into images. It is run only if DL1a images
@@ -111,3 +112,26 @@ ShowerProcessor:
         - [Positive width, "parameters.hillas.width.value > 0"]
         - [enough pixels, "parameters.morphology.n_pixels > 3"]
         - [not clipped, "parameters.leakage.intensity_width_2 < 0.8"]
+
+# The MuonProcessor analyses DL1 images for muons and writes muon parameters.
+# It is only run if DataWriter.write_muon_parameters=True.
+MuonProcessor:
+  # Choose which events should be considered for muon analysis
+  ImageParameterQuery:
+    # quality critieria should be a list of pairs of [name, filter_expression]
+    # The filter expression should contain the variable `dl1_params` which are the
+    # image parameters.
+    quality_criteria:
+      - ["min_pixels", "dl1_params.morphology.n_pixels > 100"]
+      - ["min_intensity", "dl1_params.hillas.intensity > 500"]
+
+  # Choose which rings should be processed to the actual muon parameters
+  RingQuery:
+    # quality critieria should be a list of pairs of [name, filter_expression]
+    # The filter expression can either contain the variable `parameters` 
+    # (muon ring parameters like containment), `ring` (muon ring position and radius) 
+    # or `mask` (pixel mask masking the rings pixels)
+    quality_criteria:
+      - ["radius_not_nan", "np.isfinite(ring.radius.value)"]
+      - ["min_pixels", "np.count_nonzero(mask) > 50"]
+      - ["ring_containment", "parameters.containment > 0.5"]
diff --git a/ctapipe/tools/merge.py b/ctapipe/tools/merge.py
index 7fd6d1fd..7ebeaed5 100644
--- a/ctapipe/tools/merge.py
+++ b/ctapipe/tools/merge.py
@@ -15,7 +15,7 @@ from ctapipe.utils.arrays import recarray_drop_columns
 from ..core import Provenance, Tool, traits
 from ..core.traits import Bool, CInt, Set, Unicode, flag
 from ..instrument import SubarrayDescription
-from ..io import HDF5EventSource, HDF5TableWriter, get_hdf5_datalevels
+from ..io import HDF5EventSource, get_hdf5_datalevels
 from ..io import metadata as meta
 
 PROV = Provenance()
@@ -37,6 +37,7 @@ optional_nodes = {
     "/simulation/event/telescope/impact",
     "/dl1/event/telescope/parameters",
     "/dl1/event/telescope/images",
+    "/dl1/event/telescope/muon",
     "/dl2/event/telescope/geometry",
     "/dl2/event/telescope/impact",
     "/dl2/event/telescope/energy",
@@ -61,6 +62,7 @@ simulation_nodes = {
 nodes_with_tels = {
     "/dl1/monitoring/telescope/pointing",
     "/dl1/event/telescope/parameters",
+    "/dl1/event/telescope/muon",
     "/dl1/event/telescope/images",
     "/simulation/event/telescope/parameters",
     "/simulation/event/telescope/images",
@@ -72,6 +74,7 @@ image_nodes = {
 parameter_nodes = {
     "/simulation/event/telescope/parameters",
     "/dl1/event/telescope/parameters",
+    "/dl1/event/telescope/muon",
 }
 
 SIMULATED_IMAGE_GROUP = "/simulation/event/telescope/images"
@@ -262,7 +265,12 @@ class MergeTool(Tool):
                     self.usable_nodes.remove(node)
 
         # create output file with subarray from first file
-        self.first_subarray = SubarrayDescription.from_hdf(self.input_files[0])
+        self.first_subarray = SubarrayDescription.from_hdf(
+            self.input_files[0],
+            # focal length choice has no effect here, so use EQUIVALENT
+            # to support merging files that do not have EFFECTIVE
+            focal_length_choice="EQUIVALENT",
+        )
         if self.allowed_tels:
             self.first_subarray = self.first_subarray.select_subarray(
                 tel_ids=self.allowed_tels
@@ -286,7 +294,12 @@ class MergeTool(Tool):
             )
             return True
 
-        current_subarray = SubarrayDescription.from_hdf(file_path)
+        current_subarray = SubarrayDescription.from_hdf(
+            file_path,
+            # focal length choice has no effect here, so use EQUIVALENT
+            # to support merging files that do not have EFFECTIVE
+            focal_length_choice="EQUIVALENT",
+        )
         if self.allowed_tels:
             current_subarray = current_subarray.select_subarray(
                 tel_ids=self.allowed_tels
@@ -469,7 +482,6 @@ class MergeTool(Tool):
 
     def finish(self):
         datalevels = [d.name for d in get_hdf5_datalevels(self.output_file)]
-        self.output_file.close()
 
         activity = PROV.current_activity.provenance
         process_type_ = "Observation"
@@ -500,11 +512,8 @@ class MergeTool(Tool):
         )
 
         headers = reference.to_dict()
-
-        with HDF5TableWriter(
-            self.output_path, parent=self, mode="a", add_prefix=True
-        ) as writer:
-            meta.write_to_hdf5(headers, writer.h5file)
+        meta.write_to_hdf5(headers, self.output_file)
+        self.output_file.close()
 
 
 def main():
diff --git a/ctapipe/tools/muon_reconstruction.py b/ctapipe/tools/muon_reconstruction.py
deleted file mode 100644
index f9666226..00000000
--- a/ctapipe/tools/muon_reconstruction.py
+++ /dev/null
@@ -1,280 +0,0 @@
-import numpy as np
-from astropy.coordinates import SkyCoord
-from tqdm.auto import tqdm
-
-from ctapipe.image.muon import (
-    MuonIntensityFitter,
-    MuonRingFitter,
-    intensity_ratio_inside_ring,
-    mean_squared_error,
-    ring_completeness,
-    ring_containment,
-)
-
-from ..calib import CameraCalibrator
-from ..containers import MuonParametersContainer, TelEventIndexContainer
-from ..coordinates import CameraFrame, TelescopeFrame
-from ..core import Provenance, Tool, traits
-from ..image.cleaning import TailcutsImageCleaner
-from ..instrument import CameraGeometry
-from ..io import EventSource, HDF5TableWriter
-
-
-class MuonAnalysis(Tool):
-    """
-    Detect and extract muon ring parameters, and write the muon ring and
-    intensity parameters to an output table.
-
-    The resulting output can be read e.g. using for example
-    `pandas.read_hdf(filename, 'dl1/event/telescope/parameters/muon')`
-    """
-
-    name = "ctapipe-reconstruct-muons"
-    description = traits.Unicode(__doc__)
-
-    output = traits.Path(
-        directory_ok=False, allow_none=False, help="HDF5 output file name"
-    ).tag(config=True)
-
-    completeness_threshold = traits.FloatTelescopeParameter(
-        default_value=30.0, help="Threshold for calculating the ``ring_completeness``"
-    ).tag(config=True)
-
-    ratio_width = traits.FloatTelescopeParameter(
-        default_value=1.5,
-        help=(
-            "Ring width for intensity ratio"
-            " computation as multiple of pixel diameter"
-        ),
-    ).tag(config=True)
-
-    min_pixels = traits.IntTelescopeParameter(
-        help=(
-            "Minimum number of pixels after cleaning and ring finding"
-            "required to process an event"
-        ),
-        default_value=100,
-    ).tag(config=True)
-
-    pedestal = traits.FloatTelescopeParameter(
-        help="Pedestal noise rms", default_value=1.1
-    ).tag(config=True)
-
-    classes = [
-        CameraCalibrator,
-        TailcutsImageCleaner,
-        EventSource,
-        MuonRingFitter,
-        MuonIntensityFitter,
-    ]
-
-    aliases = {
-        ("i", "input"): "EventSource.input_url",
-        ("o", "output"): "MuonAnalysis.output",
-        ("m", "max-events"): "EventSource.max_events",
-        ("t", "allowed-tels"): "EventSource.allowed_tels",
-    }
-
-    def setup(self):
-        self.check_output(self.output)
-        self.source = EventSource(parent=self)
-        subarray = self.source.subarray
-
-        self.calib = CameraCalibrator(subarray=subarray, parent=self)
-        self.ring_fitter = MuonRingFitter(parent=self)
-        self.intensity_fitter = MuonIntensityFitter(subarray=subarray, parent=self)
-        self.cleaning = TailcutsImageCleaner(parent=self, subarray=subarray)
-        self.writer = HDF5TableWriter(
-            self.output, "", add_prefix=True, parent=self, mode="w"
-        )
-        self.pixels_in_tel_frame = {}
-        self.field_of_view = {}
-        self.pixel_widths = {}
-
-        for p in ["min_pixels", "pedestal", "ratio_width", "completeness_threshold"]:
-            getattr(self, p).attach_subarray(self.source.subarray)
-
-    def start(self):
-        for event in tqdm(self.source, desc="Processing events: "):
-            self.process_array_event(event)
-
-    def process_array_event(self, event):
-        self.calib(event)
-
-        for tel_id, dl1 in event.dl1.tel.items():
-            self.process_telescope_event(event.index, tel_id, dl1)
-
-        self.writer.write(
-            "sim/event/subarray/shower", [event.index, event.simulation.shower]
-        )
-
-    def process_telescope_event(self, event_index, tel_id, dl1):
-        event_id = event_index.event_id
-
-        if self.source.subarray.tel[tel_id].optics.n_mirrors != 1:
-            self.log.warn(
-                f"Skipping non-single mirror telescope {tel_id}"
-                " set --allowed_tels to get rid of this warning"
-            )
-            return
-
-        self.log.debug(f"Processing event {event_id}, telescope {tel_id}")
-        image = dl1.image
-        if dl1.image_mask is None:
-            dl1.image_mask = self.cleaning(tel_id, image)
-
-        if np.count_nonzero(dl1.image_mask) <= self.min_pixels.tel[tel_id]:
-            self.log.debug(
-                f"Skipping event {event_id}-{tel_id}:"
-                f" has less then {self.min_pixels.tel[tel_id]} pixels after cleaning"
-            )
-            return
-
-        x, y = self.get_pixel_coords(tel_id)
-
-        # iterative ring fit.
-        # First use cleaning pixels, then only pixels close to the ring
-        # three iterations seems to be enough for most rings
-        mask = dl1.image_mask
-        for i in range(3):
-            ring = self.ring_fitter(x, y, image, mask)
-            dist = np.sqrt((x - ring.center_x) ** 2 + (y - ring.center_y) ** 2)
-            mask = np.abs(dist - ring.radius) / ring.radius < 0.4
-
-        if np.count_nonzero(mask) <= self.min_pixels.tel[tel_id]:
-            self.log.debug(
-                f"Skipping event {event_id}-{tel_id}:"
-                f" Less then {self.min_pixels.tel[tel_id]} pixels on ring"
-            )
-            return
-
-        if np.isnan(
-            [ring.radius.value, ring.center_x.value, ring.center_y.value]
-        ).any():
-            self.log.debug(
-                f"Skipping event {event_id}-{tel_id}: Ring fit did not succeed"
-            )
-            return
-
-        parameters = self.calculate_muon_parameters(tel_id, image, dl1.image_mask, ring)
-        # intensity_fitter does not support a mask yet, set ignored pixels to 0
-        image[~mask] = 0
-
-        result = self.intensity_fitter(
-            tel_id,
-            ring.center_x,
-            ring.center_y,
-            ring.radius,
-            image,
-            pedestal=self.pedestal.tel[tel_id],
-        )
-
-        self.log.info(
-            f"Muon fit: r={ring.radius:.2f}"
-            f", width={result.width:.4f}"
-            f", efficiency={result.optical_efficiency:.2%}"
-        )
-
-        tel_event_index = TelEventIndexContainer(**event_index, tel_id=tel_id)
-
-        self.writer.write(
-            "dl1/event/telescope/parameters/muons",
-            [tel_event_index, ring, parameters, result],
-        )
-
-    def calculate_muon_parameters(self, tel_id, image, clean_mask, ring):
-        fov_radius = self.get_fov(tel_id)
-        x, y = self.get_pixel_coords(tel_id)
-
-        # add ring containment, not filled in fit
-        containment = ring_containment(
-            ring.radius, ring.center_x, ring.center_y, fov_radius
-        )
-
-        completeness = ring_completeness(
-            x,
-            y,
-            image,
-            ring.radius,
-            ring.center_x,
-            ring.center_y,
-            threshold=self.completeness_threshold.tel[tel_id],
-        )
-
-        pixel_width = self.get_pixel_width(tel_id)
-        intensity_ratio = intensity_ratio_inside_ring(
-            x[clean_mask],
-            y[clean_mask],
-            image[clean_mask],
-            ring.radius,
-            ring.center_x,
-            ring.center_y,
-            width=self.ratio_width.tel[tel_id] * pixel_width,
-        )
-
-        mse = mean_squared_error(
-            x[clean_mask],
-            y[clean_mask],
-            image[clean_mask],
-            ring.radius,
-            ring.center_x,
-            ring.center_y,
-        )
-
-        return MuonParametersContainer(
-            containment=containment,
-            completeness=completeness,
-            intensity_ratio=intensity_ratio,
-            mean_squared_error=mse,
-        )
-
-    def get_fov(self, tel_id):
-        """Guesstimate fov radius for telescope with id `tel_id`"""
-        # memoize fov calculation
-        if tel_id not in self.field_of_view:
-            cam = self.source.subarray.tel[tel_id].camera.geometry
-            border = cam.get_border_pixel_mask()
-
-            x, y = self.get_pixel_coords(tel_id)
-            self.field_of_view[tel_id] = np.sqrt(x[border] ** 2 + y[border] ** 2).mean()
-
-        return self.field_of_view[tel_id]
-
-    def get_pixel_width(self, tel_id):
-        """Guesstimate fov radius for telescope with id `tel_id`"""
-        # memoize fov calculation
-        if tel_id not in self.pixel_widths:
-            x, y = self.get_pixel_coords(tel_id)
-            self.pixel_widths[tel_id] = CameraGeometry.guess_pixel_width(x, y)
-
-        return self.pixel_widths[tel_id]
-
-    def get_pixel_coords(self, tel_id):
-        """Get pixel coords in telescope frame for telescope with id `tel_id`"""
-        # memoize transformation
-        if tel_id not in self.pixels_in_tel_frame:
-            telescope = self.source.subarray.tel[tel_id]
-            cam = telescope.camera.geometry
-            camera_frame = CameraFrame(
-                focal_length=telescope.optics.equivalent_focal_length,
-                rotation=cam.cam_rotation,
-            )
-            cam_coords = SkyCoord(x=cam.pix_x, y=cam.pix_y, frame=camera_frame)
-            tel_coord = cam_coords.transform_to(TelescopeFrame())
-            self.pixels_in_tel_frame[tel_id] = tel_coord
-
-        coords = self.pixels_in_tel_frame[tel_id]
-        return coords.fov_lon, coords.fov_lat
-
-    def finish(self):
-        Provenance().add_output_file(self.output, role="muon_efficiency_parameters")
-        self.writer.close()
-
-
-def main():
-    tool = MuonAnalysis()
-    tool.run()
-
-
-if __name__ == "__main__":
-    main()
diff --git a/ctapipe/tools/process.py b/ctapipe/tools/process.py
index da86c012..0a82f7eb 100644
--- a/ctapipe/tools/process.py
+++ b/ctapipe/tools/process.py
@@ -11,6 +11,7 @@ from ..core import QualityQuery, Tool
 from ..core.traits import Bool, classes_with_traits, flag
 from ..image import ImageCleaner, ImageModifier, ImageProcessor
 from ..image.extractor import ImageExtractor
+from ..image.muon import MuonProcessor
 from ..instrument import SoftwareTrigger
 from ..io import (
     DataLevel,
@@ -36,8 +37,14 @@ __all__ = ["ProcessorTool"]
 
 class ProcessorTool(Tool):
     """
-    Process data from lower-data levels up to DL1, including both image
-    extraction and optinally image parameterization
+    Process data from lower-data levels up to DL1 and DL2, including image
+    extraction and optionally image parameterization as well as muon analysis
+    and shower reconstruction.
+
+    Note that the muon analysis and shower reconstruction both depend on
+    parametrized images and therefore compute image parameters even if
+    DataWriter.write_parameters=False in case these are not already present
+    in the input file.
     """
 
     name = "ctapipe-process"
@@ -126,6 +133,12 @@ class ProcessorTool(Tool):
             "DataWriter.write_index_tables",
             "generate PyTables index tables for the parameter and image datasets",
         ),
+        **flag(
+            "write-muon-parameters",
+            "DataWriter.write_muon_parameters",
+            "store DL1/Event/Telescope muon parameters in output",
+            "don't store DL1/Event/Telescope muon parameters in output",
+        ),
         "camera-frame": (
             {"ImageProcessor": {"use_telescope_frame": False}},
             "Use camera frame for image parameters instead of telescope frame",
@@ -137,6 +150,7 @@ class ProcessorTool(Tool):
             CameraCalibrator,
             DataWriter,
             ImageProcessor,
+            MuonProcessor,
             ShowerProcessor,
             metadata.Instrument,
             metadata.Contact,
@@ -176,6 +190,7 @@ class ProcessorTool(Tool):
         self.write = self.enter_context(
             DataWriter(event_source=self.event_source, parent=self)
         )
+        self.process_muons = MuonProcessor(subarray=subarray, parent=self)
 
         # add ml reco classes if model paths were supplied via cli and not already configured
         reco_aliases = {
@@ -229,10 +244,15 @@ class ProcessorTool(Tool):
         if DataLevel.DL1_PARAMETERS in self.event_source.datalevels:
             return False
 
-        return self.write.write_parameters or self.should_compute_dl2
+        return (
+            self.write.write_parameters
+            or self.should_compute_dl2
+            or self.should_compute_muon_parameters
+        )
 
     @property
     def should_calibrate(self):
+        """returns true if data should be calibrated"""
         if self.force_recompute_dl1:
             return True
 
@@ -247,6 +267,14 @@ class ProcessorTool(Tool):
 
         return False
 
+    @property
+    def should_compute_muon_parameters(self):
+        """returns true if we should compute muon parameters info"""
+        if self.write.write_muon_parameters:
+            return True
+
+        return False
+
     def _write_processing_statistics(self):
         """write out the event selection stats, etc."""
         # NOTE: don't remove this, not part of DataWriter
@@ -279,6 +307,9 @@ class ProcessorTool(Tool):
         """
         self.log.info("(re)compute DL1: %s", self.should_compute_dl1)
         self.log.info("(re)compute DL2: %s", self.should_compute_dl2)
+        self.log.info(
+            "compute muon parameters: %s", self.should_compute_muon_parameters
+        )
         self.event_source.subarray.info(printer=self.log.info)
 
         for event in tqdm(
@@ -306,6 +337,9 @@ class ProcessorTool(Tool):
             if self.should_compute_dl1:
                 self.process_images(event)
 
+            if self.should_compute_muon_parameters:
+                self.process_muons(event)
+
             if self.should_compute_dl2:
                 self.process_shower(event)
 
diff --git a/docs/changes/2168.api.rst b/docs/changes/2168.api.rst
new file mode 100644
index 00000000..bb7c6057
--- /dev/null
+++ b/docs/changes/2168.api.rst
@@ -0,0 +1,21 @@
+Migrate muon analysis into the ``ctapipe-process`` tool:
+
+1. The former ``muon_reconstruction`` tool is dropped and all functionalities are transferred 
+   into the ``ctapipe-process`` tool.
+
+2. The ``process`` tool now has a ``write_muon_parameters`` flag which defaults to ``false``.
+   Muons are only analyzed and written if the flag is set. Analyzing muons requires DL1b image 
+   parameters so even if ``write_parameters`` is ``false`` and no parameters are present in the 
+   ``EventSource`` those are computed. 
+
+3. Two instances of ``QualityQuery``, ``MuonProcessor.ImageParameterQuery`` and ``MuonProcessor.RingQuery`` 
+   are added to the muon analysis to either preselect images according to image parameters and 
+   to select images according to the initial, geometrical ring fit for further processing. 
+   Deselected events or those where the muon analysis fails are being returned and written 
+   filled with ``NaN``s instead of being ignored.
+   Base configure options for the muon analysis were added to the ``base_config.yaml``.
+
+4. The ``DataWriter`` now writes the results of a muon analysis into ``/dl1/event/telescope/muon/tel_id``,
+   given ``write_moun_parameters`` is set to ``true``.
+
+5. Muon nodes were added to the ``HDF5EventSource``, the ``TableLoader`` and the ``ctapipe-merge`` tool.
diff --git a/docs/ctapipe_api/image/muon.rst b/docs/ctapipe_api/image/muon.rst
index 625576b6..90899a54 100644
--- a/docs/ctapipe_api/image/muon.rst
+++ b/docs/ctapipe_api/image/muon.rst
@@ -8,6 +8,8 @@ Algorithms for detecting and characterizing muon images
 
 .. automodapi:: ctapipe.image.muon
 
+.. automodapi:: ctapipe.image.muon.processor
+
 .. automodapi:: ctapipe.image.muon.intensity_fitter
 
 .. automodapi:: ctapipe.image.muon.ring_fitter
diff --git a/setup.cfg b/setup.cfg
index c73a2528..4c8698c4 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -89,7 +89,6 @@ console_scripts =
     ctapipe-info = ctapipe.tools.info:main
     ctapipe-dump-triggers = ctapipe.tools.dump_triggers:main
     ctapipe-dump-instrument=ctapipe.tools.dump_instrument:main
-    ctapipe-reconstruct-muons = ctapipe.tools.muon_reconstruction:main
     ctapipe-display-dl1 = ctapipe.tools.display_dl1:main
     ctapipe-process = ctapipe.tools.process:main
     ctapipe-merge = ctapipe.tools.merge:main


Test patch:
diff --git a/ctapipe/conftest.py b/ctapipe/conftest.py
index dec57bff..0efc616e 100644
--- a/ctapipe/conftest.py
+++ b/ctapipe/conftest.py
@@ -443,14 +443,41 @@ def dl1_muon_file(dl1_tmp_path):
             f"--input={infile}",
             f"--output={output}",
             "--write-images",
-            "--DataWriter.write_parameters=False",
-            "--DataWriter.Contact.name=αℓℓ the äüöß",
+            "--no-write-parameters",
             "--SimTelEventSource.focal_length_choice=EQUIVALENT",
         ]
         assert run_tool(ProcessorTool(), argv=argv, cwd=dl1_tmp_path) == 0
         return output
 
 
+@pytest.fixture(scope="session")
+def dl1_muon_output_file(dl1_tmp_path, dl1_muon_file):
+    """
+    DL1 file containing images, parameters and muon ring parameters
+    from a muon simulation set.
+    """
+    from ctapipe.tools.process import ProcessorTool
+
+    output = dl1_tmp_path / "muon_output.dl1.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"--input={dl1_muon_file}",
+            f"--output={output}",
+            "--no-write-images",
+            "--no-write-parameters",
+            "--write-muon-parameters",
+            "--HDF5EventSource.focal_length_choice=EQUIVALENT",
+            "--max-events=30",
+        ]
+        assert run_tool(ProcessorTool(), argv=argv, cwd=dl1_tmp_path) == 0
+        return output
+
+
 @pytest.fixture(scope="session")
 def dl1_proton_file(dl1_tmp_path, prod5_proton_simtel_path):
     """
diff --git a/ctapipe/image/muon/tests/test_processor.py b/ctapipe/image/muon/tests/test_processor.py
new file mode 100644
index 00000000..e2a7211e
--- /dev/null
+++ b/ctapipe/image/muon/tests/test_processor.py
@@ -0,0 +1,34 @@
+"""
+Tests for MuonProcessor functionality
+"""
+import numpy as np
+
+from ctapipe.image import ImageProcessor
+from ctapipe.image.muon import MuonProcessor
+from ctapipe.io import EventSource
+
+
+def test_processor(dl1_muon_file):
+    """Test that the MuonProcessor component analyses events or
+    at least provides defaults."""
+
+    with EventSource(dl1_muon_file, focal_length_choice="EQUIVALENT") as source:
+        image_processor = ImageProcessor(source.subarray)
+        muon_processor = MuonProcessor(source.subarray)
+
+        efficiencies = []
+        for event in source:
+            image_processor(event)
+            muon_processor(event)
+            for tel_id in event.dl1.tel:
+                efficiencies.append(
+                    event.muon.tel[tel_id].efficiency.optical_efficiency
+                )
+
+        assert len(efficiencies) > 0  # Assert there were events analyzed
+        assert np.any(
+            np.isfinite(efficiencies)
+        )  # Assert at least some were completely analyzed
+        assert np.all(
+            np.logical_or(np.isfinite(efficiencies), np.isnan(efficiencies))
+        )  # Assert all were at least provided defaults
diff --git a/ctapipe/image/muon/tests/test_ring_fitter.py b/ctapipe/image/muon/tests/test_ring_fitter.py
index f71af2ed..6886be6e 100644
--- a/ctapipe/image/muon/tests/test_ring_fitter.py
+++ b/ctapipe/image/muon/tests/test_ring_fitter.py
@@ -42,6 +42,6 @@ def test_MuonRingFitter(method, prod5_mst_flashcam):
     print(fit_result)
     print(center_xs, center_ys, radius)
 
-    assert u.isclose(fit_result.center_x, center_xs, 5e-2)
-    assert u.isclose(fit_result.center_y, center_ys, 5e-2)
+    assert u.isclose(fit_result.center_fov_lon, center_xs, 5e-2)
+    assert u.isclose(fit_result.center_fov_lat, center_ys, 5e-2)
     assert u.isclose(fit_result.radius, radius, 5e-2)
diff --git a/ctapipe/io/tests/test_hdf5eventsource.py b/ctapipe/io/tests/test_hdf5eventsource.py
index 029922ce..bae1eb13 100644
--- a/ctapipe/io/tests/test_hdf5eventsource.py
+++ b/ctapipe/io/tests/test_hdf5eventsource.py
@@ -13,11 +13,15 @@ def test_is_not_compatible(prod5_gamma_simtel_path):
     assert not HDF5EventSource.is_compatible(prod5_gamma_simtel_path)
 
 
-@pytest.mark.parametrize("compatible_file", ["dl1_file", "dl2_only_file"])
+@pytest.mark.parametrize(
+    "compatible_file", ["dl1_file", "dl2_only_file", "dl1_muon_output_file"]
+)
 def test_is_compatible(compatible_file, request):
+    old_files = {"dl1_muon_output_file"}
+    focal_length = "EQUIVALENT" if compatible_file in old_files else "EFFECTIVE"
     file = request.getfixturevalue(compatible_file)
     assert HDF5EventSource.is_compatible(file)
-    with EventSource(input_url=file) as source:
+    with EventSource(input_url=file, focal_length_choice=focal_length) as source:
         assert isinstance(source, HDF5EventSource)
 
 
diff --git a/ctapipe/io/tests/test_table_loader.py b/ctapipe/io/tests/test_table_loader.py
index 843267f0..095cc3ea 100644
--- a/ctapipe/io/tests/test_table_loader.py
+++ b/ctapipe/io/tests/test_table_loader.py
@@ -87,6 +87,25 @@ def test_telescope_events_for_tel_id(dl1_file):
     assert not table_loader.h5file.isopen
 
 
+def test_telescope_muon_events_for_tel_id(dl1_muon_output_file):
+    """Test loading muon data for a single telescope"""
+    from ctapipe.io.tableloader import TableLoader
+
+    with TableLoader(
+        dl1_muon_output_file,
+        load_dl1_muons=True,
+        load_dl1_parameters=False,
+        focal_length_choice="EQUIVALENT",
+    ) as table_loader:
+        table = table_loader.read_telescope_events([1])
+        assert "muonring_radius" in table.colnames
+        assert "muonparameters_containment" in table.colnames
+        assert "muonefficiency_optical_efficiency" in table.colnames
+        assert np.all(table["tel_id"] == 1)
+
+    assert not table_loader.h5file.isopen
+
+
 def test_load_instrument(dl1_file):
     """Test joining instrument data onto telescope events"""
     from ctapipe.io.tableloader import TableLoader
diff --git a/ctapipe/tools/tests/test_merge.py b/ctapipe/tools/tests/test_merge.py
index 4b7d57ca..b9aec27c 100644
--- a/ctapipe/tools/tests/test_merge.py
+++ b/ctapipe/tools/tests/test_merge.py
@@ -11,6 +11,7 @@ from astropy.utils.diff import report_diff_values
 from ctapipe.core import run_tool
 from ctapipe.io import TableLoader
 from ctapipe.io.astropy_helpers import read_table
+from ctapipe.io.tests.test_astropy_helpers import assert_table_equal
 from ctapipe.tools.process import ProcessorTool
 
 try:
@@ -27,7 +28,7 @@ def run_stage1(input_path, cwd, output_path=None):
             tempfile.NamedTemporaryFile(suffix=".dl1.h5", dir=cwd).name
         ).absolute()
 
-    ret = run_tool(
+    run_tool(
         ProcessorTool(),
         argv=[
             f"--config={config}",
@@ -39,20 +40,20 @@ def run_stage1(input_path, cwd, output_path=None):
             "--max-events=5",
         ],
         cwd=cwd,
+        raises=True,
     )
-    assert ret == 0, "Running stage1 failed"
 
 
 def test_simple(tmp_path, dl1_file, dl1_proton_file):
     from ctapipe.tools.merge import MergeTool
 
     output = tmp_path / "merged_simple.dl1.h5"
-    ret = run_tool(
+    run_tool(
         MergeTool(),
         argv=[str(dl1_file), str(dl1_proton_file), f"--output={output}", "--overwrite"],
         cwd=tmp_path,
+        raises=True,
     )
-    assert ret == 0
     run_stage1(output, cwd=tmp_path)
 
 
@@ -67,7 +68,7 @@ def test_pattern(tmp_path: Path, dl1_file, dl1_proton_file):
         shutil.copy(f, tmp_path)
 
     output = tmp_path / "merged_pattern.dl1.h5"
-    ret = run_tool(
+    run_tool(
         tool=MergeTool(),
         argv=[
             "-i",
@@ -78,8 +79,8 @@ def test_pattern(tmp_path: Path, dl1_file, dl1_proton_file):
             "--overwrite",
         ],
         cwd=tmp_path,
+        raises=True,
     )
-    assert ret == 0
     run_stage1(output, cwd=tmp_path)
 
 
@@ -88,7 +89,7 @@ def test_skip_images(tmp_path, dl1_file, dl1_proton_file):
 
     # create a second file so we can test the patterns
     output = tmp_path / "merged_no_images.dl1.h5"
-    ret = run_tool(
+    run_tool(
         MergeTool(),
         argv=[
             str(dl1_file),
@@ -98,8 +99,8 @@ def test_skip_images(tmp_path, dl1_file, dl1_proton_file):
             "--overwrite",
         ],
         cwd=tmp_path,
+        raises=True,
     )
-    assert ret == 0
 
     with tables.open_file(output, "r") as f:
         assert "images" not in f.root.dl1.event.telescope
@@ -124,8 +125,7 @@ def test_allowed_tels(tmp_path, dl1_file, dl1_proton_file):
     for tel_id in allowed_tels:
         argv.append(f"--allowed-tels={tel_id}")
 
-    ret = run_tool(MergeTool(), argv=argv, cwd=tmp_path)
-    assert ret == 0
+    run_tool(MergeTool(), argv=argv, cwd=tmp_path, raises=True)
 
     s = SubarrayDescription.from_hdf(output)
     assert s.tel.keys() == allowed_tels
@@ -143,15 +143,15 @@ def test_dl2(tmp_path, dl2_shower_geometry_file, dl2_proton_geometry_file):
     from ctapipe.tools.merge import MergeTool
 
     output = tmp_path / "merged.dl2.h5"
-    ret = run_tool(
+    run_tool(
         MergeTool(),
         argv=[
             f"--output={output}",
             str(dl2_shower_geometry_file),
             str(dl2_proton_geometry_file),
         ],
+        raises=True,
     )
-    assert ret == 0, f"Running merge for dl2 files failed with exit code {ret}"
 
     table1 = read_table(
         dl2_shower_geometry_file, "/dl2/event/subarray/geometry/HillasReconstructor"
@@ -190,3 +190,26 @@ def test_dl2(tmp_path, dl2_shower_geometry_file, dl2_proton_geometry_file):
     assert "true_impact_distance" in tel_events.colnames
     # regression test for #2051
     assert "HillasReconstructor_tel_impact_distance" in tel_events.colnames
+
+
+def test_muon(tmp_path, dl1_muon_output_file):
+    from ctapipe.tools.merge import MergeTool
+
+    output = tmp_path / "muon_merged.dl2.h5"
+    run_tool(
+        MergeTool(),
+        argv=[
+            f"--output={output}",
+            str(dl1_muon_output_file),
+            str(dl1_muon_output_file),
+        ],
+        raises=True,
+    )
+
+    table = read_table(output, "/dl1/event/telescope/muon/tel_001")
+    input_table = read_table(dl1_muon_output_file, "/dl1/event/telescope/muon/tel_001")
+
+    n_input = len(input_table)
+    assert len(table) == 2 * n_input
+    assert_table_equal(table[:n_input], input_table)
+    assert_table_equal(table[n_input:], input_table)
diff --git a/ctapipe/tools/tests/test_process.py b/ctapipe/tools/tests/test_process.py
index 6884b51e..e129fa56 100644
--- a/ctapipe/tools/tests/test_process.py
+++ b/ctapipe/tools/tests/test_process.py
@@ -6,6 +6,7 @@ Test ctapipe-process on a few different use cases
 
 from subprocess import CalledProcessError
 
+import astropy.units as u
 import numpy as np
 import pandas as pd
 import pytest
@@ -14,13 +15,14 @@ from numpy.testing import assert_allclose, assert_array_equal
 
 from ctapipe.core import run_tool
 from ctapipe.instrument.subarray import SubarrayDescription
-from ctapipe.io import TableLoader, read_table
+from ctapipe.io import EventSource, TableLoader, read_table
 from ctapipe.io.tests.test_event_source import DummyEventSource
 from ctapipe.tools.process import ProcessorTool
 from ctapipe.tools.quickstart import CONFIGS_TO_WRITE, QuickStartTool
 from ctapipe.utils import get_dataset_path, resource_file
 
 GAMMA_TEST_LARGE = get_dataset_path("gamma_test_large.simtel.gz")
+LST_MUONS = get_dataset_path("lst_muons.simtel.zst")
 
 
 @pytest.mark.parametrize(
@@ -72,20 +74,18 @@ def test_stage_1_dl1(tmp_path, dl1_image_file, dl1_parameters_file):
 
     # DL1A file as input
     dl1b_from_dl1a_file = tmp_path / "dl1b_fromdl1a.dl1.h5"
-    assert (
-        run_tool(
-            ProcessorTool(),
-            argv=[
-                f"--config={config}",
-                f"--input={dl1_image_file}",
-                f"--output={dl1b_from_dl1a_file}",
-                "--camera-frame",
-                "--write-parameters",
-                "--overwrite",
-            ],
-            cwd=tmp_path,
-        )
-        == 0
+    run_tool(
+        ProcessorTool(),
+        argv=[
+            f"--config={config}",
+            f"--input={dl1_image_file}",
+            f"--output={dl1b_from_dl1a_file}",
+            "--camera-frame",
+            "--write-parameters",
+            "--overwrite",
+        ],
+        cwd=tmp_path,
+        raises=True,
     )
 
     # check tables were written
@@ -177,18 +177,16 @@ def test_stage_2_from_simtel(tmp_path):
     config = resource_file("stage2_config.json")
     output = tmp_path / "test_stage2_from_simtel.DL2.h5"
 
-    assert (
-        run_tool(
-            ProcessorTool(),
-            argv=[
-                f"--config={config}",
-                "--input=dataset://gamma_prod5.simtel.zst",
-                f"--output={output}",
-                "--overwrite",
-            ],
-            cwd=tmp_path,
-        )
-        == 0
+    run_tool(
+        ProcessorTool(),
+        argv=[
+            f"--config={config}",
+            "--input=dataset://gamma_prod5.simtel.zst",
+            f"--output={output}",
+            "--overwrite",
+        ],
+        cwd=tmp_path,
+        raises=True,
     )
 
     # check tables were written
@@ -210,18 +208,16 @@ def test_stage_2_from_dl1_images(tmp_path, dl1_image_file):
     config = resource_file("stage2_config.json")
     output = tmp_path / "test_stage2_from_dl1image.DL2.h5"
 
-    assert (
-        run_tool(
-            ProcessorTool(),
-            argv=[
-                f"--config={config}",
-                f"--input={dl1_image_file}",
-                f"--output={output}",
-                "--overwrite",
-            ],
-            cwd=tmp_path,
-        )
-        == 0
+    run_tool(
+        ProcessorTool(),
+        argv=[
+            f"--config={config}",
+            f"--input={dl1_image_file}",
+            f"--output={output}",
+            "--overwrite",
+        ],
+        cwd=tmp_path,
+        raises=True,
     )
 
     # check tables were written
@@ -235,18 +231,16 @@ def test_stage_2_from_dl1_params(tmp_path, dl1_parameters_file):
     config = resource_file("stage2_config.json")
     output = tmp_path / "test_stage2_from_dl1param.DL2.h5"
 
-    assert (
-        run_tool(
-            ProcessorTool(),
-            argv=[
-                f"--config={config}",
-                f"--input={dl1_parameters_file}",
-                f"--output={output}",
-                "--overwrite",
-            ],
-            cwd=tmp_path,
-        )
-        == 0
+    run_tool(
+        ProcessorTool(),
+        argv=[
+            f"--config={config}",
+            f"--input={dl1_parameters_file}",
+            f"--output={output}",
+            "--overwrite",
+        ],
+        cwd=tmp_path,
+        raises=True,
     )
 
     # check tables were written
@@ -260,20 +254,18 @@ def test_ml_preprocessing_from_simtel(tmp_path):
     config = resource_file("ml_preprocessing_config.json")
     output = tmp_path / "test_ml_preprocessing.DL1DL2.h5"
 
-    assert (
-        run_tool(
-            ProcessorTool(),
-            argv=[
-                f"--config={config}",
-                f"--input={GAMMA_TEST_LARGE}",
-                f"--output={output}",
-                "--max-events=5",
-                "--overwrite",
-                "--SimTelEventSource.focal_length_choice=EQUIVALENT",
-            ],
-            cwd=tmp_path,
-        )
-        == 0
+    run_tool(
+        ProcessorTool(),
+        argv=[
+            f"--config={config}",
+            f"--input={GAMMA_TEST_LARGE}",
+            f"--output={output}",
+            "--max-events=5",
+            "--overwrite",
+            "--SimTelEventSource.focal_length_choice=EQUIVALENT",
+        ],
+        cwd=tmp_path,
+        raises=True,
     )
 
     # check tables were written
@@ -294,19 +286,17 @@ def test_image_modifications(tmp_path, dl1_image_file):
     noise_config = resource_file("image_modification_config.json")
 
     dl1_modified = tmp_path / "dl1_modified.dl1.h5"
-    assert (
-        run_tool(
-            ProcessorTool(),
-            argv=[
-                f"--config={noise_config}",
-                f"--input={dl1_image_file}",
-                f"--output={dl1_modified}",
-                "--write-parameters",
-                "--overwrite",
-            ],
-            cwd=tmp_path,
-        )
-        == 0
+    run_tool(
+        ProcessorTool(),
+        argv=[
+            f"--config={noise_config}",
+            f"--input={dl1_image_file}",
+            f"--output={dl1_modified}",
+            "--write-parameters",
+            "--overwrite",
+        ],
+        cwd=tmp_path,
+        raises=True,
     )
     modified_images = read_table(dl1_modified, "/dl1/event/telescope/images/tel_025")
     # Test that significantly more light is recorded (bias in dim pixels)
@@ -433,6 +423,52 @@ def test_read_from_simtel_and_dl1(prod5_proton_simtel_path, tmp_path):
     )
 
 
+def test_muon_reconstruction_simtel(tmp_path):
+    """ensure processor tool generates expected output when used to analyze muons"""
+    muon_simtel_output_file = tmp_path / "muon_reco_on_simtel.h5"
+    run_tool(
+        ProcessorTool(),
+        argv=[
+            f"--input={LST_MUONS}",
+            f"--output={muon_simtel_output_file}",
+            "--SimTelEventSource.focal_length_choice=EQUIVALENT",
+            "--overwrite",
+            "--write-muon-parameters",
+        ],
+        cwd=tmp_path,
+        raises=True,
+    )
+
+    table = read_table(muon_simtel_output_file, "/dl1/event/telescope/muon/tel_001")
+    assert len(table) > 20
+    assert np.count_nonzero(np.isfinite(table["muonring_radius"])) > 0
+    assert np.all(
+        np.logical_or(
+            np.isfinite(table["muonring_radius"]),
+            np.isnan(table["muonring_radius"]),
+        )
+    )
+
+    with EventSource(
+        muon_simtel_output_file, focal_length_choice="EQUIVALENT"
+    ) as source:
+        radius = table["muonring_radius"].quantity
+        efficiency = table["muonefficiency_optical_efficiency"]
+        completeness = table["muonparameters_completeness"]
+
+        for event in source:
+            muon = event.muon.tel[1]
+            assert u.isclose(muon.ring.radius, radius[event.count], equal_nan=True)
+            assert np.isclose(
+                muon.parameters.completeness, completeness[event.count], equal_nan=True
+            )
+            assert np.isclose(
+                muon.efficiency.optical_efficiency,
+                efficiency[event.count],
+                equal_nan=True,
+            )
+
+
 def test_plugin_help(capsys):
     ProcessorTool().print_help(classes=True)
     captured = capsys.readouterr()
diff --git a/ctapipe/tools/tests/test_tools.py b/ctapipe/tools/tests/test_tools.py
index 747e3d78..c31c7736 100644
--- a/ctapipe/tools/tests/test_tools.py
+++ b/ctapipe/tools/tests/test_tools.py
@@ -5,73 +5,20 @@ import subprocess
 import sys
 
 import matplotlib as mpl
-import numpy as np
 
 # pylint: disable=C0103,C0116,C0415
 import pytest
-import tables
 
 from ctapipe.core import run_tool
 from ctapipe.core.tool import ToolConfigurationError
 from ctapipe.utils import get_dataset_path
 
 GAMMA_TEST_LARGE = get_dataset_path("gamma_test_large.simtel.gz")
-LST_MUONS = get_dataset_path("lst_muons.simtel.zst")
 PROD5B_PATH = get_dataset_path(
     "gamma_20deg_0deg_run2___cta-prod5-paranal_desert-2147m-Paranal-dark_cone10-100evts.simtel.zst"
 )
 
 
-def test_muon_reconstruction_simtel(tmp_path):
-    from ctapipe.tools.muon_reconstruction import MuonAnalysis
-
-    muon_simtel_output_file = tmp_path / "muon_reco_on_simtel.h5"
-    assert (
-        run_tool(
-            MuonAnalysis(),
-            argv=[
-                f"--input={LST_MUONS}",
-                f"--output={muon_simtel_output_file}",
-                "--SimTelEventSource.focal_length_choice=EQUIVALENT",
-                "--overwrite",
-            ],
-            cwd=tmp_path,
-        )
-        == 0
-    )
-
-    with tables.open_file(muon_simtel_output_file) as t:
-        table = t.root.dl1.event.telescope.parameters.muons[:]
-        assert len(table) > 20
-        assert np.count_nonzero(np.isnan(table["muonring_radius"])) == 0
-
-
-def test_muon_reconstruction_dl1(tmp_path, dl1_muon_file):
-    from ctapipe.tools.muon_reconstruction import MuonAnalysis
-
-    muon_dl1_output_file = tmp_path / "muon_reco_on_dl1a.h5"
-    assert (
-        run_tool(
-            MuonAnalysis(),
-            argv=[
-                f"--input={dl1_muon_file}",
-                f"--output={muon_dl1_output_file}",
-                "--HDF5EventSource.focal_length_choice=EQUIVALENT",
-                "--overwrite",
-            ],
-            cwd=tmp_path,
-        )
-        == 0
-    )
-
-    with tables.open_file(muon_dl1_output_file) as t:
-        table = t.root.dl1.event.telescope.parameters.muons[:]
-        assert len(table) > 20
-        assert np.count_nonzero(np.isnan(table["muonring_radius"])) == 0
-
-    assert run_tool(MuonAnalysis(), ["--help-all"]) == 0
-
-
 def test_display_dl1(tmp_path, dl1_image_file, dl1_parameters_file):
     from ctapipe.tools.display_dl1 import DisplayDL1Calib