Repository: cta-observatory/ctapipe
Pull request #2736 (2025-04-09 15:00:10)

Muon fit improvement
The muon ring fitter needs to be improved.
The existing toy tests for muons are overly simplified.
Realistic simulations reveal significant discrepancies between the true ring and the fitted result (see the attached PNG). The fitted parameters are shown in green.

![Screenshot from 2025-04-09 16-55-44](https://github.com/user-attachments/assets/8076ce27-31c1-4ce4-ab29-847132f90681)
  
The second major component will be fitting the muon impact point while accounting for the muon’s inclination.

I stated with :
Enhance the toy model for muon ring simulation by incorporating user-defined asymmetry in both the x and y directions.
Update the muon ring fitter test to include this ring asymmetry.


Patch:
diff --git a/docs/api-reference/io/index.rst b/docs/api-reference/io/index.rst
index fc9b9218..9d4bd7f3 100644
--- a/docs/api-reference/io/index.rst
+++ b/docs/api-reference/io/index.rst
@@ -133,13 +133,11 @@ Writing Output Files
 ====================
 
 The `DataWriter` Component allows one to write a series of events (stored in
-`ctapipe.containers.ArrayEventContainer`) to an HDF5 format file
-following the data format described in :ref:`dataformat`.
-This includes all related datasets such as the instrument and simulation
-configuration information, simulated shower and image information, observed
-images and parameters and reconstruction information.
-
-It can be used in an event loop like:
+`ctapipe.containers.ArrayEventContainer`) to a standardized HDF5 format file
+following the data model (see :ref:`datamodels`). This includes all related datasets
+such as the instrument and simulation configuration information, simulated
+shower and image information, observed images and parameters and reconstruction
+information. It can be used in an event loop like:
 
 .. code-block:: python
 
diff --git a/docs/changes/2736.api.rst b/docs/changes/2736.api.rst
new file mode 100644
index 00000000..610f2205
--- /dev/null
+++ b/docs/changes/2736.api.rst
@@ -0,0 +1,12 @@
+The algorithms for muon ring fitting have been improved
+
+- Upgrade the taubin_circle_fit with weights to avoid bias toward rings with larger radii
+- Add the optional initial parameters : ring center and radius.
+- Add a combined ring fitting method (kundu_chaudhuri_taubin).
+- This new combined method is set to be the default method.
+- Updated the test with an additional fitting method, each method is tested separately.
+- Realistic muon ring simulation for all camera types implemented in the test with toy model.
+- Include the parameter fit errors in the tolerance metric.
+- Add error fields to the muon container.
+- The API for MuonFitter has been modified.
+- Minor documentation updates.
diff --git a/docs/user-guide/data_format/index.rst b/docs/user-guide/data_models/index.rst
similarity index 87%
rename from docs/user-guide/data_format/index.rst
rename to docs/user-guide/data_models/index.rst
index 8b82dc9d..408bc04e 100644
--- a/docs/user-guide/data_format/index.rst
+++ b/docs/user-guide/data_models/index.rst
@@ -1,7 +1,7 @@
-.. _dataformat:
+.. _datamodels:
 
 ***********
-Data Format
+Data Models
 ***********
 
 The ``ctapipe`` output files are using the HDF5 format files,
@@ -13,16 +13,8 @@ generated by ctapipe (e.g. the ``ctapipe-process`` tool which uses the
 
 The tables are written with `pytables <https://www.pytables.org>`_ (not ``h5py``).
 
-The in-memory data format is defined by the `~ctapipe.containers` structure.
-
-To read this data, there are two high-level APIs available:
-
-- Bulk-reading into astropy `~astropy.table.Table` objects using `~ctapipe.io.TableLoader`
-- Event-wise using the `~ctapipe.io.HDF5EventSource`
-
-
-DL1 Data Format
-===============
+DL1 Data Model
+==============
 
 This describes data that change per-event.
 The following datasets will be written to the group ``/dl1/event/`` in the  output file:
@@ -42,22 +34,22 @@ The following datasets will be written to the group ``/dl1/event/`` in the  outp
     * - ``telescope/``
       - Per-telescope Per-event information
       - (group)
-    * - ``telescope/images/tel_{TEL_ID:03d}``
-      - tables of telescope images (one per telescope)
-      - :py:class:`~ctapipe.containers.TelEventIndexContainer`, :py:class:`~ctapipe.containers.DL1CameraContainer`
     * - ``telescope/parameters/tel_{TEL_ID:03d}``
       - tables of image parameters (one per telescope)
       - :py:class:`~ctapipe.containers.TelEventIndexContainer`, :py:class:`~ctapipe.containers.ImageParametersContainer`
+    * - ``telescope/images/tel_{TEL_ID:03d}``
+      - tables of telescope images (one per telescope)
+      - :py:class:`~ctapipe.containers.TelEventIndexContainer`, :py:class:`~ctapipe.containers.DL1CameraContainer`
 
 
-DL2 Data Format
-===============
+DL2 Data Model
+==============
 
 This describes data that change per-event. The following datasets will be
-written to the group ``/dl2/event/subarray/<algorithm>/`` and or
-``/dl2/event/telescope/<algorithm>/``, one for each reconstruction algorithm in the
-output file, where ``<algorithm>`` is the identifier of the algorithm
-(e.g. "HillasReconstructor"):
+written to the group ``/dl2/event/stereo/<algorithm>/`` and or
+``/dl2/event/mono/<algorithm>/``, one for each reconstruction algorithm in the
+output file, where ``<algorithm>`` is the identifier of the algorithm (e.g.
+"Hillas"):
 
 .. list-table::
     :widths: 25 50 25
@@ -77,8 +69,8 @@ output file, where ``<algorithm>`` is the identifier of the algorithm
       - :py:class:`~ctapipe.containers.EventIndexContainer`, :py:class:`~ctapipe.containers.ParticleClassificationContainer`
 
 
-Simulation Data Format
-======================
+Simulation Data Model
+=====================
 
 .. list-table::
     :header-rows: 1
@@ -100,8 +92,8 @@ Simulation Data Format
       - :py:class:`~ctapipe.containers.SimulatedShowerDistribution`
 
 
-Configuration Data Format
-=========================
+Configuration Data Model
+========================
 
 The output file should also contain serializations of the instrument, observation (if
 applicable), simulation (if applicable) configuration information, written to the
diff --git a/docs/user-guide/index.rst b/docs/user-guide/index.rst
index dd649280..0474279f 100644
--- a/docs/user-guide/index.rst
+++ b/docs/user-guide/index.rst
@@ -9,6 +9,6 @@ User Guide
 
    getting-started
    tools/index
-   data_format/index
+   data_models/index
    examples_tutorials
    FAQ
diff --git a/src/ctapipe/containers.py b/src/ctapipe/containers.py
index 89b33ac0..8784a590 100644
--- a/src/ctapipe/containers.py
+++ b/src/ctapipe/containers.py
@@ -1164,6 +1164,13 @@ class MuonRingContainer(Container):
         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_fov_lon_err = Field(
+        nan * u.deg, "center (fov_lon) of the fitted muon ring error", unit=u.deg
+    )
+    center_fov_lat_err = Field(
+        nan * u.deg, "center (fov_lat) of the fitted muon ring error", unit=u.deg
+    )
+    radius_err = Field(nan * u.deg, "radius of the fitted muon ring error", unit=u.deg)
     center_phi = Field(
         nan * u.deg, "Angle of ring center within camera plane", unit=u.deg
     )
diff --git a/src/ctapipe/image/muon/__init__.py b/src/ctapipe/image/muon/__init__.py
index 4d352823..b1a2475d 100644
--- a/src/ctapipe/image/muon/__init__.py
+++ b/src/ctapipe/image/muon/__init__.py
@@ -4,13 +4,18 @@ from .features import (
     ring_completeness,
     ring_containment,
 )
-from .fitting import kundu_chaudhuri_circle_fit
+from .fitting import kundu_chaudhuri_circle_fit, taubin_circle_fit
 from .intensity_fitter import MuonIntensityFitter  # noqa: F401
 from .processor import MuonProcessor  # noqa: F401
-from .ring_fitter import MuonRingFitter  # noqa: F401
+from .ring_fitter import (
+    MuonRingFitter,  # noqa: F401
+    kundu_chaudhuri_taubin,
+)
 
 __all__ = [
     "kundu_chaudhuri_circle_fit",
+    "taubin_circle_fit",
+    "kundu_chaudhuri_taubin",
     "mean_squared_error",
     "intensity_ratio_inside_ring",
     "ring_completeness",
diff --git a/src/ctapipe/image/muon/fitting.py b/src/ctapipe/image/muon/fitting.py
index 64d8e4f8..0e2dee7e 100644
--- a/src/ctapipe/image/muon/fitting.py
+++ b/src/ctapipe/image/muon/fitting.py
@@ -1,3 +1,4 @@
+import astropy.units as u
 import numpy as np
 from astropy.units import Quantity
 
@@ -12,11 +13,11 @@ except ModuleNotFoundError:
     Minuit = None
 
 
-def kundu_chaudhuri_circle_fit(x, y, weights):
+def kundu_chaudhuri_circle_fit(x, y, weights, nan_errors_flag=False):
     """
-    Fast and reliable analytical circle fitting method previously used
-    in the H.E.S.S. experiment for muon identification
+    Fast and reliable analytical circle fitting method.
 
+    Previously used in the H.E.S.S. experiment for muon identification.
     Implementation based on :cite:p:`chaudhuri93`
 
     Parameters
@@ -27,7 +28,23 @@ def kundu_chaudhuri_circle_fit(x, y, weights):
         y coordinates of the points
     weights: array-like
         weights of the points
-
+    nan_errors_flag: bool
+        The flag defines whether errors are set to NaN.
+
+    Returns
+    -------
+    radius : astropy.units.Quantity
+        Fitted radius of the circle.
+    center_x : astropy.units.Quantity
+        Fitted x-coordinate of the circle center.
+    center_y : astropy.units.Quantity
+        Fitted y-coordinate of the circle center.
+    radius_err : astropy.units.Quantity
+        Fitted radius of the circle error.
+    center_x_err : astropy.units.Quantity
+        Fitted x-coordinate of the circle center error.
+    center_y_err : astropy.units.Quantity
+        Fitted y-coordinate of the circle center error.
     """
 
     weights_sum = np.sum(weights)
@@ -50,22 +67,79 @@ def kundu_chaudhuri_circle_fit(x, y, weights):
         np.sum(weights * ((center_x - x) ** 2 + (center_y - y) ** 2)) / weights_sum
     )
 
-    return radius, center_x, center_y
+    if nan_errors_flag:
+        radius_err = np.nan
+        center_x_err = np.nan
+        center_y_err = np.nan
+    else:
+        radius_err, center_x_err, center_y_err = naive_circle_fit_error_calculator(
+            x, y, weights, radius, center_x, center_y
+        )
+
+    return radius, center_x, center_y, radius_err, center_x_err, center_y_err
 
 
-def taubin_circle_fit(x, y, mask):
+def taubin_circle_fit(
+    x, y, mask, weights=None, r_initial=None, xc_initial=None, yc_initial=None
+):
     """
-    reference : Barcelona_Muons_TPA_final.pdf (slide 6)
+    Perform a Taubin circle fit with weights (optional).
+
+    The minimized loss function in this method tends to
+    maximize the radius of the ring, whereas using a simple
+    ring equation systematically results in a smaller radius.
+    Adding weights mitigates both effects and yields a more accurate fit.
 
     Parameters
     ----------
-    x: array-like or astropy quantity
-        x coordinates of the points
-    y: array-like or astropy quantity
-        y coordinates of the points
-    mask: array-like boolean
-        true for pixels surviving the cleaning
+    x : array-like or astropy.units.Quantity
+        x-coordinates of the points.
+    y : array-like or astropy.units.Quantity
+        y-coordinates of the points.
+    mask : array-like of bool
+        Boolean mask indicating which pixels survive the cleaning process.
+    weights : array-like of float, optional
+        Weights for the points. If not provided, all points are assigned equal weights (1).
+    r_initial : astropy.units.Quantity, optional
+        Initial guess for the radius of the circle. If not provided, it defaults to 1.1 deg.
+        1.1 deg. is the approximate Cherenkov photon angle produced by muons in the atmosphere
+        at La Palma altitude (2426 m a.s.l.), with momentum greater than 15 GeV.
+    xc_initial : astropy.units.Quantity, optional
+        Initial guess for the x-coordinate of the circle center. Defaults to 0.
+    yc_initial : astropy.units.Quantity, optional
+        Initial guess for the y-coordinate of the circle center. Defaults to 0.
+
+    Returns
+    -------
+    radius : astropy.units.Quantity
+        Fitted radius of the circle.
+    center_x : astropy.units.Quantity
+        Fitted x-coordinate of the circle center.
+    center_y : astropy.units.Quantity
+        Fitted y-coordinate of the circle center.
+    radius_err : astropy.units.Quantity
+        Fitted radius of the circle error.
+    center_x_err : astropy.units.Quantity
+        Fitted x-coordinate of the circle center error.
+    center_y_err : astropy.units.Quantity
+        Fitted y-coordinate of the circle center error.
+
+    Raises
+    ------
+    OptionalDependencyMissing
+        If the iminuit package is not installed.
+
+    Notes
+    -----
+    The Taubin circle fit minimizes a specific loss function that balances the
+    squared residuals of the points from the circle with the weights. This method
+    is particularly useful for fitting circles to noisy data.
+
+    References
+    ----------
+    - Barcelona_Muons_TPA_final.pdf (slide 6)
     """
+
     if Minuit is None:
         raise OptionalDependencyMissing("iminuit")
 
@@ -75,52 +149,136 @@ def taubin_circle_fit(x, y, mask):
     x_masked = x[mask]
     y_masked = y[mask]
 
-    R = x.max()  # x.max() just happens to be identical with R in many cases.
+    max_fov = 2 * x.max()
 
-    taubin_r_initial = R / 2
-    taubin_error = R * 0.1
-    xc = 0
-    yc = 0
+    if weights is None:
+        weights_masked = np.ones(len(x_masked))
+    else:
+        weights_masked = weights[mask]
+
+    if original_unit.is_equivalent(u.deg):
+        r_initial = (1.1 * u.deg).to(original_unit) if r_initial is None else r_initial
+    else:
+        r_initial = max_fov / 4.0 * original_unit if r_initial is None else r_initial
+
+    xc_initial = 0 * original_unit if xc_initial is None else xc_initial
+    yc_initial = 0 * original_unit if yc_initial is None else yc_initial
 
     # minimization method
     fit = Minuit(
-        make_taubin_loss_function(x_masked, y_masked), xc=xc, yc=yc, r=taubin_r_initial
+        make_loss_function(x_masked, y_masked, weights_masked),
+        xc=xc_initial.to_value(original_unit),
+        yc=yc_initial.to_value(original_unit),
+        r=r_initial.to_value(original_unit),
     )
     fit.errordef = Minuit.LEAST_SQUARES
 
+    # set initial parameters uncertainty to a big value
+    taubin_error = max_fov * 0.1
     fit.errors["xc"] = taubin_error
     fit.errors["yc"] = taubin_error
     fit.errors["r"] = taubin_error
 
-    fit.limits["xc"] = (-2 * R, 2 * R)
-    fit.limits["yc"] = (-2 * R, 2 * R)
-    fit.limits["r"] = (0, R)
+    # set wide rage for the minimisation
+    fit.limits["xc"] = (-max_fov, max_fov)
+    fit.limits["yc"] = (-max_fov, max_fov)
+    fit.limits["r"] = (0, max_fov)
 
     fit.migrad()
 
     radius = Quantity(fit.values["r"], original_unit)
     center_x = Quantity(fit.values["xc"], original_unit)
     center_y = Quantity(fit.values["yc"], original_unit)
+    radius_err = Quantity(fit.errors["r"], original_unit)
+    center_x_err = Quantity(fit.errors["xc"], original_unit)
+    center_y_err = Quantity(fit.errors["yc"], original_unit)
 
-    return radius, center_x, center_y
+    return radius, center_x, center_y, radius_err, center_x_err, center_y_err
 
 
-def make_taubin_loss_function(x, y):
+def make_loss_function(x, y, w):
     """closure around taubin_loss_function to make
     surviving pixel positions availaboe inside.
 
     x, y: positions of pixels surviving the cleaning
         should not be quantities
+    w : array-like of float, weights for the points
     """
 
     def taubin_loss_function(xc, yc, r):
         """taubin fit formula
         reference : Barcelona_Muons_TPA_final.pdf (slide 6)
         """
-        upper_term = (((x - xc) ** 2 + (y - yc) ** 2 - r**2) ** 2).sum()
 
-        lower_term = ((x - xc) ** 2 + (y - yc) ** 2).sum()
+        distance_squared = (x - xc) ** 2 + (y - yc) ** 2
+        upper_term = ((w * (distance_squared - r**2)) ** 2).sum()
+        lower_term = (w * distance_squared).sum()
 
         return np.abs(upper_term) / np.abs(lower_term)
 
     return taubin_loss_function
+
+
+def naive_circle_fit_error_calculator(x, y, weights, radius, center_x, center_y):
+    """
+    Naive (simplified) error calculator for circular data with weights.
+
+    In this naive approach, we assume zero correlation between the radius,
+    center_x, and center_y. However, the error in the radius is twice as
+    small as that of center_x and center_y, which are assumed to have equal errors.
+
+    Parameters
+    ----------
+    x: array-like or astropy quantity
+        x coordinates of the points
+    y: array-like or astropy quantity
+        y coordinates of the points
+    weights: array-like
+        weights of the points
+    radius : astropy.units.Quantity
+        Fitted radius of the circle.
+    center_x : astropy.units.Quantity
+        Fitted x-coordinate of the circle center.
+    center_y : astropy.units.Quantity
+        Fitted y-coordinate of the circle center.
+
+    Returns
+    -------
+    radius_err : astropy.units.Quantity
+        Fitted radius of the circle error.
+    center_x_err : astropy.units.Quantity
+        Fitted x-coordinate of the circle center error.
+    center_y_err : astropy.units.Quantity
+        Fitted y-coordinate of the circle center error.
+    """
+
+    weights_sum = np.sum(weights)
+    radius_squared = (x - center_x) ** 2 + (y - center_y) ** 2
+    delta = radius_squared - radius**2
+    partial_derivative = np.sqrt(radius_squared + radius**2 / 4)
+
+    delta_weighted_mean = np.average(delta, weights=weights)
+    partial_derivative_weighted_mean = np.average(partial_derivative, weights=weights)
+
+    delta_weighted_variance = np.average(
+        (delta - delta_weighted_mean) ** 2,
+        weights=weights,
+    )
+    partial_derivative_weighted_variance = np.average(
+        (partial_derivative - partial_derivative_weighted_mean) ** 2,
+        weights=weights,
+    )
+
+    parameter_err = (
+        np.sqrt(delta_weighted_variance)
+        / np.sqrt(partial_derivative_weighted_variance)
+        / weights_sum
+        / 2
+        / np.sqrt(2)
+    )
+
+    radius_err = parameter_err / 2
+    center_x_err = parameter_err
+    center_y_err = parameter_err
+
+    return radius_err, center_x_err, center_y_err
diff --git a/src/ctapipe/image/muon/processor.py b/src/ctapipe/image/muon/processor.py
index 4a877b92..a96a8ba4 100644
--- a/src/ctapipe/image/muon/processor.py
+++ b/src/ctapipe/image/muon/processor.py
@@ -182,7 +182,7 @@ class MuonProcessor(TelescopeComponent):
         # 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)
+            ring = self.ring_fitter(geometry, image, mask)
             dist = np.sqrt(
                 (fov_lon - ring.center_fov_lon) ** 2
                 + (fov_lat - ring.center_fov_lat) ** 2
diff --git a/src/ctapipe/image/muon/ring_fitter.py b/src/ctapipe/image/muon/ring_fitter.py
index 8f8d49ee..0c0760c3 100644
--- a/src/ctapipe/image/muon/ring_fitter.py
+++ b/src/ctapipe/image/muon/ring_fitter.py
@@ -22,7 +22,26 @@ def taubin(fov_lon, fov_lat, weights, mask):
     return taubin_circle_fit(fov_lon, fov_lat, mask)
 
 
-FIT_METHOD_BY_NAME = {m.__name__: m for m in [kundu_chaudhuri, taubin]}
+def kundu_chaudhuri_taubin(fov_lon, fov_lat, weights, mask):
+    """taubin_circle_fit with fov_lon, fov_lat, weights, mask interface
+    with initial parameters provided by kundu_chaudhuri"""
+    taubin_r_initial, xc_initial, yc_initial, _, _, _ = kundu_chaudhuri_circle_fit(
+        fov_lon[mask], fov_lat[mask], weights[mask], nan_errors_flag=True
+    )
+    return taubin_circle_fit(
+        fov_lon,
+        fov_lat,
+        mask,
+        weights,
+        taubin_r_initial,
+        xc_initial,
+        yc_initial,
+    )
+
+
+FIT_METHOD_BY_NAME = {
+    m.__name__: m for m in [kundu_chaudhuri, taubin, kundu_chaudhuri_taubin]
+}
 
 __all__ = ["MuonRingFitter"]
 
@@ -32,22 +51,47 @@ class MuonRingFitter(Component):
 
     fit_method = traits.CaselessStrEnum(
         list(FIT_METHOD_BY_NAME.keys()),
-        default_value=list(FIT_METHOD_BY_NAME.keys())[0],
+        default_value="kundu_chaudhuri_taubin",
     ).tag(config=True)
 
-    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")
+    def __call__(self, geom, image, clean_mask):
         """
+        Perform a circle fit to ``image`` with the chosen ``fit_method``.
+
+        Parameters
+        -----------
+        geom: CameraGeometry
+            Defines the pixel coordinates.
+            Must be in the `ctapipe.coordinates.TelescopeFrame`
+        image: np.ndarray[np.float32]
+            Image intensity values
+        clean_mask: np.array([],dtype=bool)
+            Boolean mask of clean pixels, where True means the pixel contains signal.
+            This can be generated using a `ctapipe.image.ImageCleaner`.
+
+        Returns
+        --------
+        MuonRingContainer:
+            Results of the ring fit.
+        """
+
         fit_function = FIT_METHOD_BY_NAME[self.fit_method]
-        radius, center_fov_lon, center_fov_lat = fit_function(
-            fov_lon, fov_lat, img, mask
-        )
+        (
+            radius,
+            center_fov_lon,
+            center_fov_lat,
+            radius_err,
+            center_fov_lon_err,
+            center_fov_lat_err,
+        ) = fit_function(geom.pix_x, geom.pix_y, image, clean_mask)
 
         return MuonRingContainer(
             center_fov_lon=center_fov_lon,
             center_fov_lat=center_fov_lat,
             radius=radius,
+            center_fov_lon_err=center_fov_lon_err,
+            center_fov_lat_err=center_fov_lat_err,
+            radius_err=radius_err,
             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/src/ctapipe/image/toymodel.py b/src/ctapipe/image/toymodel.py
index ae482ae9..ed95d116 100644
--- a/src/ctapipe/image/toymodel.py
+++ b/src/ctapipe/image/toymodel.py
@@ -24,7 +24,9 @@ from abc import ABCMeta, abstractmethod
 import astropy.units as u
 import numpy as np
 from numpy.random import default_rng
+from scipy.integrate import quad
 from scipy.ndimage import convolve1d
+from scipy.special import ellipe
 from scipy.stats import multivariate_normal, norm, skewnorm
 
 from ctapipe.calib.camera.gainselection import GainChannel
@@ -375,20 +377,70 @@ class SkewedGaussian(ImageModel):
 class RingGaussian(ImageModel):
     """A shower image consisting of a ring with gaussian radial profile.
 
-    Simplified model for a muon ring.
+    Simplified model for a muon ring, optionally with asymmetry based on the muon impact.
+
+    Parameters
+    ----------
+    x : u.Quantity
+        x-coordinate of the ring center
+    y : u.Quantity
+        y-coordinate of the ring center
+    radius : u.Quantity
+        ring radius
+    sigma : u.Quantity
+        std. dev. along the radial axis, i.e. the width of point spread function
+    rho : float
+        Ratio of the muon impact parameter to the mirror radius. A value of
+        0 means a homogeneously illuminated ring, values > 1 result in only
+        partially illuminated rings.
+    phi0 : u.Quantity[angle]
+        Position angle of the muon impact on the mirror.
     """
 
-    def __init__(self, x, y, radius, sigma):
+    def __init__(
+        self,
+        x,
+        y,
+        radius,
+        sigma,
+        rho=0.0,
+        phi0=0 * u.deg,
+    ):
         self.unit = x.unit
         self.x = x
         self.y = y
         self.sigma = sigma
         self.radius = radius
-        self.dist = norm(
+        self.r_dist = norm(
             self.radius.to_value(self.unit), self.sigma.to_value(self.unit)
         )
+        self.rho = rho
+        self.phi0 = phi0.to(u.rad)
+
+        if rho >= 1:
+            self.phi_max = np.arcsin(1 / self.rho)
+            integral, _ = quad(self._inner_term, 0, self.phi_max)
+            self.norm = 1.0 / (2.0 * integral)
+        else:
+            self.norm = 0.25 / ellipe(self.rho**2)
 
     def pdf(self, x, y):
         """2d probability for photon electrons in the camera plane."""
-        r = np.sqrt((x - self.x) ** 2 + (y - self.y) ** 2)
-        return self.dist.pdf(r)
+        dx = (x - self.x).to_value(self.unit)
+        dy = (y - self.y).to_value(self.unit)
+        r = np.sqrt(dx**2 + dy**2)
+        phi = np.arctan2(dy, dx)
+        return self.r_dist.pdf(r) * self._pdf_phi(phi)
+
+    def _inner_term(self, phi):
+        return np.sqrt(1 - self.rho**2 * np.sin(phi) ** 2)
+
+    def _pdf_phi(self, phi):
+        phi = phi + self.phi0.to_value(u.rad)
+        if self.rho >= 1:
+            mask = np.abs(phi) < self.phi_max
+            pdf = np.zeros(phi.shape)
+            pdf[mask] = self._inner_term(phi[mask])
+            return self.norm * pdf
+        else:
+            return self.norm * (self._inner_term(phi) + self.rho * np.cos(phi))


Test patch:
diff --git a/src/ctapipe/image/muon/tests/test_muon_fitting.py b/src/ctapipe/image/muon/tests/test_muon_fitting.py
index b3f81994..9b64ea52 100644
--- a/src/ctapipe/image/muon/tests/test_muon_fitting.py
+++ b/src/ctapipe/image/muon/tests/test_muon_fitting.py
@@ -18,7 +18,7 @@ def test_kundu_chaudhuri():
 
         weights = np.ones_like(x)
 
-        fit_radius, fit_x, fit_y = kundu_chaudhuri_circle_fit(x, y, weights)
+        fit_radius, fit_x, fit_y, _, _, _ = kundu_chaudhuri_circle_fit(x, y, weights)
 
         assert np.isclose(fit_x, center_x)
         assert np.isclose(fit_y, center_y)
@@ -37,7 +37,7 @@ def test_kundu_chaudhuri_with_units():
 
     weights = np.ones_like(x)
 
-    fit_radius, fit_x, fit_y = kundu_chaudhuri_circle_fit(x, y, weights)
+    fit_radius, fit_x, fit_y, _, _, _ = kundu_chaudhuri_circle_fit(x, y, weights)
 
     assert fit_x.unit == center_x.unit
     assert fit_y.unit == center_y.unit
diff --git a/src/ctapipe/image/muon/tests/test_ring_fitter.py b/src/ctapipe/image/muon/tests/test_ring_fitter.py
index 7d5b3b25..d5dbdfe3 100644
--- a/src/ctapipe/image/muon/tests/test_ring_fitter.py
+++ b/src/ctapipe/image/muon/tests/test_ring_fitter.py
@@ -1,49 +1,320 @@
+from collections import namedtuple
+
 import astropy.units as u
+import numpy as np
 import pytest
 
+from ctapipe.coordinates import TelescopeFrame
 from ctapipe.image import tailcuts_clean, toymodel
 from ctapipe.image.muon import MuonRingFitter
 
 
-def test_MuonRingFitter_has_methods():
+def test_muon_ring_fitter_has_methods():
     # just to make sure, the test below is running for at least 2 methods
     # basically making sure, we do not test no method at all and always pass
     assert len(MuonRingFitter.fit_method.values) >= 2
 
 
-@pytest.mark.parametrize("method", MuonRingFitter.fit_method.values)
-def test_MuonRingFitter(method, prod5_mst_flashcam):
+parameter_names = [
+    "tel_fixture_name",
+    "method",
+    "center_x",
+    "center_y",
+    "impact_rho",
+    "impact_angle",
+    "intensity",
+    "nsb_level_pe",
+    "picture_thresh",
+    "boundary_thresh",
+]
+Parameters = namedtuple("MuonTestParams", parameter_names)
+
+
+@pytest.mark.parametrize(
+    parameter_names,
+    [
+        Parameters(
+            tel_fixture_name="prod5_lst",
+            method=MuonRingFitter.fit_method.values[0],
+            center_x=-0.3 * u.deg,
+            center_y=0.4 * u.deg,
+            impact_rho=0.5,
+            impact_angle=45 * u.deg,
+            intensity=2000,
+            nsb_level_pe=3,
+            picture_thresh=7,
+            boundary_thresh=5,
+        ),
+        Parameters(
+            tel_fixture_name="prod5_lst",
+            method=MuonRingFitter.fit_method.values[1],
+            center_x=-0.3 * u.deg,
+            center_y=0.4 * u.deg,
+            impact_rho=0.8,
+            impact_angle=45 * u.deg,
+            intensity=3000,
+            nsb_level_pe=3,
+            picture_thresh=7,
+            boundary_thresh=5,
+        ),
+        Parameters(
+            tel_fixture_name="prod5_lst",
+            method=MuonRingFitter.fit_method.values[2],
+            center_x=-0.3 * u.deg,
+            center_y=0.4 * u.deg,
+            impact_rho=0.8,
+            impact_angle=45 * u.deg,
+            intensity=3000,
+            nsb_level_pe=3,
+            picture_thresh=7,
+            boundary_thresh=5,
+        ),
+        Parameters(
+            tel_fixture_name="prod5_mst_flashcam",
+            method=MuonRingFitter.fit_method.values[0],
+            center_x=-1.3 * u.deg,
+            center_y=1.4 * u.deg,
+            impact_rho=0.5,
+            impact_angle=45 * u.deg,
+            intensity=1500,
+            nsb_level_pe=2,
+            picture_thresh=7,
+            boundary_thresh=5,
+        ),
+        Parameters(
+            tel_fixture_name="prod5_mst_flashcam",
+            method=MuonRingFitter.fit_method.values[1],
+            center_x=-1.3 * u.deg,
+            center_y=1.4 * u.deg,
+            impact_rho=0.8,
+            impact_angle=45 * u.deg,
+            intensity=1500,
+            nsb_level_pe=2,
+            picture_thresh=7,
+            boundary_thresh=5,
+        ),
+        Parameters(
+            tel_fixture_name="prod5_mst_flashcam",
+            method=MuonRingFitter.fit_method.values[2],
+            center_x=-1.3 * u.deg,
+            center_y=1.4 * u.deg,
+            impact_rho=0.8,
+            impact_angle=45 * u.deg,
+            intensity=1500,
+            nsb_level_pe=2,
+            picture_thresh=7,
+            boundary_thresh=5,
+        ),
+        Parameters(
+            tel_fixture_name="prod5_mst_nectarcam",
+            method=MuonRingFitter.fit_method.values[0],
+            center_x=-1.3 * u.deg,
+            center_y=1.4 * u.deg,
+            impact_rho=0.5,
+            impact_angle=45 * u.deg,
+            intensity=1500,
+            nsb_level_pe=2,
+            picture_thresh=7,
+            boundary_thresh=5,
+        ),
+        Parameters(
+            tel_fixture_name="prod5_mst_nectarcam",
+            method=MuonRingFitter.fit_method.values[1],
+            center_x=-1.3 * u.deg,
+            center_y=1.4 * u.deg,
+            impact_rho=0.8,
+            impact_angle=45 * u.deg,
+            intensity=1500,
+            nsb_level_pe=2,
+            picture_thresh=7,
+            boundary_thresh=5,
+        ),
+        Parameters(
+            tel_fixture_name="prod5_mst_nectarcam",
+            method=MuonRingFitter.fit_method.values[2],
+            center_x=-1.3 * u.deg,
+            center_y=1.4 * u.deg,
+            impact_rho=0.8,
+            impact_angle=45 * u.deg,
+            intensity=1500,
+            nsb_level_pe=2,
+            picture_thresh=7,
+            boundary_thresh=5,
+        ),
+        Parameters(
+            tel_fixture_name="prod5_sst",
+            method=MuonRingFitter.fit_method.values[0],
+            center_x=-2.5 * u.deg,
+            center_y=1.3 * u.deg,
+            impact_rho=0.3,
+            impact_angle=45 * u.deg,
+            intensity=1200,
+            nsb_level_pe=0,
+            picture_thresh=2,
+            boundary_thresh=1,
+        ),
+        Parameters(
+            tel_fixture_name="prod5_sst",
+            method=MuonRingFitter.fit_method.values[1],
+            center_x=-2.5 * u.deg,
+            center_y=1.3 * u.deg,
+            impact_rho=0.3,
+            impact_angle=45 * u.deg,
+            intensity=1200,
+            nsb_level_pe=0,
+            picture_thresh=2,
+            boundary_thresh=1,
+        ),
+        Parameters(
+            tel_fixture_name="prod5_sst",
+            method=MuonRingFitter.fit_method.values[2],
+            center_x=-2.5 * u.deg,
+            center_y=1.3 * u.deg,
+            impact_rho=0.3,
+            impact_angle=45 * u.deg,
+            intensity=1200,
+            nsb_level_pe=0,
+            picture_thresh=2,
+            boundary_thresh=1,
+        ),
+    ],
+)
+def test_muon_ring_fitter(
+    request,
+    tel_fixture_name,
+    method,
+    center_x,
+    center_y,
+    impact_rho,
+    impact_angle,
+    intensity,
+    nsb_level_pe,
+    picture_thresh,
+    boundary_thresh,
+):
     """test MuonRingFitter"""
     pytest.importorskip("iminuit")
 
-    # flashCam example
-    center_xs = 0.3 * u.m
-    center_ys = 0.6 * u.m
-    radius = 0.3 * u.m
-    width = 0.05 * u.m
+    # Dynamically retrieve the fixture for the specified camera
+    tel = request.getfixturevalue(tel_fixture_name)
+    geom = tel.camera.geometry.transform_to(TelescopeFrame())
+
+    radius = 1.1 * u.deg
+    width = 0.07 * radius
+    min_error = 0.05 * radius
 
     muon_model = toymodel.RingGaussian(
-        x=center_xs,
-        y=center_ys,
+        x=center_x,
+        y=center_y,
         radius=radius,
         sigma=width,
+        rho=impact_rho,
+        phi0=impact_angle,
     )
 
-    # testing with flashcam
-    geom = prod5_mst_flashcam.camera.geometry
     charge, _, _ = muon_model.generate_image(
         geom,
-        intensity=1000,
-        nsb_level_pe=5,
+        intensity=intensity,
+        nsb_level_pe=nsb_level_pe,
     )
-    survivors = tailcuts_clean(geom, charge, 10, 12)
+    survivors = tailcuts_clean(geom, charge, picture_thresh, boundary_thresh)
 
     muonfit = MuonRingFitter(fit_method=method)
-    fit_result = muonfit(geom.pix_x, geom.pix_y, charge, survivors)
+    fit_result = muonfit(geom, charge, survivors)
+
+    assert u.isclose(
+        fit_result.center_fov_lon,
+        center_x,
+        atol=(max(fit_result.center_fov_lon_err, min_error)),
+    )
+    assert u.isclose(
+        fit_result.center_fov_lat,
+        center_y,
+        atol=(max(fit_result.center_fov_lat_err, min_error)),
+    )
+    assert u.isclose(
+        fit_result.radius, radius, atol=(max(fit_result.radius_err, min_error))
+    )
+
+
+parameter_names = [
+    "tel_fixture_name",
+    "method",
+    "center_x",
+    "center_y",
+    "radius",
+    "width",
+]
+Parameters = namedtuple("MuonTestParams", parameter_names)
+
+
+@pytest.mark.parametrize(
+    parameter_names,
+    [
+        Parameters(
+            tel_fixture_name="prod5_lst",
+            method=MuonRingFitter.fit_method.values[0],
+            center_x=-0.3 * u.deg,
+            center_y=0.4 * u.deg,
+            radius=1.1 * u.deg,
+            width=0.1 * u.deg,
+        ),
+        Parameters(
+            tel_fixture_name="prod5_lst",
+            method=MuonRingFitter.fit_method.values[1],
+            center_x=-0.3 * u.deg,
+            center_y=0.4 * u.deg,
+            radius=1.1 * u.deg,
+            width=0.1 * u.deg,
+        ),
+        Parameters(
+            tel_fixture_name="prod5_lst",
+            method=MuonRingFitter.fit_method.values[2],
+            center_x=-0.3 * u.deg,
+            center_y=0.4 * u.deg,
+            radius=1.1 * u.deg,
+            width=0.1 * u.deg,
+        ),
+    ],
+)
+def test_muon_ring_fitter_error_calculator(
+    request,
+    tel_fixture_name,
+    method,
+    center_x,
+    center_y,
+    radius,
+    width,
+):
+    """test MuonRingFitter error_calculator"""
+    pytest.importorskip("iminuit")
+
+    # Dynamically retrieve the fixture for the specified camera
+    tel = request.getfixturevalue(tel_fixture_name)
+    geom = tel.camera.geometry.transform_to(TelescopeFrame())
+
+    muon_model = toymodel.RingGaussian(
+        x=center_x,
+        y=center_y,
+        radius=radius,
+        sigma=width,
+        rho=0.5,
+        phi0=40 * u.deg,
+    )
+
+    charge, _, _ = muon_model.generate_image(
+        geom,
+        intensity=3000,
+        nsb_level_pe=0,
+    )
+    survivors = tailcuts_clean(geom, charge, 7, 5)
 
-    print(fit_result)
-    print(center_xs, center_ys, radius)
+    muonfit = MuonRingFitter(fit_method=method)
+    fit_result = muonfit(geom, charge, survivors)
 
-    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)
+    assert np.isfinite(fit_result.center_fov_lon_err)
+    assert np.isfinite(fit_result.center_fov_lat_err)
+    assert np.isfinite(fit_result.radius_err)
+    assert fit_result.center_fov_lon_err > 0
+    assert fit_result.center_fov_lat_err > 0
+    assert fit_result.radius_err > 0