merged main

python/CMakeLists.txt

@@ -39,7 +39,10 @@ set(PYTHON_FILES
     aare/ClusterFinder.py
     aare/ClusterVector.py
     aare/Cluster.py
+    aare/FastPedestal.py
+    aare/Pedestal.py
     aare/calibration.py
+    aare/factory.py
     aare/experimental.py
     aare/RawFile.py
     aare/transform.py

python/aare/Cluster.py

@@ -1,6 +1,7 @@
+# SPDX-License-Identifier: MPL-2.0
 from . import _aare 
 import numpy as np
-from .ClusterFinder import _type_to_char
+from .factory import _type_to_char
 
 
 def Cluster(x : int, y : int, data, cluster_size=(3,3), dtype = np.int32):
@@ -21,4 +22,4 @@ def Cluster(x : int, y : int, data, cluster_size=(3,3), dtype = np.int32):
     except AttributeError:
         raise ValueError(f"Unsupported combination of type and cluster size: {dtype}/{cluster_size} when requesting {class_name}")
    
-    return cls(x, y, data) 
+    return cls(x, y, data)

python/aare/ClusterFinder.py

@@ -1,21 +1,10 @@
 # SPDX-License-Identifier: MPL-2.0
 from . import _aare
 import numpy as np
+from .factory import _type_to_char
 
 _supported_cluster_sizes = [(2,2), (3,3), (5,5), (7,7), (9,9),]
 
-def _type_to_char(dtype):
-    if dtype == np.int32:
-        return 'i'
-    elif dtype == np.float32:
-        return 'f'
-    elif dtype == np.float64:
-        return 'd'
-    elif dtype == np.int16:
-        return 'i16'
-    else:
-        raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32, np.float32, and np.float64 are supported.")
-
 def _get_class(name, cluster_size, dtype):
     """
     Helper function to get the class based on the name, cluster size, and dtype.
@@ -29,24 +18,55 @@ def _get_class(name, cluster_size, dtype):
 
 
 
-def ClusterFinder(image_size, cluster_size=(3,3), n_sigma=5, dtype = np.int32, capacity = 1024):
+def ClusterFinder(image_size, *, cluster_size=(3,3), n_sigma=5, dtype = np.int32, capacity = 1024, min_pedestal_samples = 1000):
     """
     Factory function to create a ClusterFinder object. Provides a cleaner syntax for 
     the templated ClusterFinder in C++.
+
+    Parameters
+    ----------
+    image_size : tuple
+        The size of the image as a tuple (height, width).
+    cluster_size : tuple, optional
+        The size of the cluster to find as a tuple (height, width). Default is (3,3).
+    n_sigma : int, optional
+        Multiplier of the standard deviation used as a threshold to identify potential photon pixels. Default is 5.
+    dtype : data-type, optional
+        The data type of the image. Default is np.int32.
+    capacity : int, optional
+        The maximum number of clusters than can be stored before reallocating. Default is 1024.
+    min_pedestal_samples : int, optional
+        The minimum number of pedestal samples to accumulate before using the pedestal. Default is 1000.
     """
     cls = _get_class("ClusterFinder", cluster_size, dtype)
-    return cls(image_size, n_sigma=n_sigma, capacity=capacity)
+    return cls(image_size, n_sigma=n_sigma, capacity=capacity, min_pedestal_samples=min_pedestal_samples)
 
 
-
-def ClusterFinderMT(image_size, cluster_size = (3,3), dtype=np.int32, n_sigma=5, capacity = 1024, n_threads = 3): 
+def ClusterFinderMT(image_size, *, cluster_size = (3,3), dtype=np.int32, n_sigma=5, capacity = 1024, n_threads = 3, min_pedestal_samples = 1000): 
     """ 
     Factory function to create a ClusterFinderMT object. Provides a cleaner syntax for 
     the templated ClusterFinderMT in C++.
+
+    Parameters
+    ----------
+    image_size : tuple
+        The size of the image as a tuple (height, width).
+    cluster_size : tuple, optional
+        The size of the cluster to find as a tuple (height, width). Default is (3,3).
+    n_sigma : int, optional
+        Multiplier of the standard deviation used as a threshold to identify potential photon pixels. Default is 5.
+    dtype : data-type, optional
+        The data type of the image. Default is np.int32.
+    capacity : int, optional
+        The maximum number of clusters than can be stored before reallocating. Default is 1024.
+    n_threads : int, optional
+        The number of threads to use for processing. Default is 3.
+    min_pedestal_samples : int, optional
+        The minimum number of pedestal samples to accumulate before using the pedestal. Default is 1000.
     """
 
     cls = _get_class("ClusterFinderMT", cluster_size, dtype)
-    return cls(image_size, n_sigma=n_sigma, capacity=capacity, n_threads=n_threads)
+    return cls(image_size, n_sigma=n_sigma, capacity=capacity, min_pedestal_samples=min_pedestal_samples, n_threads=n_threads)
 
 
 def ClusterCollector(clusterfindermt, dtype=np.int32): 
@@ -69,19 +89,55 @@ def ClusterFileSink(clusterfindermt, cluster_file, dtype=np.int32):
 
 
 def ClusterFile(fname, cluster_size=(3,3), dtype=np.int32, chunk_size = 1000, mode = "r"):
-    """
-    Factory function to create a ClusterFile object. Provides a cleaner syntax for
-    the templated ClusterFile in C++.
+    """Create a reader or writer for a binary cluster file.
+
+    Parameters
+    ----------
+    fname : path-like
+        Cluster file to open.
+    cluster_size : tuple[int, int], default=(3, 3)
+        Cluster dimensions stored in the file.
+    dtype : numpy dtype, default=numpy.int32
+        Data type of the cluster values stored in the file.
+    chunk_size : int, default=1000
+        Maximum number of selected clusters returned by ``chunks()`` and
+        default iteration. Must be positive when iterating over chunks.
+    mode : {"r", "w", "a"}, default="r"
+        Open for reading, truncate and write, or append, respectively.
+
+    Returns
+    -------
+    ClusterFile
+        The compiled ClusterFile specialization matching ``cluster_size`` and
+        ``dtype``.
+
+    Notes
+    -----
+    The file format contains no cluster shape or data-type metadata. Supplying
+    values that do not match the file causes its bytes to be interpreted
+    incorrectly. Use ``frames()`` to iterate over complete frames, including
+    empty or fully filtered frames with their stored frame numbers. Use
+    ``chunks()`` or ``chunks(chunk_size)`` to iterate over selected clusters
+    in batches. Chunks may combine frames, so their frame number is not
+    reliable per-cluster metadata.
+
+    Iterators consume the current file position without rewinding; use one
+    traversal at a time. Each result owns its storage and remains valid after
+    advancing the iterator or closing the file.
+
+    Examples
+    --------
 
     .. code-block:: python
 
         from aare import ClusterFile
-        
-        with ClusterFile("clusters.clust", cluster_size=(3,3), dtype=np.int32) as cf:
-            # cf is now a ClusterFile_Cluster3x3i object but you don't need to know that.
-            for clusters in cf:
-                # Loop over clusters in chunks of 1000 
-                # The type of clusters will be a ClusterVector_Cluster3x3i in this case
+
+        with ClusterFile(
+            "clusters.clust", cluster_size=(3, 3), dtype=np.int32
+        ) as cf:
+            for clusters in cf.chunks():
+                # Process clusters in chunks of at most 1000.
+                ...
 
     """
 

python/aare/ClusterVector.py

@@ -7,16 +7,33 @@ from .ClusterFinder import _get_class
 
 def ClusterVector(cluster_size=(3,3), dtype = np.int32):
     """
-    Factory function to create a ClusterVector object. Provides a cleaner syntax for
-    the templated ClusterVector in C++.
+    Create an empty ClusterVector for a supported cluster size and pixel dtype.
 
-    .. code-block:: python
+    Parameters
+    ----------
+    cluster_size : tuple[int, int], default=(3, 3)
+        Cluster dimensions in the x and y directions.
+    dtype : numpy.dtype, default=numpy.int32
+        Pixel storage type. The cluster size and dtype combination must have a
+        compiled binding.
 
-        from aare import ClusterVector
-        
-        ClusterVector(cluster_size=(3,3), dtype=np.float64)
+    Returns
+    -------
+    ClusterVector
+        An empty vector with frame number 0 and space reserved for at least
+        1024 clusters.
+
+    Raises
+    ------
+    ValueError
+        If the requested size and dtype combination is unavailable.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from aare import ClusterVector
+    >>> clusters = ClusterVector(cluster_size=(3, 3), dtype=np.float64)
     """
 
     cls = _get_class("ClusterVector", cluster_size, dtype)
     return cls()
-

python/aare/FastPedestal.py

@@ -0,0 +1,56 @@
+# SPDX-License-Identifier: MPL-2.0
+
+import numpy as np
+
+from .factory import _get_typed_class
+
+
+def _get_fast_pedestal_class(dtype):
+    return _get_typed_class("FastPedestal", dtype)
+
+
+def FastPedestal(rows, cols, n_samples=1000, dtype=np.float64):
+    """Create an empty per-pixel running pedestal.
+
+    This factory hides the dtype suffix used by the templated C++ bindings.
+    Call ``add_init_frame()`` exactly ``n_samples`` times before using the
+    statistics or calling ``push_ema()``. Subsequent frames have weight
+    ``1 / n_samples`` in the running mean and population variance.
+
+    Args:
+        rows: Number of image rows.
+        cols: Number of image columns.
+        n_samples: Initialization frame count and steady-state update-weight
+            denominator.
+        dtype: Output dtype for the mean and standard deviation.
+            Supported values are ``np.float64``, ``np.float32``, and
+            ``np.int16``.
+    """
+    cls = _get_fast_pedestal_class(dtype)
+    return cls(rows, cols, n_samples)
+
+
+def from_file(filename, n_samples=1000, skip_first=0, dtype=np.float64):
+    """Create a FastPedestal from frames in a file.
+
+    After ignoring ``skip_first`` frames, the next ``n_samples`` frames
+    initialize the pedestal. Every remaining frame is then applied as a
+    steady-state update. Input frames are read as uint16 data.
+
+    Args:
+        filename: Input image file.
+        n_samples: Number of frames used for initialization.
+        skip_first: Number of leading frames to ignore.
+        dtype: Output dtype for the mean and standard deviation.
+
+    Raises:
+        RuntimeError: If fewer than ``n_samples`` frames remain after
+            ``skip_first`` or ``n_samples`` is zero.
+    """
+    cls = _get_fast_pedestal_class(dtype)
+    return cls.from_file(
+        filename, n_samples=n_samples, skip_first=skip_first
+    )
+
+
+FastPedestal.from_file = from_file

python/aare/Pedestal.py

@@ -0,0 +1,26 @@
+# SPDX-License-Identifier: MPL-2.0
+
+import numpy as np
+
+from .factory import _get_typed_class
+
+
+def Pedestal(rows, cols, n_samples=1000, dtype=np.float64):
+    """Create an empty per-pixel running pedestal.
+
+    This factory hides the dtype suffix used by the templated C++ bindings.
+    Call ``push()`` to update the statistics and cached mean for each frame.
+    Statistics are available during initialization and are zero for empty pixels.
+    Internal moments and variance always use double precision.
+
+    Args:
+        rows: Number of image rows.
+        cols: Number of image columns.
+        n_samples: Number of samples accumulated before switching to
+            steady-state updates with weight ``1 / n_samples``.
+        dtype: Output dtype for the mean and standard deviation.
+            Supported values are ``np.float64``, ``np.float32``, and
+            ``np.int16``.
+    """
+    cls = _get_typed_class("Pedestal", dtype)
+    return cls(rows, cols, n_samples)

python/aare/RawFile.py

@@ -43,7 +43,7 @@ class RawFile(_aare.RawFile):
         Returns:
             int: Number of frames in file.
         """
-        return super().frames_in_file
+        return self.total_frames
 
     def __enter__(self):
         return self

python/aare/__init__.py

@@ -5,6 +5,16 @@ from . import _aare
 from . import transform
 from . import experimental
 
+from ._aare import (
+    FastPedestal_d,
+    FastPedestal_f,
+    FastPedestal_i16,
+    Pedestal_d,
+    Pedestal_f,
+    Pedestal_i16,
+    ClusterFinder_Cluster3x3i,
+    VarClusterFinder,
+)
 from ._aare import (
     File,
     JungfrauDataFile,
@@ -12,20 +22,25 @@ from ._aare import (
     RawSubFile,
 )
 from ._aare import Pedestal_d, Pedestal_f, ClusterFinder_Cluster3x3i, VarClusterFinder
-from ._aare import DetectorType, ReadoutMode 
+from ._aare import DetectorType, FrameDiscardPolicy, ReadoutMode
 from ._aare import hitmap
 from ._aare import ROI
 from ._aare import corner 
 
+from ._aare import UDPPortPosition
+
 # from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i
 
 from ._version import __version__
+from .FastPedestal import FastPedestal
+from .Pedestal import Pedestal
 from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, ClusterFileSink, ClusterFile
 from .ClusterVector import ClusterVector
 from .Cluster import Cluster
 
 from ._aare import Gaussian, RisingScurve, FallingScurve, Pol1, Pol2, GaussianErfcPlateau, GaussianChargeSharing, GaussianChargeSharingKb
 from ._aare import fit
+from ._aare import Minimizer
 from ._aare import Interpolator
 from ._aare import calculate_eta2, calculate_eta3, calculate_cross_eta3, calculate_full_eta2
 from ._aare import reduce_to_2x2, reduce_to_3x3

python/aare/_version.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 # This file is used to get the version of the package from the VERSION file
 from importlib.metadata import PackageNotFoundError, version
 from pathlib import Path

python/aare/calibration.py

@@ -5,9 +5,12 @@ def load_calibration(fname, hg0=False):
     """
     Load calibration data from a file.
     
-    Parameters:
-    fname (str): Path to the calibration file.
-    hg0 (bool): If True, load HG0 calibration data instead of G0.
+    Parameters
+    ----------
+    fname : str
+        Path to the calibration file.
+    hg0 : bool
+        If True, load HG0 calibration data instead of G0.
 
     """
     gains = 3

python/aare/factory.py

@@ -0,0 +1,36 @@
+# SPDX-License-Identifier: MPL-2.0
+
+import numpy as np
+
+from . import _aare
+
+
+_TYPE_TO_CHAR = {
+    np.dtype(np.int32): "i",
+    np.dtype(np.float32): "f",
+    np.dtype(np.float64): "d",
+    np.dtype(np.int16): "i16",
+}
+
+
+def _type_to_char(dtype):
+    """Return the suffix used by bindings instantiated for ``dtype``."""
+    try:
+        return _TYPE_TO_CHAR[np.dtype(dtype)]
+    except (KeyError, TypeError):
+        supported = ", ".join(str(dtype) for dtype in _TYPE_TO_CHAR)
+        raise ValueError(
+            f"Unsupported dtype: {dtype}. Supported dtypes are {supported}."
+        ) from None
+
+
+def _get_typed_class(name, dtype):
+    """Return a bound class named ``<name>_<dtype suffix>``."""
+    class_name = f"{name}_{_type_to_char(dtype)}"
+    try:
+        return getattr(_aare, class_name)
+    except AttributeError:
+        raise ValueError(
+            f"Unsupported dtype for {name}: {dtype} "
+            f"(binding {class_name} is not available)."
+        ) from None

python/src/bind_Cluster.hpp

@@ -10,7 +10,6 @@
 #include <pybind11/stl_bind.h>
 
 namespace py = pybind11;
-using pd_type = double;
 
 using namespace aare;
 
@@ -25,20 +24,26 @@ void define_Cluster(py::module &m, const std::string &typestr) {
     py::class_<Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>>(
         m, class_name.c_str(), py::buffer_protocol())
 
+        // conversion ok: small input, accept lists and any numeric dtype
         .def(py::init([](CoordType x, CoordType y,
                          py::array_t<Type, py::array::forcecast> data) {
-            py::buffer_info buf_info = data.request();
+            constexpr py::ssize_t n_pixels = ClusterSizeX * ClusterSizeY;
+            if (data.ndim() != 1 || data.size() != n_pixels) {
+                throw py::value_error(fmt::format(
+                    "Cluster data must be a 1D array of {} values, got {}D "
+                    "array with {} values",
+                    n_pixels, data.ndim(), data.size()));
+            }
             Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> cluster;
             cluster.x = x;
             cluster.y = y;
             auto r = data.template unchecked<1>(); // no bounds checks
-            for (py::ssize_t i = 0; i < data.size(); ++i) {
+            for (py::ssize_t i = 0; i < n_pixels; ++i) {
                 cluster.data[i] = r(i);
             }
             return cluster;
         }))
 
-        // TODO! Review if to keep or not
         .def_property_readonly(
             "data",
             [](Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> &c)
@@ -61,6 +66,26 @@ void define_Cluster(py::module &m, const std::string &typestr) {
         .def_readonly("y",
                       &Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>::y)
 
+        .def("__repr__",
+             [class_name](const Cluster<Type, ClusterSizeX, ClusterSizeY,
+                                        CoordType> &self) {
+                 fmt::memory_buffer data;
+                 fmt::format_to(std::back_inserter(data), "[[");
+                 for (size_t i = 0; i < self.data.size(); ++i) {
+                     if (i != 0) {
+                         fmt::format_to(std::back_inserter(data),
+                                        i % ClusterSizeX == 0 ? "], [" : ", ");
+                     }
+                     fmt::format_to(std::back_inserter(data), "{}",
+                                    self.data[i]);
+                 }
+                 fmt::format_to(std::back_inserter(data), "]]");
+
+                 return fmt::format("{}(x={}, y={}, data={})", class_name,
+                                    self.x, self.y,
+                                    fmt::string_view(data.data(), data.size()));
+             })
+
         .def(
             "max_sum_2x2",
             [](Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> &self) {
@@ -68,8 +93,16 @@ void define_Cluster(py::module &m, const std::string &typestr) {
                 return py::make_tuple(max_sum.sum,
                                       static_cast<int>(max_sum.index));
             },
-            R"(calculates sum of 2x2 subcluster with highest energy and index relative to cluster center 0: top_left, 1: top_right, 2: bottom_left, 3: bottom_right
-            )");
+            R"doc(
+            Return the highest-sum center-adjacent 2x2 subcluster.
+
+            Returns
+            -------
+            tuple
+                ``(sum, index)``, where index is 0 for top-left, 1 for
+                top-right, 2 for bottom-left, or 3 for bottom-right relative to
+                the cluster center.
+            )doc");
 }
 
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
@@ -81,7 +114,14 @@ void reduce_to_3x3(py::module &m) {
         [](const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
             return reduce_to_3x3(cl);
         },
-        py::return_value_policy::move, R"(Reduce cluster to 3x3 subcluster)");
+        py::return_value_policy::move, py::arg("cluster"), R"doc(
+        Return the 3x3 block around the cluster's center index.
+
+        Both input dimensions must be at least 3, and at least one must be
+        greater than 3.
+        The input coordinates are preserved and output data is stored in
+        row-major order.
+        )doc");
 }
 
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
@@ -93,16 +133,12 @@ void reduce_to_2x2(py::module &m) {
         [](const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
             return reduce_to_2x2(cl);
         },
-        py::return_value_policy::move,
-        R"(
-        Reduce cluster to 2x2 subcluster by taking the 2x2 subcluster with 
-        the highest photon energy.
+        py::return_value_policy::move, py::arg("cluster"), R"doc(
+        Return the highest-sum center-adjacent 2x2 block.
 
-        RETURN: 
-        
-        reduced cluster (cluster is filled in row major ordering starting at the top left. Thus for a max subcluster in the top left corner the photon hit is at the fourth position.)
-        
-        )");
+        The input coordinates are preserved and output data is stored in
+        row-major order.
+        )doc");
 }
 
-#pragma GCC diagnostic pop
+#pragma GCC diagnostic pop

python/src/bind_ClusterCollector.hpp

@@ -6,6 +6,8 @@
 #include "aare/ClusterVector.hpp"
 #include "aare/NDView.hpp"
 #include "aare/Pedestal.hpp"
+
+#include "module_config.hpp"
 #include "np_helper.hpp"
 
 #include <cstdint>
@@ -15,7 +17,6 @@
 #include <pybind11/stl_bind.h>
 
 namespace py = pybind11;
-using pd_type = double;
 
 using namespace aare;
 
@@ -30,7 +31,7 @@ void define_ClusterCollector(py::module &m, const std::string &typestr) {
     using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
 
     py::class_<ClusterCollector<ClusterType>>(m, class_name.c_str())
-        .def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *>())
+        .def(py::init<ClusterFinderMT<ClusterType, uint16_t, pd_type> *>())
         .def("stop", &ClusterCollector<ClusterType>::stop)
         .def(
             "steal_clusters",

python/src/bind_ClusterFile.hpp

@@ -11,6 +11,7 @@
 #include <pybind11/stl.h>
 #include <pybind11/stl/filesystem.h>
 #include <string>
+#include <utility>
 
 // Disable warnings for unused parameters, as we ignore some
 // in the __exit__ method
@@ -20,18 +21,90 @@
 namespace py = pybind11;
 using namespace ::aare;
 
+template <typename Range> class ClusterFileIterator {
+    Range m_range;
+    std::optional<typename Range::Iterator> m_iterator;
+    bool m_done{false};
+
+  public:
+    explicit ClusterFileIterator(Range range) : m_range(std::move(range)) {}
+
+    typename Range::Iterator::value_type next() {
+        if (m_done) {
+            throw py::stop_iteration();
+        }
+        // Advance before returning the next result, never after yielding it.
+        if (m_iterator) {
+            ++*m_iterator;
+        } else {
+            m_iterator.emplace(m_range.begin());
+        }
+        if (*m_iterator == m_range.end()) {
+            m_done = true;
+            throw py::stop_iteration();
+        }
+        return std::move(**m_iterator);
+    }
+};
+
+template <typename Iterator>
+void define_cluster_file_iterator(py::module &m, const std::string &name) {
+    py::class_<Iterator>(m, name.c_str())
+        .def(
+            "__iter__", [](Iterator &self) -> Iterator & { return self; },
+            py::return_value_policy::reference_internal)
+        .def("__next__", &Iterator::next);
+}
+
 template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY,
           typename CoordType = uint16_t>
 void define_ClusterFile(py::module &m, const std::string &typestr) {
 
     using ClusterType = Cluster<Type, CoordSizeX, CoordSizeY, CoordType>;
+    using File = ClusterFile<ClusterType>;
+    using FrameIterator = ClusterFileIterator<typename File::FrameRange>;
+    using ChunkIterator = ClusterFileIterator<typename File::ChunkRange>;
 
     auto class_name = fmt::format("ClusterFile_{}", typestr);
+    define_cluster_file_iterator<FrameIterator>(m,
+                                                class_name + "_FrameIterator");
+    define_cluster_file_iterator<ChunkIterator>(m,
+                                                class_name + "_ChunkIterator");
 
-    py::class_<ClusterFile<ClusterType>>(m, class_name.c_str())
+    py::class_<ClusterFile<ClusterType>>(
+        m, class_name.c_str(),
+        "Read and write legacy binary cluster files. The format contains no "
+        "cluster type or shape metadata, so this class must match the file.")
         .def(py::init<const std::filesystem::path &, size_t,
                       const std::string &>(),
-             py::arg(), py::arg("chunk_size") = 1000, py::arg("mode") = "r")
+             py::arg("fname"), py::arg("chunk_size") = 1000,
+             py::arg("mode") = "r",
+             "Open a cluster file. Mode must be 'r' to read, 'w' to truncate "
+             "and write, or 'a' to append.")
+        .def(
+            "frames", [](File &self) { return FrameIterator(self.frames()); },
+            py::keep_alive<0, 1>(),
+            "Iterate over complete frames from the current position, including "
+            "empty and fully filtered frames with their stored frame numbers. "
+            "Results own their storage. The iterator keeps the file alive but "
+            "shares its cursor; use only one traversal at a time.")
+        .def(
+            "chunks", [](File &self) { return ChunkIterator(self.chunks()); },
+            py::keep_alive<0, 1>(),
+            "Iterate from the current position using the constructor's chunk "
+            "size, which must be positive. Results own their storage. The "
+            "iterator keeps the file alive but shares its cursor; use only "
+            "one traversal at a time.")
+        .def(
+            "chunks",
+            [](File &self, size_t chunk_size) {
+                return ChunkIterator(self.chunks(chunk_size));
+            },
+            py::arg("chunk_size"), py::keep_alive<0, 1>(),
+            "Iterate over chunks of up to chunk_size selected clusters. The "
+            "size must be positive and does not change the constructor's "
+            "default. Chunks may span frames; their frame numbers are not "
+            "per-cluster metadata. Only the final chunk can be short.")
         .def(
             "read_clusters",
             [](ClusterFile<ClusterType> &self, size_t n_clusters) {
@@ -39,32 +112,62 @@ void define_ClusterFile(py::module &m, const std::string &typestr) {
                     self.read_clusters(n_clusters));
                 return v;
             },
-            py::return_value_policy::take_ownership, py::arg("n_clusters"))
-        .def("read_frame",
-             [](ClusterFile<ClusterType> &self) {
-                 auto v = new ClusterVector<ClusterType>(self.read_frame());
-                 return v;
-             })
-        .def("set_roi", &ClusterFile<ClusterType>::set_roi, py::arg("roi"))
-        .def("tell", &ClusterFile<ClusterType>::tell)
-        .def("estimate_n_clusters",
-             &ClusterFile<ClusterType>::estimate_n_clusters)
+            py::return_value_policy::take_ownership, py::arg("n_clusters"),
+            "Read up to n_clusters without preserving frame boundaries. The "
+            "result may combine frames, so its frame number is not reliable "
+            "per-cluster metadata.")
         .def(
-            "set_noise_map",
-            [](ClusterFile<ClusterType> &self, py::array_t<int32_t> noise_map) {
+            "read_frame",
+            [](ClusterFile<ClusterType> &self) -> py::object {
+                auto clusters = self.read_frame();
+                if (!clusters) {
+                    return py::none();
+                }
+                return py::cast(
+                    new ClusterVector<ClusterType>(std::move(*clusters)),
+                    py::return_value_policy::take_ownership);
+            },
+            "Read and return the next complete frame with its frame number, "
+            "or None at end of file.")
+        .def("set_roi", &ClusterFile<ClusterType>::set_roi, py::arg("roi"),
+             "Select clusters whose centers lie within the half-open ROI.")
+        .def("tell", &ClusterFile<ClusterType>::tell,
+             "Return the current byte position in the file.")
+        .def("estimate_n_clusters",
+             &ClusterFile<ClusterType>::estimate_n_clusters,
+             "Estimate the number of clusters from the file size. Frame "
+             "headers can make this larger than the actual count.")
+        .def(
+            "set_noise_map", // conversions ok
+            [](ClusterFile<ClusterType> &self,
+               py::array_t<int32_t, py::array::c_style | py::array::forcecast>
+                   noise_map) {
                 auto view = make_view_2d(noise_map);
                 self.set_noise_map(view);
             },
-            py::arg("noise_map"))
+            py::arg("noise_map"),
+            "Set a two-dimensional, C-contiguous int32 noise map indexed as "
+            "[y, x]. The map must cover every cluster center coordinate.")
 
-        .def("set_gain_map",
-             [](ClusterFile<ClusterType> &self, py::array_t<double> gain_map) {
-                 auto view = make_view_2d(gain_map);
-                 self.set_gain_map(view);
-             })
+        .def(
+            "set_gain_map", // conversions ok
+            [](ClusterFile<ClusterType> &self,
+               py::array_t<double, py::array::c_style | py::array::forcecast>
+                   gain_map) {
+                auto view = make_view_2d(gain_map);
+                self.set_gain_map(view);
+            },
+            py::arg("gain_map"),
+            "Set a two-dimensional, C-contiguous float64 gain map in "
+            "ADU/energy, indexed as [y, x]. Clusters whose complete footprint "
+            "extends beyond the map are retained with all data values set to "
+            "zero.")
 
-        .def("close", &ClusterFile<ClusterType>::close)
-        .def("write_frame", &ClusterFile<ClusterType>::write_frame)
+        .def("close", &ClusterFile<ClusterType>::close,
+             "Close the file. Calling close more than once is safe.")
+        .def("write_frame", &ClusterFile<ClusterType>::write_frame,
+             py::arg("clusters"),
+             "Write one ClusterVector, including its frame number.")
         .def("__enter__", [](ClusterFile<ClusterType> &self) { return &self; })
         .def("__exit__",
              [](ClusterFile<ClusterType> &self,
@@ -74,14 +177,8 @@ void define_ClusterFile(py::module &m, const std::string &typestr) {
                  self.close();
              })
         .def("__iter__", [](ClusterFile<ClusterType> &self) { return &self; })
-        .def("__next__", [](ClusterFile<ClusterType> &self) {
-            auto v = new ClusterVector<ClusterType>(
-                self.read_clusters(self.chunk_size()));
-            if (v->size() == 0) {
-                throw py::stop_iteration();
-            }
-            return v;
-        });
+        .def("__next__",
+             [](File &self) { return ChunkIterator(self.chunks()).next(); });
 }
 
 #pragma GCC diagnostic pop

python/src/bind_ClusterFileSink.hpp

@@ -15,8 +15,6 @@
 #include <pybind11/stl_bind.h>
 
 namespace py = pybind11;
-using pd_type = double;
-
 using namespace aare;
 
 #pragma GCC diagnostic push
@@ -29,6 +27,8 @@ void define_ClusterFileSink(py::module &m, const std::string &typestr) {
 
     using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
 
+    // TODO! adapt to set pedestal type (needs templating of ClusterFileSink)
+    // or maybe access through base class?
     py::class_<ClusterFileSink<ClusterType>>(m, class_name.c_str())
         .def(py::init<ClusterFinderMT<ClusterType, uint16_t, double> *,
                       const std::filesystem::path &>())

python/src/bind_ClusterFinder.hpp

@@ -6,6 +6,8 @@
 #include "aare/ClusterVector.hpp"
 #include "aare/NDView.hpp"
 #include "aare/Pedestal.hpp"
+
+#include "module_config.hpp"
 #include "np_helper.hpp"
 
 #include <cstdint>
@@ -15,7 +17,6 @@
 #include <pybind11/stl_bind.h>
 
 namespace py = pybind11;
-using pd_type = double;
 
 using namespace aare;
 
@@ -31,8 +32,10 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) {
 
     py::class_<ClusterFinder<ClusterType, uint16_t, pd_type>>(
         m, class_name.c_str())
-        .def(py::init<Shape<2>, pd_type, size_t>(), py::arg("image_size"),
-             py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000)
+        .def(py::init<Shape<2>, pd_type, size_t, size_t>(),
+             py::arg("image_size"), py::arg("n_sigma") = 5.0,
+             py::arg("capacity") = 1'000'000,
+             py::arg("min_pedestal_samples") = 1000)
 
         .def_property(
             "nSigma",
@@ -40,14 +43,22 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) {
             &ClusterFinder<ClusterType, uint16_t, pd_type>::set_nSigma,
             R"(number of sigma above the pedestal to consider a photon during cluster finding.)")
 
-        .def("push_pedestal_frame",
-             [](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
-                py::array_t<uint16_t> frame) {
-                 auto view = make_view_2d(frame);
-                 self.push_pedestal_frame(view);
-             })
+        .def(
+            "push_pedestal_frame",
+            [](ClusterFinder<ClusterType, uint16_t, pd_type> &self,
+               py::array_t<uint16_t> frame) {
+                auto view = make_view_2d(frame);
+                self.push_pedestal_frame(view);
+            },
+            py::arg("frame").noconvert())
         .def("clear_pedestal",
              &ClusterFinder<ClusterType, uint16_t, pd_type>::clear_pedestal)
+        .def("update_threshold",
+             &ClusterFinder<ClusterType, uint16_t, pd_type>::update_threshold)
+        .def_property_readonly(
+            "pedestal_ready",
+            &ClusterFinder<ClusterType, uint16_t, pd_type>::pedestal_ready,
+            R"(True once enough pedestal frames have been pushed to find clusters.)")
         .def_property_readonly(
             "pedestal",
             [](ClusterFinder<ClusterType, uint16_t, pd_type> &self) {
@@ -79,7 +90,7 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) {
                 self.find_clusters(view, frame_number);
                 return;
             },
-            py::arg(), py::arg("frame_number") = 0);
+            py::arg("frame").noconvert(), py::arg("frame_number") = 0);
 }
 
 #pragma GCC diagnostic pop

python/src/bind_ClusterFinderMT.hpp

@@ -6,6 +6,8 @@
 #include "aare/ClusterVector.hpp"
 #include "aare/NDView.hpp"
 #include "aare/Pedestal.hpp"
+
+#include "module_config.hpp"
 #include "np_helper.hpp"
 
 #include <cstdint>
@@ -15,7 +17,6 @@
 #include <pybind11/stl_bind.h>
 
 namespace py = pybind11;
-using pd_type = double;
 
 using namespace aare;
 
@@ -31,24 +32,27 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) {
 
     py::class_<ClusterFinderMT<ClusterType, uint16_t, pd_type>>(
         m, class_name.c_str())
-        .def(py::init<Shape<2>, pd_type, size_t, size_t>(),
+        .def(py::init<Shape<2>, pd_type, size_t, size_t, size_t, size_t>(),
              py::arg("image_size"), py::arg("n_sigma") = 5.0,
-             py::arg("capacity") = 2048, py::arg("n_threads") = 3)
-        .def("push_pedestal_frame",
-             [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
-                py::array_t<uint16_t> frame) {
-                 auto view = make_view_2d(frame);
-                 self.push_pedestal_frame(view);
-             })
+             py::arg("capacity") = 2048, py::arg("n_threads") = 3,
+             py::arg("queue_depth") = 16,
+             py::arg("min_pedestal_samples") = 1000)
+        .def(
+            "push_pedestal_frame",
+            [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
+               py::array_t<uint16_t> frame) {
+                auto view = make_view_2d(frame);
+                self.push_pedestal_frame(view);
+            },
+            py::arg("frame").noconvert())
         .def(
             "find_clusters",
             [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self,
                py::array_t<uint16_t> frame, uint64_t frame_number) {
                 auto view = make_view_2d(frame);
                 self.find_clusters(view, frame_number);
-                return;
             },
-            py::arg(), py::arg("frame_number") = 0)
+            py::arg("frame").noconvert(), py::arg("frame_number") = 0)
         .def_property_readonly(
             "cluster_size",
             [](ClusterFinderMT<ClusterType, uint16_t, pd_type> &self) {
@@ -56,6 +60,12 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) {
             })
         .def("clear_pedestal",
              &ClusterFinderMT<ClusterType, uint16_t, pd_type>::clear_pedestal)
+        .def("update_threshold",
+             &ClusterFinderMT<ClusterType, uint16_t, pd_type>::update_threshold)
+        .def(
+            "pedestal_ready",
+            &ClusterFinderMT<ClusterType, uint16_t, pd_type>::pedestal_ready,
+            R"(waits for the queued frames, then returns True if all cluster finders have enough pedestal frames to find clusters.)")
         .def("sync", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::sync)
         .def("stop", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::stop)
         .def("start", &ClusterFinderMT<ClusterType, uint16_t, pd_type>::start)

python/src/bind_ClusterVector.hpp

@@ -15,7 +15,6 @@
 #include <pybind11/stl_bind.h>
 
 namespace py = pybind11;
-using pd_type = double;
 
 using namespace aare;
 
@@ -30,37 +29,72 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
 
     py::class_<ClusterVector<
         Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>, void>>(
-        m, class_name.c_str(),
+        m, class_name.c_str(), R"doc(
+        A contiguous, move-only container of fixed-size clusters.
+
+        The class supports the Python buffer protocol, so ``numpy.array`` can
+        either copy its data or create a zero-copy view. A zero-copy view is
+        valid only while the ClusterVector's underlying allocation and size
+        remain unchanged.
+        )doc",
         py::buffer_protocol())
 
-        .def(py::init()) // TODO change!!!
+        .def(py::init(), R"doc(
+        Create an empty ClusterVector with frame number 0 and space reserved
+        for at least 1024 clusters.
+        )doc")
 
         .def(
             "__call__",
             [](ClusterVector<ClusterType> &self, py::array_t<bool> mask) {
                 return self(make_view_1d(mask));
             },
-            py::arg("mask"), R"(
-            Create a copy of the clustervector and apply a boolean mask to the ClusterVector.
+            py::arg("mask").noconvert(), R"doc(
+            Return a filtered copy of this ClusterVector.
 
             Parameters
             ----------
+            mask : numpy.ndarray
+                One-dimensional, writable, C-contiguous array with dtype
+                ``numpy.bool_`` and one element per cluster.
 
-            mask : 1d boolean numpy array
-                Mask to apply to the ClusterVector. Must be the same length as the number of clusters in the ClusterVector. 
+            Returns
+            -------
+            ClusterVector
+                Selected clusters in their original order. The frame number is
+                preserved.
+            )doc")
 
-            )")
+        .def(
+            "push_back",
+            [](ClusterVector<ClusterType> &self, const ClusterType &cluster) {
+                self.push_back(cluster);
+            },
+            py::arg("cluster"), R"doc(
+            Append one cluster.
 
-        .def("push_back",
-             [](ClusterVector<ClusterType> &self, const ClusterType &cluster) {
-                 self.push_back(cluster);
-             })
+            Notes
+            -----
+            Do not call this method while a zero-copy NumPy view of the
+            ClusterVector exists. Reallocation invalidates the view, while an
+            append without reallocation leaves its shape unchanged.
+            )doc")
 
-        .def("sum",
-             [](ClusterVector<ClusterType> &self) {
-                 auto *vec = new std::vector<Type>(self.sum());
-                 return return_vector(vec);
-             })
+        .def(
+            "sum",
+            [](ClusterVector<ClusterType> &self) {
+                auto *vec = new std::vector<Type>(self.sum());
+                return return_vector(vec);
+            },
+            R"doc(
+            Return the sum of all pixels in each cluster.
+
+            Returns
+            -------
+            numpy.ndarray
+                One value per cluster, in container order and with the cluster
+                pixel dtype.
+            )doc")
         .def(
             "sum_2x2",
             [](ClusterVector<ClusterType> &self) {
@@ -69,24 +103,45 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
 
                 return return_vector(vec);
             },
-            R"(calculates sum of 2x2 subcluster with highest energy and index relative to cluster center 0: top_left, 1: top_right, 2: bottom_left, 3: bottom_right
-            )")
-        .def_property_readonly("size", &ClusterVector<ClusterType>::size)
-        .def("empty", &ClusterVector<ClusterType>::empty)
-        .def("item_size", &ClusterVector<ClusterType>::item_size)
-        .def_property_readonly("fmt",
-                               [typestr](ClusterVector<ClusterType> &self) {
-                                   return fmt_format<ClusterType>;
-                               })
+            R"doc(
+            Return the highest-sum center-adjacent 2x2 subcluster for each
+            cluster.
+
+            Returns
+            -------
+            numpy.ndarray
+                Structured array with ``sum`` and ``index`` fields. Indices are
+                0 for top-left, 1 for top-right, 2 for bottom-left, and 3 for
+                bottom-right, relative to the cluster center.
+            )doc")
+        .def_property_readonly("size", &ClusterVector<ClusterType>::size,
+                               "Number of stored clusters.")
+        .def("empty", &ClusterVector<ClusterType>::empty,
+             "Return True when no clusters are stored.")
+        .def("item_size", &ClusterVector<ClusterType>::item_size,
+             "Return the size in bytes of one stored cluster, including "
+             "padding.")
+        .def_property_readonly(
+            "fmt",
+            [typestr](ClusterVector<ClusterType> &self) {
+                return fmt_format<ClusterType>;
+            },
+            "PEP 3118 format string for one stored cluster.")
 
         .def_property_readonly("cluster_size_x",
-                               &ClusterVector<ClusterType>::cluster_size_x)
+                               &ClusterVector<ClusterType>::cluster_size_x,
+                               "Cluster size in the x dimension.")
         .def_property_readonly("cluster_size_y",
-                               &ClusterVector<ClusterType>::cluster_size_y)
+                               &ClusterVector<ClusterType>::cluster_size_y,
+                               "Cluster size in the y dimension.")
         .def_property_readonly("capacity",
-                               &ClusterVector<ClusterType>::capacity)
+                               &ClusterVector<ClusterType>::capacity,
+                               "Number of clusters that fit without "
+                               "reallocation.")
         .def_property("frame_number", &ClusterVector<ClusterType>::frame_number,
-                      &ClusterVector<ClusterType>::set_frame_number)
+                      &ClusterVector<ClusterType>::set_frame_number,
+                      "Signed 32-bit frame number; 0 can indicate clusters "
+                      "from multiple frames.")
         .def_buffer(
             [typestr](ClusterVector<ClusterType> &self) -> py::buffer_info {
                 return py::buffer_info(
@@ -100,31 +155,38 @@ void define_ClusterVector(py::module &m, const std::string &typestr) {
             });
 
     // Free functions using ClusterVector
-    m.def("hitmap",
-          [](std::array<size_t, 2> image_size, ClusterVector<ClusterType> &cv) {
-              // Create a numpy array to hold the hitmap
-              // The shape of the array is (image_size[0], image_size[1])
-              // note that the python array is passed as [row, col] which
-              // is the opposite of the clusters [x,y]
-              py::array_t<int32_t> hitmap(image_size);
-              auto r = hitmap.mutable_unchecked<2>();
+    m.def(
+        "hitmap",
+        [](std::array<size_t, 2> image_size, ClusterVector<ClusterType> &cv) {
+            // Create a numpy array to hold the hitmap
+            // The shape of the array is (image_size[0], image_size[1])
+            // note that the python array is passed as [row, col] which
+            // is the opposite of the clusters [x,y]
+            py::array_t<int32_t> hitmap(image_size);
+            auto r = hitmap.mutable_unchecked<2>();
 
-              // Initialize hitmap to 0
-              for (py::ssize_t i = 0; i < r.shape(0); i++)
-                  for (py::ssize_t j = 0; j < r.shape(1); j++)
-                      r(i, j) = 0;
+            // Initialize hitmap to 0
+            for (py::ssize_t i = 0; i < r.shape(0); i++)
+                for (py::ssize_t j = 0; j < r.shape(1); j++)
+                    r(i, j) = 0;
 
-              // Loop over the clusters and increment the hitmap
-              // Skip out of bound clusters
-              for (const auto &cluster : cv) {
-                  auto x = cluster.x;
-                  auto y = cluster.y;
-                  if (x < image_size[1] && y < image_size[0])
-                      r(cluster.y, cluster.x) += 1;
-              }
+            // Loop over the clusters and increment the hitmap
+            // Skip out of bound clusters
+            for (const auto &cluster : cv) {
+                auto x = cluster.x;
+                auto y = cluster.y;
+                if (x < image_size[1] && y < image_size[0])
+                    r(cluster.y, cluster.x) += 1;
+            }
 
-              return hitmap;
-          });
+            return hitmap;
+        },
+        py::arg("image_size"), py::arg("clusters"), R"doc(
+        Count cluster centers into an ``int32`` image whose shape is given by
+        ``image_size`` as ``(rows, columns)``. Element ``[y, x]`` contains the
+        number of cluster centers at that coordinate. Out-of-bounds centers are
+        ignored.
+        )doc");
 }
 
 template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
@@ -137,15 +199,12 @@ void define_2x2_reduction(py::module &m) {
             return new ClusterVector<Cluster<Type, 2, 2, CoordType>>(
                 reduce_to_2x2(cv));
         },
-        R"(
-        Reduce cluster to 2x2 subcluster by taking the 2x2 subcluster with 
-        the highest photon energy.
-
-        Parameters
-
-        cv : ClusterVector (clusters are filled in row-major ordering starting at the top left. Thus for a max subcluster in the top left corner the photon hit is at the fourth position.)
-
-        )",
+        R"doc(
+        Reduce every cluster to its highest-sum center-adjacent 2x2 block.
+        Returns a new ClusterVector; cluster order, coordinates, and frame
+        number are preserved. Input pixel data is interpreted in row-major
+        order.
+        )doc",
         py::arg("clustervector"));
 }
 
@@ -160,13 +219,13 @@ void define_3x3_reduction(py::module &m) {
             return new ClusterVector<Cluster<Type, 3, 3, CoordType>>(
                 reduce_to_3x3(cv));
         },
-        R"(
-        Reduce cluster to 3x3 subcluster 
-
-        Parameters
-        
-        cv : ClusterVector   
-        )",
+        R"doc(
+        Reduce every cluster to the 3x3 block around its center index.
+        Both input dimensions must be at least 3, and at least one must be
+        greater than 3.
+        Returns a new ClusterVector; cluster order, coordinates, and frame
+        number are preserved.
+        )doc",
         py::arg("clustervector"));
 }
 

python/src/bind_Defs.hpp

@@ -1,3 +1,4 @@
+// SPDX-License-Identifier: MPL-2.0
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 
@@ -31,4 +32,76 @@ void define_defs_bindings(py::module &m) {
     auto jungfrau = py::class_<Jungfrau>(m, "Jungfrau");
     jungfrau.attr("nRows") = Jungfrau::nRows;
     jungfrau.attr("nCols") = Jungfrau::nCols;
+    py::class_<ROI>(m, "ROI")
+        .def(py::init<>())
+        .def(py::init<ssize_t, ssize_t, ssize_t, ssize_t>(), py::arg("xmin"),
+             py::arg("xmax"), py::arg("ymin"), py::arg("ymax"))
+        .def_readwrite("xmin", &ROI::xmin)
+        .def_readwrite("xmax", &ROI::xmax)
+        .def_readwrite("ymin", &ROI::ymin)
+        .def_readwrite("ymax", &ROI::ymax)
+
+        .def(
+            "slice",
+            [](const ROI &self) {
+                return std::make_tuple(py::slice(self.ymin, self.ymax, 1),
+                                       py::slice(self.xmin, self.xmax, 1));
+            },
+            R"doc(
+            The slice can be used to index into a 2D array to extract the region of interest.
+
+            Returns
+            -------
+            tuple of slice
+                A tuple containing two slice objects for the x and y dimensions.
+        )doc")
+
+        .def("__str__",
+             [](const ROI &self) {
+                 return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}",
+                                    self.xmin, self.xmax, self.ymin, self.ymax);
+             })
+        .def("__repr__",
+             [](const ROI &self) {
+                 return fmt::format(
+                     "<ROI: xmin: {} xmax: {} ymin: {} ymax: {}>", self.xmin,
+                     self.xmax, self.ymin, self.ymax);
+             })
+        .def("__iter__",
+             [](const ROI &self) {
+                 return py::make_iterator(&self.xmin, &self.ymax + 1); // NOLINT
+             })
+
+        .def("__eq__", [](const ROI &self, const ROI &other) {
+            return self.xmin == other.xmin && self.xmax == other.xmax &&
+                   self.ymin == other.ymin && self.ymax == other.ymax;
+        });
+
+    py::enum_<DetectorType>(m, "DetectorType")
+        .value("Jungfrau", DetectorType::Jungfrau)
+        .value("Eiger", DetectorType::Eiger)
+        .value("Mythen3", DetectorType::Mythen3)
+        .value("Moench", DetectorType::Moench)
+        .value("Moench03", DetectorType::Moench03)
+        .value("Moench03_old", DetectorType::Moench03_old)
+        .value("ChipTestBoard", DetectorType::ChipTestBoard)
+        .value("Unknown", DetectorType::Unknown);
+
+    py::enum_<UDPPortPosition>(m, "UDPPortPosition")
+        .value("LEFT", UDPPortPosition::LEFT)
+        .value("RIGHT", UDPPortPosition::RIGHT)
+        .value("TOP", UDPPortPosition::TOP)
+        .value("BOTTOM", UDPPortPosition::BOTTOM);
+
+    py::enum_<TimingMode>(m, "TimingMode")
+        .value("AUTO_TIMING", TimingMode::AUTO_TIMING)
+        .value("TRIGGER_EXPOSURE", TimingMode::TRIGGER_EXPOSURE)
+        .value("GATED", TimingMode::GATED)
+        .value("BURST_TRIGGER", TimingMode::BURST_TRIGGER)
+        .value("TRIGGER_GATED", TimingMode::TRIGGER_GATED);
+
+    py::enum_<FrameDiscardPolicy>(m, "FrameDiscardPolicy")
+        .value("NoDiscard", FrameDiscardPolicy::NoDiscard)
+        .value("Discard", FrameDiscardPolicy::Discard)
+        .value("DiscardPartial", FrameDiscardPolicy::DiscardPartial);
 }

python/src/bind_Eta.hpp

@@ -1,6 +1,9 @@
+// SPDX-License-Identifier: MPL-2.0
 #include "aare/CalculateEta.hpp"
+#include "aare/to_string.hpp"
 
 #include <cstdint>
+#include <fmt/format.h>
 // #include <pybind11/native_enum.h> only for version 3
 #include <pybind11/pybind11.h>
 
@@ -18,7 +21,11 @@ void define_eta(py::module &m, const std::string &typestr) {
         .def_readwrite("c", &Eta2<T>::c,
                        "eta corner value cTopLeft, cTopRight, "
                        "cBottomLeft, cBottomRight")
-        .def_readwrite("sum", &Eta2<T>::sum, "photon energy of cluster");
+        .def_readwrite("sum", &Eta2<T>::sum, "photon energy of cluster")
+        .def("__repr__", [class_name](const Eta2<T> &self) {
+            return fmt::format("{}(x={}, y={}, c={}, sum={})", class_name,
+                               self.x, self.y, to_string(self.c), self.sum);
+        });
 }
 
 void define_corner_enum(py::module &m) {
@@ -101,4 +108,4 @@ void register_calculate_3x3eta(py::module &m) {
         },
         R"(calculates eta3x3 taking into account cross pixels in cluster)",
         py::arg("cluster"));
-}
+}

python/src/bind_FastPedestal.hpp

@@ -0,0 +1,135 @@
+// SPDX-License-Identifier: MPL-2.0
+
+#include "aare/FastPedestal.hpp"
+#include "np_helper.hpp"
+
+#include <cstdint>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+
+template <typename SUM_TYPE>
+void define_fast_pedestal_bindings(py::module &m, const std::string &name) {
+
+    py::class_<FastPedestal<SUM_TYPE>>(
+        m, name.c_str(),
+        "Maintain a per-pixel running mean and population standard deviation.",
+        py::buffer_protocol())
+        .def(py::init<uint32_t, uint32_t, uint32_t>(), py::arg("rows"),
+             py::arg("cols"), py::arg("n_samples"),
+             "Construct an empty pedestal. It becomes ready after n_samples "
+             "calls to add_init_frame().")
+        .def(py::init<uint32_t, uint32_t>(), py::arg("rows"), py::arg("cols"),
+             "Construct an empty pedestal with n_samples=1000.")
+
+        .def(
+            "mean",
+            [](FastPedestal<SUM_TYPE> &self) {
+                auto mean = new NDArray<SUM_TYPE, 2>{};
+                *mean = self.mean();
+                return return_image_data(mean);
+            },
+            "Return a copy of the cached mean. The pedestal must be ready.")
+        .def(
+            "std",
+            [](FastPedestal<SUM_TYPE> &self) {
+                auto standard_deviation = new NDArray<SUM_TYPE, 2>{};
+                *standard_deviation = self.std();
+                return return_image_data(standard_deviation);
+            },
+            "Return the population standard deviation as a NumPy array. The "
+            "pedestal must be ready.")
+        .def(
+            "view",
+            [](py::object self_py) {
+                return py::module_::import("numpy").attr("asarray")(self_py);
+            },
+            "Return a non-owning, non-writable NumPy view of the cached mean. "
+            "The pedestal must be ready.")
+
+        // We need to buffer protocol to allow for numpy operations using the
+        // pedestal mean
+        .def_buffer([](FastPedestal<SUM_TYPE> &self) {
+            auto mean = self.view();
+            return py::buffer_info(
+                const_cast<SUM_TYPE *>(mean.data()), sizeof(SUM_TYPE),
+                py::format_descriptor<SUM_TYPE>::format(), 2,
+                {static_cast<py::ssize_t>(mean.shape(0)),
+                 static_cast<py::ssize_t>(mean.shape(1))},
+                {static_cast<py::ssize_t>(mean.strides()[0] * sizeof(SUM_TYPE)),
+                 static_cast<py::ssize_t>(mean.strides()[1] *
+                                          sizeof(SUM_TYPE))},
+                true);
+        })
+        // Subtracting a FastPedestal from a NumPy array
+        .def(
+            "__array_ufunc__",
+            [](py::object self, py::object ufunc, const std::string &method,
+               py::args inputs, py::kwargs kwargs) -> py::object {
+                if (method != "__call__" || inputs.size() != 2 ||
+                    inputs[1].ptr() != self.ptr() ||
+                    py::cast<std::string>(ufunc.attr("__name__")) !=
+                        "subtract") {
+                    return py::reinterpret_borrow<py::object>(
+                        Py_NotImplemented);
+                }
+
+                auto mean =
+                    py::module_::import("builtins").attr("memoryview")(self);
+                return ufunc(inputs[0], mean, **kwargs);
+            },
+            "Support subtracting a FastPedestal from a NumPy array.")
+        .def("clear", py::overload_cast<>(&FastPedestal<SUM_TYPE>::clear),
+             "Reset all statistics and initialization state to zero.")
+        .def_property_readonly("rows", &FastPedestal<SUM_TYPE>::rows,
+                               "Number of image rows.")
+        .def_property_readonly("cols", &FastPedestal<SUM_TYPE>::cols,
+                               "Number of image columns.")
+        .def_property_readonly("cur_samples",
+                               &FastPedestal<SUM_TYPE>::cur_samples,
+                               "Number of initialization frames accumulated. "
+                               "Steady-state pushes do not change it.")
+        .def_property_readonly(
+            "ready", &FastPedestal<SUM_TYPE>::ready,
+            "Whether n_samples initialization frames have been accumulated.")
+        .def_property_readonly(
+            "n_samples", &FastPedestal<SUM_TYPE>::n_samples,
+            "Initialization frame count and steady-state update-weight "
+            "denominator.")
+        .def(
+            "clone",
+            [](FastPedestal<SUM_TYPE> &pedestal) {
+                return FastPedestal<SUM_TYPE>(pedestal);
+            },
+            "Return an independent copy of the pedestal and its state.")
+        .def(
+            "push_ema",
+            [](FastPedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> frame) {
+                pedestal.push_ema(make_view_2d(frame));
+            },
+            py::arg("frame").noconvert(),
+            "Update exponential moving average. The pedstal must  "
+            "already be ready for this update.")
+        .def(
+            "add_init_frame",
+            [](FastPedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> frame) {
+                pedestal.add_init_frame(make_view_2d(frame));
+            },
+            py::arg("frame").noconvert(),
+            "Accumulate one uint16 initialization frame. Call exactly "
+            "n_samples times to make the pedestal ready.")
+        .def_static(
+            "from_file",
+            [](const std::filesystem::path &filename, uint32_t n_samples,
+               uint32_t skip_first) {
+                return FastPedestal<SUM_TYPE>::template from_file<uint16_t>(
+                    filename, n_samples, skip_first);
+            },
+            py::arg("filename"), py::arg("n_samples") = 1000,
+            py::arg("skip_first") = 0,
+            "Create a pedestal from a uint16 file. Skip skip_first frames, use "
+            "the next n_samples for initialization, then apply every remaining "
+            "frame as a steady-state update.");
+}

python/src/bind_Interpolator.hpp

@@ -101,9 +101,7 @@ void define_interpolation_bindings(py::module &m) {
     auto interpolator =
         py::class_<aare::Interpolator>(m, "Interpolator")
             .def(py::init(
-                [](py::array_t<double,
-                               py::array::c_style | py::array::forcecast>
-                       etacube,
+                [](py::array_t<double> etacube,
                    py::array_t<double> xbins, py::array_t<double> ybins,
                    py::array_t<double> ebins) {
                     return Interpolator(
@@ -124,9 +122,9 @@ void define_interpolation_bindings(py::module &m) {
                 ebins: 
                     bin edges of photon energy
                 )doc",
-                py::arg("etacube"),
-                py::arg("xbins"), py::arg("ybins"),
-                py::arg("ebins"))
+                py::arg("etacube").noconvert(),
+                py::arg("xbins").noconvert(), py::arg("ybins").noconvert(),
+                py::arg("ebins").noconvert())
 
             .def(py::init(
                 [](py::array_t<double> xbins, py::array_t<double> ybins,
@@ -146,13 +144,12 @@ void define_interpolation_bindings(py::module &m) {
                     bin edges of etay
                 ebins: 
                     bin edges of photon energy
-                )", py::arg("xbins"),
-                py::arg("ybins"), py::arg("ebins"))
+                )", py::arg("xbins").noconvert(),
+                py::arg("ybins").noconvert(), py::arg("ebins").noconvert())
             .def(
                 "rosenblatttransform",
                 [](Interpolator &self,
-                   py::array_t<double,
-                               py::array::c_style | py::array::forcecast>
+                   py::array_t<double>
                        etacube) {
                     return self.rosenblatttransform(make_view_3d(etacube));
                 },
@@ -162,7 +159,7 @@ void define_interpolation_bindings(py::module &m) {
                 etacube: 
                     joint distribution of eta_x, eta_y and photon energy (**Note:** for the joint distribution first dimension is eta_x, second: eta_y, third: energy bins.)
                 )",
-                py::arg("etacube"))
+                py::arg("etacube").noconvert())
             .def("get_ietax",
                  [](Interpolator &self) {
                      auto *ptr = new NDArray<double, 3>{};

python/src/bind_Pedestal.hpp

@@ -0,0 +1,121 @@
+// SPDX-License-Identifier: MPL-2.0
+
+#include "aare/Pedestal.hpp"
+#include "np_helper.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+
+template <typename SUM_TYPE>
+void define_pedestal_bindings(py::module &m, const std::string &name) {
+
+    py::class_<Pedestal<SUM_TYPE>>(
+        m, name.c_str(),
+        "Maintain a per-pixel running mean and population standard deviation. "
+        "Statistics are available during initialization and are zero for empty "
+        "pixels.",
+        py::buffer_protocol())
+        .def(py::init<uint32_t, uint32_t, uint32_t>(), py::arg("rows"),
+             py::arg("cols"), py::arg("n_samples"),
+             "Construct an empty pedestal. Each pixel accumulates n_samples "
+             "values before switching to exponential updates.")
+        .def(py::init<uint32_t, uint32_t>(), py::arg("rows"), py::arg("cols"),
+             "Construct an empty pedestal with n_samples=1000.")
+        .def(
+            "mean",
+            [](Pedestal<SUM_TYPE> &self) {
+                auto mea = new NDArray<SUM_TYPE, 2>{};
+                *mea = self.mean();
+                return return_image_data(mea);
+            },
+            "Return a copy of the cached mean. Empty pixels return zero.")
+        .def(
+            "view",
+            [](py::object self_py) {
+                return py::module_::import("numpy").attr("asarray")(self_py);
+            },
+            "Return a non-owning, non-writable NumPy view of the cached mean.")
+        .def(
+            "std",
+            [](Pedestal<SUM_TYPE> &self) {
+                auto std = new NDArray<SUM_TYPE, 2>{};
+                *std = self.std();
+                return return_image_data(std);
+            },
+            "Return the population standard deviation as a NumPy array. "
+            "Empty pixels return zero.")
+        .def(
+            "__array_ufunc__",
+            [](py::object self, py::object ufunc, const std::string &method,
+               py::args inputs, py::kwargs kwargs) -> py::object {
+                if (method != "__call__" || inputs.size() != 2 ||
+                    inputs[1].ptr() != self.ptr() ||
+                    py::cast<std::string>(ufunc.attr("__name__")) !=
+                        "subtract") {
+                    return py::reinterpret_borrow<py::object>(
+                        Py_NotImplemented);
+                }
+
+                auto mean =
+                    py::module_::import("builtins").attr("memoryview")(self);
+                return ufunc(inputs[0], mean, **kwargs);
+            },
+            "Support subtracting a Pedestal from a NumPy array.")
+        .def("clear", py::overload_cast<>(&Pedestal<SUM_TYPE>::clear),
+             "Reset all statistics and per-pixel sample counts to zero.")
+        .def_property_readonly("rows", &Pedestal<SUM_TYPE>::rows,
+                               "Number of image rows.")
+        .def_property_readonly("cols", &Pedestal<SUM_TYPE>::cols,
+                               "Number of image columns.")
+        .def_property_readonly(
+            "n_samples", &Pedestal<SUM_TYPE>::n_samples,
+            "Initialization sample count per pixel and steady-state "
+            "update-weight denominator.")
+        .def(
+            "clone",
+            [&](Pedestal<SUM_TYPE> &pedestal) {
+                return Pedestal<SUM_TYPE>(pedestal);
+            },
+            "Return an independent copy of the pedestal and its state.")
+        // TODO! add push for other data types
+        .def(
+            "push",
+            [](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> f) {
+                auto v = make_view_2d(f);
+                pedestal.push(v);
+            },
+            py::arg("frame").noconvert(),
+            "Accumulate or exponentially update every pixel from a "
+            "C-contiguous uint16 frame matching the pedestal shape. After "
+            "n_samples values per pixel, new values have weight 1 / n_samples.")
+        .def(
+            "push_with_threshold",
+            [](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> f,
+               py::array_t<SUM_TYPE> threshold) {
+                auto frame_view = make_view_2d(f);
+                auto threshold_view = make_view_2d(threshold);
+                pedestal.push_with_threshold(frame_view, threshold_view);
+            },
+            py::arg("frame").noconvert(), py::arg("threshold").noconvert(),
+            "Push only pixels where abs(frame - mean) is strictly less than "
+            "threshold. Both arrays must be C-contiguous with the pedestal "
+            "shape; frame must be uint16 and threshold must use the output "
+            "dtype. Rejected pixels keep their statistics and sample counts.")
+        .def_buffer([](Pedestal<SUM_TYPE> &self) {
+            auto mean = self.view();
+            return py::buffer_info(
+                const_cast<SUM_TYPE *>(mean.data()), sizeof(SUM_TYPE),
+                py::format_descriptor<SUM_TYPE>::format(), 2,
+                {static_cast<py::ssize_t>(mean.shape(0)),
+                 static_cast<py::ssize_t>(mean.shape(1))},
+                {static_cast<py::ssize_t>(mean.strides()[0] * sizeof(SUM_TYPE)),
+                 static_cast<py::ssize_t>(mean.strides()[1] *
+                                          sizeof(SUM_TYPE))},
+                true);
+        });
+}

python/src/bind_PedestalTrackingPixelHistogram.hpp

@@ -68,8 +68,7 @@ void define_pedestal_tracking_pixel_histogram_bindings(py::module &m) {
         .def(
             "push_pedestal_no_update",
             [](PedestalTrackingPixelHistogram &self,
-               py::array_t<PedestalTrackingPixelHistogram::FrameType, 0>
-                   frame) {
+               py::array_t<PedestalTrackingPixelHistogram::FrameType> frame) {
                 auto view = make_view_2d(frame);
                 self.push_pedestal_no_update(view);
             },
@@ -146,8 +145,7 @@ void define_pedestal_tracking_pixel_histogram_bindings(py::module &m) {
         .def(
             "fill_async",
             [](PedestalTrackingPixelHistogram &self,
-               py::array_t<PedestalTrackingPixelHistogram::FrameType, 0>
-                   image) {
+               py::array_t<PedestalTrackingPixelHistogram::FrameType> image) {
                 // Copy the numpy buffer into an owned NDArray while we
                 // still hold the GIL so we don't depend on the array's
                 // backing storage outliving this call.
@@ -195,9 +193,13 @@ void define_pedestal_tracking_pixel_histogram_bindings(py::module &m) {
              Args:
                  file_path: Path to the file to fill from
                  max_frames: Maximum number of frames to fill from the file (default: -1)
+                 reader_threads: Number of parallel file reader workers (default: 2)
+                 reader_chunk_size: Frames claimed by each reader worker per batch (default: 4)
              )",
              py::call_guard<py::gil_scoped_release>(), py::arg("fname"),
-             py::arg("max_frames") = -1, py::arg("verbose") = false)
+             py::arg("max_frames") = -1, py::arg("verbose") = false,
+             py::arg("reader_threads") = std::size_t{2},
+             py::arg("reader_chunk_size") = std::size_t{4})
         .def("process_pedestal_file",
              &PedestalTrackingPixelHistogram::process_pedestal_file,
              R"(

python/src/bind_PixelHistogram.hpp

@@ -46,7 +46,7 @@ void define_pixel_histogram_binding(py::module &m, const char *class_name,
 
         .def(
             "fill_async",
-            [](Hist &self, py::array_t<double, 0> image) {
+            [](Hist &self, py::array_t<double> image) {
                 // Copy the numpy buffer into an owned NDArray while we
                 // still hold the GIL so we don't depend on the array's
                 // backing storage outliving this call.

python/src/bind_RawFile.hpp

@@ -11,6 +11,7 @@
 
 #include <cstdint>
 #include <filesystem>
+#include <fmt/format.h>
 #include <pybind11/iostream.h>
 #include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
@@ -21,14 +22,23 @@
 namespace py = pybind11;
 using namespace ::aare;
 
+namespace {
+std::runtime_error raw_file_read_error(RawFile &file,
+                                       const std::string &message) {
+    return std::runtime_error(fmt::format(
+        "Error reading frame index {} from file '{}': {}", file.tell(),
+        file.master().master_fname().string(), message));
+}
+} // namespace
+
 void define_raw_file_io_bindings(py::module &m) {
     py::class_<RawFile>(m, "RawFile")
         .def(py::init<const std::filesystem::path &>())
         .def("read_frame",
              [](RawFile &self) {
                  if (self.n_modules_in_roi().size() > 1) {
-                     throw std::runtime_error(
-                         "File contains multiple ROIs - use read_ROIs()");
+                     throw raw_file_read_error(
+                         self, "File contains multiple ROIs - use read_ROIs()");
                  }
 
                  std::vector<size_t> shape;
@@ -51,16 +61,19 @@ void define_raw_file_io_bindings(py::module &m) {
             "read_n",
             [](RawFile &self, size_t n_frames) {
                 if (self.n_modules_in_roi().size() > 1) {
-                    throw std::runtime_error(
+                    throw raw_file_read_error(
+                        self,
                         "File contains multiple ROIs - use read_n_ROIs() to "
                         "read a specific ROI or use read_ROIs and "
                         "read one frame at a time.");
                 }
                 // adjust for actual frames left in the file
                 n_frames =
-                    std::min(n_frames, self.total_frames() - self.tell());
+                    std::min(n_frames, self.tell() < self.total_frames()
+                                           ? self.total_frames() - self.tell()
+                                           : size_t{0});
                 if (n_frames == 0) {
-                    throw std::runtime_error("No frames left in file");
+                    throw raw_file_read_error(self, "No frames left in file");
                 }
                 std::vector<size_t> shape{n_frames, self.rows(), self.cols()};
 
@@ -89,12 +102,13 @@ void define_raw_file_io_bindings(py::module &m) {
             "read_roi",
             [](RawFile &self, const size_t roi_index) {
                 if (self.num_rois() == 0) {
-                    throw std::runtime_error(LOCATION + "No ROIs defined.");
+                    throw raw_file_read_error(self,
+                                              LOCATION + "No ROIs defined.");
                 }
 
                 if (roi_index >= self.num_rois()) {
-                    throw std::runtime_error(LOCATION +
-                                             "ROI index out of range.");
+                    throw raw_file_read_error(
+                        self, LOCATION + "ROI index out of range.");
                 }
 
                 // return headers from all subfiles
@@ -141,7 +155,8 @@ void define_raw_file_io_bindings(py::module &m) {
             "read_rois",
             [](RawFile &self) {
                 if (self.num_rois() == 0) {
-                    throw std::runtime_error(LOCATION + "No ROIs defined.");
+                    throw raw_file_read_error(self,
+                                              LOCATION + "No ROIs defined.");
                 }
 
                 size_t number_of_ROIs = self.num_rois();
@@ -196,19 +211,22 @@ void define_raw_file_io_bindings(py::module &m) {
             "read_n_with_roi",
             [](RawFile &self, const size_t num_frames, const size_t roi_index) {
                 if (self.num_rois() == 0) {
-                    throw std::runtime_error(LOCATION + "No ROIs defined.");
+                    throw raw_file_read_error(self,
+                                              LOCATION + "No ROIs defined.");
                 }
 
                 if (roi_index >= self.num_rois()) {
-                    throw std::runtime_error(LOCATION +
-                                             "ROI index out of range.");
+                    throw raw_file_read_error(
+                        self, LOCATION + "ROI index out of range.");
                 }
 
                 // adjust for actual frames left in the file
                 size_t n_frames =
-                    std::min(num_frames, self.total_frames() - self.tell());
+                    std::min(num_frames, self.tell() < self.total_frames()
+                                             ? self.total_frames() - self.tell()
+                                             : size_t{0});
                 if (n_frames == 0) {
-                    throw std::runtime_error("No frames left in file");
+                    throw raw_file_read_error(self, "No frames left in file");
                 }
                 std::vector<size_t> shape{
                     n_frames, self.roi_geometries(roi_index).pixels_y(),
@@ -306,4 +324,4 @@ void define_raw_file_io_bindings(py::module &m) {
         .def_property_readonly("n_modules", &RawFile::n_modules)
         .def_property_readonly("n_modules_in_roi", &RawFile::n_modules_in_roi)
         .def_property_readonly("num_rois", &RawFile::num_rois);
-}
+}

python/src/bind_calibration.hpp

@@ -11,38 +11,34 @@
 namespace py = pybind11;
 
 template <typename DataType>
-py::array_t<DataType> pybind_apply_calibration(
-    py::array_t<uint16_t, py::array::c_style | py::array::forcecast> data,
-    py::array_t<DataType, py::array::c_style | py::array::forcecast> pedestal,
-    py::array_t<DataType, py::array::c_style | py::array::forcecast>
-        calibration,
-    int n_threads = 4) {
+py::array_t<DataType>
+pybind_apply_calibration(py::array_t<uint16_t> data,
+                         py::array_t<DataType> pedestal,
+                         py::array_t<DataType> calibration, int n_threads = 4) {
 
     auto data_span = make_view_3d(data); // data is always 3D
     /* No pointer is passed, so NumPy will allocate the buffer */
     auto result = py::array_t<DataType>(data_span.shape());
     auto res = make_view_3d(result);
-    if (data.ndim() == 3 && pedestal.ndim() == 3 && calibration.ndim() == 3) {
+    if (pedestal.ndim() == 3 && calibration.ndim() == 3) {
         auto ped = make_view_3d(pedestal);
         auto cal = make_view_3d(calibration);
         aare::apply_calibration<DataType, 3>(res, data_span, ped, cal,
                                              n_threads);
-    } else if (data.ndim() == 3 && pedestal.ndim() == 2 &&
-               calibration.ndim() == 2) {
+    } else if (pedestal.ndim() == 2 && calibration.ndim() == 2) {
         auto ped = make_view_2d(pedestal);
         auto cal = make_view_2d(calibration);
         aare::apply_calibration<DataType, 2>(res, data_span, ped, cal,
                                              n_threads);
     } else {
         throw std::runtime_error(
-            "Invalid number of dimensions for data, pedestal or calibration");
+            "Invalid number of dimensions for pedestal or calibration");
     }
     return result;
 }
 
-py::array_t<int> pybind_count_switching_pixels(
-    py::array_t<uint16_t, py::array::c_style | py::array::forcecast> data,
-    ssize_t n_threads = 4) {
+py::array_t<int> pybind_count_switching_pixels(py::array_t<uint16_t> data,
+                                               ssize_t n_threads = 4) {
 
     auto data_span = make_view_3d(data);
     auto arr = new NDArray<int, 2>{};
@@ -51,9 +47,8 @@ py::array_t<int> pybind_count_switching_pixels(
 }
 
 template <typename T>
-py::array_t<T> pybind_calculate_pedestal(
-    py::array_t<uint16_t, py::array::c_style | py::array::forcecast> data,
-    ssize_t n_threads) {
+py::array_t<T> pybind_calculate_pedestal(py::array_t<uint16_t> data,
+                                         ssize_t n_threads) {
 
     auto data_span = make_view_3d(data);
     auto arr = new NDArray<T, 3>{};
@@ -62,9 +57,8 @@ py::array_t<T> pybind_calculate_pedestal(
 }
 
 template <typename T>
-py::array_t<T> pybind_calculate_pedestal_g0(
-    py::array_t<uint16_t, py::array::c_style | py::array::forcecast> data,
-    ssize_t n_threads) {
+py::array_t<T> pybind_calculate_pedestal_g0(py::array_t<uint16_t> data,
+                                            ssize_t n_threads) {
 
     auto data_span = make_view_3d(data);
     auto arr = new NDArray<T, 2>{};
@@ -73,10 +67,27 @@ py::array_t<T> pybind_calculate_pedestal_g0(
 }
 
 void bind_calibration(py::module &m) {
+    py::options options;
+    options.disable_function_signatures();
+
     m.def("apply_calibration", &pybind_apply_calibration<double>,
           py::arg("raw_data").noconvert(), py::kw_only(),
           py::arg("pd").noconvert(), py::arg("cal").noconvert(),
-          py::arg("n_threads") = 4);
+          py::arg("n_threads") = 4,
+          R"(
+        Apply the calibration to the raw data and return the result as a 3D array of doubles.
+        
+        Parameters
+        ----------
+        raw_data : array_like
+            3D array of shape ``(frames, rows, cols)``.
+        pd : array_like
+            Pedestal array.
+        cal : array_like
+            Calibration array.
+        n_threads : int, optional
+            Number of threads to use. Defaults to 4.
+        )");
 
     m.def("apply_calibration", &pybind_apply_calibration<float>,
           py::arg("raw_data").noconvert(), py::kw_only(),

python/src/bind_testing.hpp

@@ -0,0 +1,65 @@
+// SPDX-License-Identifier: MPL-2.0
+// Probes used by the python test-suite to exercise the helpers in
+// np_helper.hpp. Internal, no API guarantees.
+#pragma once
+
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include "np_helper.hpp"
+
+#include <algorithm>
+#include <cstdint>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <type_traits>
+
+namespace py = pybind11;
+using namespace ::aare;
+
+// Arrays are taken with the default flags (no c_style) so that nothing is
+// made contiguous by pybind11 before it reaches make_view
+template <typename T, ssize_t Ndim>
+void define_make_view_probes(py::module &m, const std::string &suffix) {
+
+    // Copy what C++ sees through the view back to python
+    m.def(("view_roundtrip_" + suffix).c_str(),
+          [](py::array_t<T> arr) {
+              auto view = make_view<Ndim>(arr);
+              return return_image_data(new NDArray<T, Ndim>(view));
+          },
+          py::arg("arr").noconvert());
+
+    // Same as above but through a read only view
+    m.def(("const_view_roundtrip_" + suffix).c_str(),
+          [](py::array_t<T> arr) {
+              auto view = make_const_view<Ndim>(arr);
+              static_assert(
+                  std::is_same_v<decltype(view), NDView<const T, Ndim>>);
+              auto out = new NDArray<T, Ndim>(view.shape());
+              std::copy(view.begin(), view.end(), out->begin());
+              return return_image_data(out);
+          },
+          py::arg("arr").noconvert());
+
+    // Write through the view
+    m.def(("view_fill_" + suffix).c_str(),
+          [](py::array_t<T> arr, T value) {
+              auto view = make_view<Ndim>(arr);
+              view = value;
+          },
+          py::arg("arr").noconvert(), py::arg("value"));
+}
+
+void define_testing_bindings(py::module &m) {
+    auto testing = m.def_submodule(
+        "testing", "Internal helpers for the test-suite, no API guarantees");
+
+    // noconvert + overloads gives dispatch on dtype
+    define_make_view_probes<double, 1>(testing, "1d");
+    define_make_view_probes<double, 2>(testing, "2d");
+    define_make_view_probes<double, 3>(testing, "3d");
+    define_make_view_probes<uint16_t, 1>(testing, "1d");
+    define_make_view_probes<uint16_t, 2>(testing, "2d");
+    define_make_view_probes<uint16_t, 3>(testing, "3d");
+}

python/src/ctb_raw_file.hpp

@@ -24,177 +24,119 @@
 namespace py = pybind11;
 using namespace ::aare;
 
+// One binding for the three ADC SAR decoders. Each is registered under
+// its original Python name so existing transforms keep working.
+template <void (*decode)(NDView<const uint8_t, 2>, NDView<uint16_t, 2>)>
+py::array_t<uint16_t> decode64to16(py::array_t<uint8_t> input) {
+    auto bytes = make_const_view_2d(input);
+    constexpr ssize_t bytes_per_word = sizeof(uint64_t);
+    auto output = py::array_t<uint16_t>(
+        {bytes.shape(0), bytes.shape(1) / bytes_per_word});
+    decode(bytes, make_view_2d(output));
+    return output;
+}
+
+constexpr auto decode64to16_doc =
+    "Decode packed 64-bit ADC samples. Takes a two-dimensional, "
+    "C-contiguous uint8 array whose row length is a multiple of 8 bytes "
+    "and returns a uint16 array with one value per 64-bit word.";
+
 void define_ctb_raw_file_io_bindings(py::module &m) {
 
-    m.def("adc_sar_05_06_07_08decode64to16", [](py::array_t<uint8_t> input) {
-        if (input.ndim() != 2) {
-            throw std::runtime_error(
-                "Only 2D arrays are supported at this moment");
-        }
+    m.def("adc_sar_05_06_07_08decode64to16",
+          &decode64to16<adc_sar_05_06_07_08decode64to16>,
+          py::arg("input").noconvert(), decode64to16_doc);
 
-        // Create a 2D output array with the same shape as the input
-        std::vector<ssize_t> shape{input.shape(0),
-                                   input.shape(1) /
-                                       static_cast<ssize_t>(bits_per_byte)};
-        py::array_t<uint16_t> output(shape);
+    m.def("adc_sar_05_decode64to16", &decode64to16<adc_sar_05_decode64to16>,
+          py::arg("input").noconvert(), decode64to16_doc);
 
-        // Create a view of the input and output arrays
-        NDView<uint64_t, 2> input_view(
-            reinterpret_cast<uint64_t *>(input.mutable_data()),
-            {output.shape(0), output.shape(1)});
-        NDView<uint16_t, 2> output_view(output.mutable_data(),
-                                        {output.shape(0), output.shape(1)});
+    m.def("adc_sar_04_decode64to16", &decode64to16<adc_sar_04_decode64to16>,
+          py::arg("input").noconvert(), decode64to16_doc);
 
-        adc_sar_05_06_07_08decode64to16(input_view, output_view);
+    m.def(
+        "apply_custom_weights",
+        [](py::array_t<uint16_t> input, py::array_t<double> weights) {
+            auto input_view = make_view_1d(input);
+            auto weights_view = make_view_1d(weights);
+            auto output = py::array_t<double>(input_view.size());
+            apply_custom_weights(input_view, make_view_1d(output),
+                                 weights_view);
+            return output;
+        },
+        py::arg("input").noconvert(), py::arg("weights").noconvert(),
+        "Apply per-bit weights to every value of a one-dimensional, "
+        "C-contiguous uint16 array and return a float64 array.");
 
-        return output;
-    });
+    m.def(
+        "expand24to32bit",
+        [](py::array_t<uint8_t> input, uint32_t offset) {
+            constexpr ssize_t bytes_per_channel = 3; // 24 bit
+            auto input_view = make_const_view_1d(input);
+            auto output =
+                py::array_t<uint32_t>(input_view.size() / bytes_per_channel);
+            aare::expand24to32bit(input_view, make_view_1d(output),
+                                  aare::BitOffset(offset));
+            return output;
+        },
+        py::arg("input").noconvert(), py::arg("offset"),
+        "Expand packed 24-bit values from a one-dimensional, C-contiguous "
+        "uint8 array into a uint32 array.");
 
-    m.def("adc_sar_05_decode64to16", [](py::array_t<uint8_t> input) {
-        if (input.ndim() != 2) {
-            throw std::runtime_error(
-                "Only 2D arrays are supported at this moment");
-        }
+    m.def(
+        "expand4to8bit",
+        [](py::array_t<uint8_t> input) {
+            auto input_view = make_const_view_1d(input);
+            auto output = py::array_t<uint8_t>(input_view.size() * 2);
+            aare::expand4to8bit(input_view, make_view_1d(output));
+            return output;
+        },
+        py::arg("input").noconvert(),
+        "Expand two 4-bit values per byte of a one-dimensional, C-contiguous "
+        "uint8 array into a uint8 array of twice the length.");
 
-        // Create a 2D output array with the same shape as the input
-        std::vector<ssize_t> shape{input.shape(0),
-                                   input.shape(1) /
-                                       static_cast<ssize_t>(bits_per_byte)};
-        py::array_t<uint16_t> output(shape);
+    m.def(
+        "decode_my302",
+        [](py::array_t<uint8_t> input, uint32_t offset) {
+            // Physical layout of the chip
+            constexpr ssize_t channels = 64;
+            constexpr ssize_t counters = 3;
+            constexpr ssize_t bytes_per_channel = 3; // 24 bit
+            constexpr ssize_t n_outputs = 2;
 
-        // Create a view of the input and output arrays
-        NDView<uint64_t, 2> input_view(
-            reinterpret_cast<uint64_t *>(input.mutable_data()),
-            {output.shape(0), output.shape(1)});
-        NDView<uint16_t, 2> output_view(output.mutable_data(),
-                                        {output.shape(0), output.shape(1)});
+            auto input_view = make_const_view_1d(input);
 
-        adc_sar_05_decode64to16(input_view, output_view);
+            ssize_t expected_size = channels * counters * bytes_per_channel;
 
-        return output;
-    });
+            // If we have an offset we need one extra byte per output
+            aare::BitOffset bitoff(offset);
+            if (bitoff.value())
+                expected_size += n_outputs;
 
-    m.def("adc_sar_04_decode64to16", [](py::array_t<uint8_t> input) {
-        if (input.ndim() != 2) {
-            throw std::runtime_error(
-                "Only 2D arrays are supported at this moment");
-        }
+            if (input_view.size() != expected_size) {
+                throw py::value_error(
+                    fmt::format("{} Expected an input size of {} bytes. Called "
+                                "with input size of {}",
+                                LOCATION, expected_size, input_view.size()));
+            }
 
-        // Create a 2D output array with the same shape as the input
-        std::vector<ssize_t> shape{input.shape(0),
-                                   input.shape(1) /
-                                       static_cast<ssize_t>(bits_per_byte)};
-        py::array_t<uint16_t> output(shape);
+            auto output = py::array_t<uint32_t>(channels * counters);
+            auto output_view = make_view_1d(output);
 
-        // Create a view of the input and output arrays
-        NDView<uint64_t, 2> input_view(
-            reinterpret_cast<uint64_t *>(input.mutable_data()),
-            {output.shape(0), output.shape(1)});
-        NDView<uint16_t, 2> output_view(output.mutable_data(),
-                                        {output.shape(0), output.shape(1)});
+            const ssize_t in_step = input_view.size() / n_outputs;
+            const ssize_t out_step = output_view.size() / n_outputs;
+            for (ssize_t i = 0; i != n_outputs; ++i) {
+                NDView<const uint8_t, 1> in_part(
+                    input_view.data() + in_step * i, {in_step});
+                NDView<uint32_t, 1> out_part(output_view.data() + out_step * i,
+                                             {out_step});
+                aare::expand24to32bit(in_part, out_part, bitoff);
+            }
 
-        adc_sar_04_decode64to16(input_view, output_view);
-
-        return output;
-    });
-
-    m.def("apply_custom_weights",
-          [](py::array_t<uint16_t, py::array::c_style | py::array::forcecast>
-                 &input,
-             py::array_t<double, py::array::c_style | py::array::forcecast>
-                 &weights) {
-              // Create new array with same shape as the input array
-              // (uninitialized values)
-              py::buffer_info buf = input.request();
-              py::array_t<double> output(buf.shape);
-
-              // Use NDViews to call into the C++ library
-              auto weights_view = make_view_1d(weights);
-              NDView<uint16_t, 1> input_view(input.mutable_data(),
-                                             {input.size()});
-              NDView<double, 1> output_view(output.mutable_data(),
-                                            {output.size()});
-
-              apply_custom_weights(input_view, output_view, weights_view);
-              return output;
-          });
-
-    m.def("expand24to32bit",
-          [](py::array_t<uint8_t, py::array::c_style | py::array::forcecast>
-                 &input,
-             uint32_t offset) {
-              aare::BitOffset bitoff(offset);
-              py::buffer_info buf = input.request();
-
-              constexpr uint32_t bytes_per_channel = 3; // 24 bit
-              py::array_t<uint32_t> output(buf.size / bytes_per_channel);
-
-              NDView<uint8_t, 1> input_view(input.mutable_data(),
-                                            {input.size()});
-              NDView<uint32_t, 1> output_view(output.mutable_data(),
-                                              {output.size()});
-
-              aare::expand24to32bit(input_view, output_view, bitoff);
-              return output;
-          });
-
-    m.def("expand4to8bit",
-          [](py::array_t<uint8_t, py::array::c_style | py::array::forcecast>
-                 &input) {
-              py::buffer_info buf = input.request();
-
-              py::array_t<uint8_t> output(buf.size * 2);
-
-              NDView<uint8_t, 1> input_view(input.mutable_data(),
-                                            {input.size()});
-              NDView<uint8_t, 1> output_view(output.mutable_data(),
-                                             {output.size()});
-
-              aare::expand4to8bit(input_view, output_view);
-              return output;
-          });
-
-    m.def("decode_my302",
-          [](py::array_t<uint8_t, py::array::c_style | py::array::forcecast>
-                 &input,
-             uint32_t offset) {
-              // Physical layout of the chip
-              constexpr size_t channels = 64;
-              constexpr size_t counters = 3;
-              constexpr size_t bytes_per_channel = 3; // 24 bit
-              constexpr int n_outputs = 2;
-
-              ssize_t expected_size = channels * counters * bytes_per_channel;
-
-              // If whe have an offset we need one extra byte per output
-              aare::BitOffset bitoff(offset);
-              if (bitoff.value())
-                  expected_size += n_outputs;
-
-              if (input.size() != expected_size) {
-                  throw std::runtime_error(fmt::format(
-                      "{} Expected an input size of {} bytes. Called "
-                      "with input size of {}",
-                      LOCATION, expected_size, input.size()));
-              }
-
-              py::buffer_info buf = input.request();
-              py::array_t<uint32_t> output(channels * counters);
-
-              for (int i = 0; i != n_outputs; ++i) {
-                  auto step = input.size() / n_outputs;
-                  auto out_step = output.size() / n_outputs;
-                  NDView<uint8_t, 1> input_view(input.mutable_data() + step * i,
-                                                {input.size() / n_outputs});
-                  NDView<uint32_t, 1> output_view(output.mutable_data() +
-                                                      out_step * i,
-                                                  {output.size() / n_outputs});
-
-                  aare::expand24to32bit(input_view, output_view, bitoff);
-              }
-
-              return output;
-          });
+            return output;
+        },
+        py::arg("input").noconvert(), py::arg("offset"),
+        "Decode a Mythen 302 readout from a one-dimensional, C-contiguous "
+        "uint8 array into a uint32 array of 64 channels times 3 counters.");
 
     py::class_<CtbRawFile>(m, "CtbRawFile")
         .def(py::init<const std::filesystem::path &>())

python/src/file.hpp

@@ -2,6 +2,7 @@
 #include "aare/CtbRawFile.hpp"
 #include "aare/File.hpp"
 #include "aare/Frame.hpp"
+#include "aare/ROI.hpp"
 #include "aare/RawFile.hpp"
 #include "aare/RawMasterFile.hpp"
 #include "aare/RawSubFile.hpp"
@@ -27,16 +28,6 @@ using namespace ::aare;
 
 void define_file_io_bindings(py::module &m) {
 
-    py::enum_<DetectorType>(m, "DetectorType")
-        .value("Jungfrau", DetectorType::Jungfrau)
-        .value("Eiger", DetectorType::Eiger)
-        .value("Mythen3", DetectorType::Mythen3)
-        .value("Moench", DetectorType::Moench)
-        .value("Moench03", DetectorType::Moench03)
-        .value("Moench03_old", DetectorType::Moench03_old)
-        .value("ChipTestBoard", DetectorType::ChipTestBoard)
-        .value("Unknown", DetectorType::Unknown);
-
     PYBIND11_NUMPY_DTYPE(DetectorHeader, frameNumber, expLength, packetNumber,
                          bunchId, timestamp, modId, row, column, reserved,
                          debug, roundRNumber, detType, version, packetMask);
@@ -170,28 +161,5 @@ void define_file_io_bindings(py::module &m) {
         .def_property_readonly("stop", &ScanParameters::stop)
         .def_property_readonly("step", &ScanParameters::step);
 
-    py::class_<ROI>(m, "ROI")
-        .def(py::init<>())
-        .def(py::init<ssize_t, ssize_t, ssize_t, ssize_t>(), py::arg("xmin"),
-             py::arg("xmax"), py::arg("ymin"), py::arg("ymax"))
-        .def_readwrite("xmin", &ROI::xmin)
-        .def_readwrite("xmax", &ROI::xmax)
-        .def_readwrite("ymin", &ROI::ymin)
-        .def_readwrite("ymax", &ROI::ymax)
-        .def("__str__",
-             [](const ROI &self) {
-                 return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}",
-                                    self.xmin, self.xmax, self.ymin, self.ymax);
-             })
-        .def("__repr__",
-             [](const ROI &self) {
-                 return fmt::format(
-                     "<ROI: xmin: {} xmax: {} ymin: {} ymax: {}>", self.xmin,
-                     self.xmax, self.ymin, self.ymax);
-             })
-        .def("__iter__", [](const ROI &self) {
-            return py::make_iterator(&self.xmin, &self.ymax + 1); // NOLINT
-        });
-
 #pragma GCC diagnostic pop
 }

python/src/fit.hpp

@@ -21,10 +21,43 @@ fit_dispatch(const aare::FitModel<Model> &model,
 
 template <typename Model> void bind_fit_model(py::module &m, const char *name) {
     using FM = aare::FitModel<Model>;
+
     auto cls = py::class_<FM>(m, name);
-    cls.def(py::init<unsigned int, unsigned int, double, bool>(),
+
+    cls.def(py::init<unsigned int, unsigned int, double, bool,
+                     aare::Minimizer>(),
             py::arg("strategy") = 0, py::arg("max_calls") = 100,
-            py::arg("tolerance") = 0.5, py::arg("compute_errors") = false)
+            py::arg("tolerance") = 0.5, py::arg("compute_errors") = false,
+            py::arg("minimizer") = aare::Minimizer::Migrad,
+            R"doc(
+            Fit configuration for this model.
+
+            Parameters
+            ----------
+            strategy : int
+                Minuit2 strategy (0 = fast, 1 = Minuit2's default). Ignored
+                by ``Minimizer.LevenbergMarquardt`` and
+                ``Minimizer.VarPro``.
+            max_calls : int
+                Work budget per pixel: Minuit2 function calls, or model
+                evaluations for ``Minimizer.LevenbergMarquardt`` and
+                ``Minimizer.VarPro`` (one per iteration, two
+                when a step crosses a limit). 0 selects Minuit's default
+                budget.
+            tolerance : float
+                EDM tolerance in Minuit's convention. Migrad,
+                LevenbergMarquardt and VarPro stop below
+                ``0.002 * tolerance``, Fumili below ``1e-4 * tolerance``.
+            compute_errors : bool
+                Also return parameter errors: from Hesse with Migrad, from
+                the linearised covariance with Fumili, LevenbergMarquardt
+                and VarPro. Fixed parameters and parameters
+                ending on a limit report 0.
+            minimizer : Minimizer
+                ``Minimizer.Migrad`` (default), ``Minimizer.Fumili``,
+                ``Minimizer.LevenbergMarquardt`` or
+                ``Minimizer.VarPro``.
+            )doc")
         .def("SetParLimits",
              py::overload_cast<unsigned int, double, double>(&FM::SetParLimits),
              py::arg("idx"), py::arg("lo"), py::arg("hi"))
@@ -59,8 +92,10 @@ template <typename Model> void bind_fit_model(py::module &m, const char *name) {
         .def_property("tolerance", &FM::tolerance, &FM::SetTolerance)
         .def_property("compute_errors", &FM::compute_errors,
                       &FM::SetComputeErrors)
+        .def_property("minimizer", &FM::minimizer, &FM::SetMinimizer)
         .def(
-            "__call__",
+            "__call__", // conversion ok, we want to be able to call with any
+                        // dtype
             [](const FM & /*self*/,
                py::array_t<double, py::array::c_style | py::array::forcecast> x,
                py::array_t<double, py::array::c_style | py::array::forcecast>
@@ -78,7 +113,7 @@ template <typename Model> void bind_fit_model(py::module &m, const char *name) {
             },
             py::arg("x"), py::arg("par"))
         .def(
-            "fit",
+            "fit", // conversion ok
             [](const FM &self,
                py::array_t<double, py::array::c_style | py::array::forcecast> x,
                py::array_t<double, py::array::c_style | py::array::forcecast> y,
@@ -88,6 +123,11 @@ template <typename Model> void bind_fit_model(py::module &m, const char *name) {
             R"doc(
             Fit this model to 1D or 3D data using Minuit2.
 
+            The minimizer is selected by the ``minimizer`` property
+            (``Minimizer.Migrad`` by default, ``Minimizer.Fumili``,
+            ``Minimizer.LevenbergMarquardt`` or
+            ``Minimizer.VarPro``).
+
             Parameters
             ----------
             x : array_like, shape (n_scan,)
@@ -202,11 +242,6 @@ fit_dispatch(const aare::FitModel<Model> &model,
                 py::array_t<double, py::array::c_style | py::array::forcecast>>(
                 y_err_obj);
 
-            if (y_err.ndim() != 3) {
-                throw std::runtime_error(
-                    "For 3D input y, y_err must also be 3D.");
-            }
-
             auto err_out =
                 new NDArray<double, 3>({y.shape(0), y.shape(1), npar}, 0.0);
             auto y_view_err = make_view_3d(y_err);
@@ -228,8 +263,23 @@ fit_dispatch(const aare::FitModel<Model> &model,
                                 "chi2"_a = return_image_data(chi2_out));
             }
         } else {
-
+            // Unweighted fit; parameter errors are still available when
+            // requested, as in the 1D path.
             NDView<double, 3> dummy_err{};
+
+            if (model.compute_errors()) {
+                auto err_out =
+                    new NDArray<double, 3>({y.shape(0), y.shape(1), npar}, 0.0);
+
+                aare::fit_3d<Model>(model, x_view, y_view, dummy_err,
+                                    par_out->view(), err_out->view(),
+                                    chi2_out->view(), n_threads);
+
+                return py::dict("par"_a = return_image_data(par_out),
+                                "par_err"_a = return_image_data(err_out),
+                                "chi2"_a = return_image_data(chi2_out));
+            }
+
             NDView<double, 3> dummy_err_out{};
 
             {
@@ -253,11 +303,6 @@ fit_dispatch(const aare::FitModel<Model> &model,
                 py::array_t<double, py::array::c_style | py::array::forcecast>>(
                 y_err_obj);
 
-            if (y_err.ndim() != 1) {
-                throw std::runtime_error(
-                    "For 1D input y, y_err must also be 1D.");
-            }
-
             auto y_view_err = make_view_1d(y_err);
             {
                 py::gil_scoped_release release; // release GIL for parallel loop
@@ -278,7 +323,82 @@ fit_dispatch(const aare::FitModel<Model> &model,
     }
 }
 
+// Resolve the model type behind a Python object and run the fit.
+py::object dispatch_any_model(
+    py::object model_obj,
+    py::array_t<double, py::array::c_style | py::array::forcecast> x,
+    py::array_t<double, py::array::c_style | py::array::forcecast> y,
+    py::object y_err_obj, int n_threads) {
+    using namespace aare::model;
+
+#define AARE_DISPATCH_MODEL(Model)                                             \
+    if (py::isinstance<aare::FitModel<Model>>(model_obj)) {                    \
+        const auto &mdl = model_obj.cast<const aare::FitModel<Model> &>();     \
+        return fit_dispatch<Model>(mdl, x, y, y_err_obj, n_threads);           \
+    }
+
+    AARE_DISPATCH_MODEL(Pol1)
+    AARE_DISPATCH_MODEL(Pol2)
+    AARE_DISPATCH_MODEL(Gaussian)
+    AARE_DISPATCH_MODEL(GaussianErfcPlateau)
+    AARE_DISPATCH_MODEL(GaussianChargeSharing)
+    AARE_DISPATCH_MODEL(GaussianChargeSharingKb)
+    AARE_DISPATCH_MODEL(RisingScurve)
+    AARE_DISPATCH_MODEL(FallingScurve)
+
+#undef AARE_DISPATCH_MODEL
+
+    throw std::runtime_error(
+        "Unknown model type. Expected Pol1, Pol2, Gaussian, "
+        "GaussianErfcPlateau, GaussianChargeSharing, GaussianChargeSharingKb, "
+        "RisingScurve or FallingScurve.");
+}
+
 void define_fit_bindings(py::module &m) {
+    // The enum must exist before it is used as a constructor default below.
+    py::enum_<aare::Minimizer>(m, "Minimizer", R"doc(
+        Minimizer used by the fit models.
+
+        ``Migrad`` (default) is Minuit2's variable-metric minimizer. With
+        ``compute_errors`` its parameter errors come from Hesse.
+
+        ``Fumili`` is Minuit2's Gauss-Newton minimizer with a trust region.
+        It takes the gradient and a linearised Hessian from the analytic
+        derivatives of the model and usually needs far fewer function
+        evaluations than Migrad. Its parameter errors come from the
+        covariance of the linearised Hessian, without an extra Hesse step.
+        Minuit2's Fumili cannot converge once a two-sided limit becomes
+        active; a pixel without a valid Fumili minimum is refitted with
+        Migrad.
+
+        ``LevenbergMarquardt`` is the built-in damped Gauss-Newton solver. It
+        uses the analytic Jacobian, reflects steps at parameter limits and
+        reuses its buffers between pixels, so data cubes are fitted without
+        per-pixel allocations. Its parameter errors are
+        ``sqrt(diag((J^T J)^-1))``; ``max_calls`` counts model evaluations
+        and ``strategy`` is ignored.
+
+        ``VarPro`` is the built-in variable projection solver.
+        Every bundled model is linear in some of its parameters (the
+        amplitude of a Gaussian, the baseline and amplitude of an S-curve,
+        every parameter of a polynomial); each trial point solves those
+        exactly and the Levenberg-Marquardt iteration runs over the
+        remaining nonlinear parameters only. It needs fewer evaluations than
+        ``LevenbergMarquardt`` and cannot be misled by poor start values of
+        the linear parameters, which it ignores. A pixel whose linear
+        solution violates a limit on a linear parameter, or that does not
+        converge, is refitted with ``LevenbergMarquardt``. Its parameter
+        errors are the same Gauss-Newton estimates; ``max_calls`` counts
+        evaluations and ``strategy`` is ignored.
+
+        With every minimizer, fixed parameters and parameters that end on a
+        limit report an error of 0 and a failed fit returns zeros.
+        )doc")
+        .value("Migrad", aare::Minimizer::Migrad)
+        .value("Fumili", aare::Minimizer::Fumili)
+        .value("LevenbergMarquardt", aare::Minimizer::LevenbergMarquardt)
+        .value("VarPro", aare::Minimizer::VarPro);
+
     // ── Bind model classes ──────────────────────────────────────────
     bind_fit_model<aare::model::Gaussian>(m, "Gaussian");
     bind_fit_model<aare::model::GaussianErfcPlateau>(m, "GaussianErfcPlateau");
@@ -293,91 +413,24 @@ void define_fit_bindings(py::module &m) {
 
     m.def(
         "fit",
-        [](py::object model_obj,
+        [](py::object model_obj, // conversion ok, we want to be able to call
+                                 // with any dtype
            py::array_t<double, py::array::c_style | py::array::forcecast> x,
            py::array_t<double, py::array::c_style | py::array::forcecast> y,
            py::object y_err_obj, int n_threads) -> py::object {
-            using namespace aare::model;
-
-            // ── Polynomial of degree 1 ───────
-            if (py::isinstance<aare::FitModel<Pol1>>(model_obj)) {
-
-                const auto &mdl =
-                    model_obj.cast<const aare::FitModel<Pol1> &>();
-                return fit_dispatch<Pol1>(mdl, x, y, y_err_obj, n_threads);
-            }
-
-            // ── Polynomial of degree 2 ───────
-            if (py::isinstance<aare::FitModel<Pol2>>(model_obj)) {
-                const auto &mdl =
-                    model_obj.cast<const aare::FitModel<Pol2> &>();
-                return fit_dispatch<Pol2>(mdl, x, y, y_err_obj, n_threads);
-            }
-
-            // ── Gaussian ───────
-            if (py::isinstance<aare::FitModel<Gaussian>>(model_obj)) {
-                const auto &mdl =
-                    model_obj.cast<const aare::FitModel<Gaussian> &>();
-                return fit_dispatch<Gaussian>(mdl, x, y, y_err_obj, n_threads);
-            }
-
-            // ── GaussianErfcPlateau ───────
-            if (py::isinstance<aare::FitModel<GaussianErfcPlateau>>(
-                    model_obj)) {
-                const auto &mdl =
-                    model_obj
-                        .cast<const aare::FitModel<GaussianErfcPlateau> &>();
-                return fit_dispatch<GaussianErfcPlateau>(mdl, x, y, y_err_obj,
-                                                         n_threads);
-            }
-
-            // ── GaussianChargeSharing ───────
-            if (py::isinstance<aare::FitModel<GaussianChargeSharing>>(
-                    model_obj)) {
-                const auto &mdl =
-                    model_obj
-                        .cast<const aare::FitModel<GaussianChargeSharing> &>();
-                return fit_dispatch<GaussianChargeSharing>(mdl, x, y, y_err_obj,
-                                                           n_threads);
-            }
-
-            // ── GaussianChargeSharingKb ───────
-            if (py::isinstance<aare::FitModel<GaussianChargeSharingKb>>(
-                    model_obj)) {
-                const auto &mdl = model_obj.cast<
-                    const aare::FitModel<GaussianChargeSharingKb> &>();
-                return fit_dispatch<GaussianChargeSharingKb>(
-                    mdl, x, y, y_err_obj, n_threads);
-            }
-
-            // ── Rising Scurve ───────
-            if (py::isinstance<aare::FitModel<RisingScurve>>(model_obj)) {
-                const auto &mdl =
-                    model_obj.cast<const aare::FitModel<RisingScurve> &>();
-                return fit_dispatch<RisingScurve>(mdl, x, y, y_err_obj,
-                                                  n_threads);
-            }
-
-            // ── Falling Scurve ───────
-            if (py::isinstance<aare::FitModel<FallingScurve>>(model_obj)) {
-                const auto &mdl =
-                    model_obj.cast<const aare::FitModel<FallingScurve> &>();
-                return fit_dispatch<FallingScurve>(mdl, x, y, y_err_obj,
-                                                   n_threads);
-            }
-
-            throw std::runtime_error(
-                "Unknown model type. Expected Pol1, Pol2, Gaussian, "
-                "RisingScurve or FallingScurve.");
+            return dispatch_any_model(model_obj, x, y, y_err_obj, n_threads);
         },
         R"(
-        Fit a model to 1D or 3D data using Minuit2.
+        Fit a model to 1D or 3D data with the minimizer selected by the
+        model's ``minimizer`` property.
  
         Parameters
         ----------
-        model : Pol1, Pol2, Gaussian, RisingScurve, or FallingScurve
-            Configured model object.  User-set limits, fixed parameters,
-            and start values take precedence over automatic estimates.
+        model : object
+            Configured model object: Pol1, Pol2, Gaussian, GaussianErfcPlateau,
+            GaussianChargeSharing, GaussianChargeSharingKb, RisingScurve or
+            FallingScurve.  User-set limits, fixed parameters, and start values
+            take precedence over automatic estimates.
         x : array_like, shape (n_scan,)
             Scan points (e.g. energy or threshold values).
         y : array_like, shape (n_scan,) or (rows, cols, n_scan)
@@ -389,15 +442,12 @@ void define_fit_bindings(py::module &m) {
  
         Returns
         -------
-        For 1D input:
-            numpy array of shape (2*npar+1,) if compute_errors else (npar+1,).
-            Layout: [params..., (errors...,) chi2].
- 
-        For 3D input:
-            dict with keys:
-              "par"     : (rows, cols, npar) fitted parameters.
-              "par_err" : (rows, cols, npar) parameter errors (if compute_errors).
-              "chi2"    : (rows, cols)       chi-squared per pixel.
+        dict
+            ``par`` holds the fitted parameters, shape (npar,) for 1D input or
+            (rows, cols, npar) for 3D input.  ``chi2`` holds the chi-squared,
+            shape (1,) or (rows, cols).  ``par_err`` holds the parameter
+            errors when ``compute_errors`` is enabled; fixed parameters and
+            parameters on a limit report 0.  A failed fit returns zeros.
         )",
         py::arg("model"), py::arg("x"), py::arg("y"),
         py::arg("y_err") = py::none(), py::arg("n_threads") = 4);

python/src/module.cpp

@@ -1,6 +1,8 @@
 // SPDX-License-Identifier: MPL-2.0
 // Files with bindings to the different classes
 
+#include "module_config.hpp"
+
 // New style file naming
 #include "bind_Cluster.hpp"
 #include "bind_ClusterCollector.hpp"
@@ -11,25 +13,28 @@
 #include "bind_ClusterVector.hpp"
 #include "bind_Defs.hpp"
 #include "bind_Eta.hpp"
+#include "bind_FastPedestal.hpp"
 #include "bind_Interpolator.hpp"
 #include "bind_MultiThreadedFileReader.hpp"
+#include "bind_Pedestal.hpp"
 #include "bind_PedestalTrackingPixelHistogram.hpp"
 #include "bind_PixelHistogram.hpp"
 #include "bind_PixelMap.hpp"
 #include "bind_RawFile.hpp"
 #include "bind_calibration.hpp"
+#include "bind_testing.hpp"
 
 // TODO! migrate the other names
 #include "ctb_raw_file.hpp"
 #include "file.hpp"
 #include "fit.hpp"
 #include "jungfrau_data_file.hpp"
-#include "pedestal.hpp"
 #include "raw_master_file.hpp"
 #include "raw_sub_file.hpp"
 #include "var_cluster.hpp"
 
 // Pybind stuff
+#include <cstdint>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 
@@ -75,11 +80,16 @@ PYBIND11_MODULE(_aare, m) {
     define_pedestal_tracking_pixel_histogram_bindings(m);
     define_pedestal_bindings<double>(m, "Pedestal_d");
     define_pedestal_bindings<float>(m, "Pedestal_f");
+    define_pedestal_bindings<int16_t>(m, "Pedestal_i16");
+    define_fast_pedestal_bindings<double>(m, "FastPedestal_d");
+    define_fast_pedestal_bindings<float>(m, "FastPedestal_f");
+    define_fast_pedestal_bindings<int16_t>(m, "FastPedestal_i16");
     define_fit_bindings(m);
     define_interpolation_bindings(m);
     define_jungfrau_data_file_io_bindings(m);
 
     bind_calibration(m);
+    define_testing_bindings(m);
 
     DEFINE_CLUSTER_BINDINGS(int, 3, 3, uint16_t, i);
     DEFINE_CLUSTER_BINDINGS(double, 3, 3, uint16_t, d);
@@ -106,6 +116,7 @@ PYBIND11_MODULE(_aare, m) {
     DEFINE_BINDINGS_CLUSTERFINDER(int, 3, 3, uint16_t, i);
     DEFINE_BINDINGS_CLUSTERFINDER(double, 3, 3, uint16_t, d);
     DEFINE_BINDINGS_CLUSTERFINDER(float, 3, 3, uint16_t, f);
+    DEFINE_BINDINGS_CLUSTERFINDER(int16_t, 3, 3, uint16_t, i16);
 
     DEFINE_BINDINGS_CLUSTERFINDER(int, 5, 5, uint16_t, i);
     DEFINE_BINDINGS_CLUSTERFINDER(double, 5, 5, uint16_t, d);
@@ -119,9 +130,6 @@ PYBIND11_MODULE(_aare, m) {
     DEFINE_BINDINGS_CLUSTERFINDER(double, 9, 9, uint16_t, d);
     DEFINE_BINDINGS_CLUSTERFINDER(float, 9, 9, uint16_t, f);
 
-    define_3x3_reduction<int, 3, 3, uint16_t>(m);
-    define_3x3_reduction<double, 3, 3, uint16_t>(m);
-    define_3x3_reduction<float, 3, 3, uint16_t>(m);
     define_3x3_reduction<int, 5, 5, uint16_t>(m);
     define_3x3_reduction<double, 5, 5, uint16_t>(m);
     define_3x3_reduction<float, 5, 5, uint16_t>(m);
@@ -132,9 +140,6 @@ PYBIND11_MODULE(_aare, m) {
     define_3x3_reduction<double, 9, 9, uint16_t>(m);
     define_3x3_reduction<float, 9, 9, uint16_t>(m);
 
-    reduce_to_3x3<int, 3, 3, uint16_t>(m);
-    reduce_to_3x3<double, 3, 3, uint16_t>(m);
-    reduce_to_3x3<float, 3, 3, uint16_t>(m);
     reduce_to_3x3<int, 5, 5, uint16_t>(m);
     reduce_to_3x3<double, 5, 5, uint16_t>(m);
     reduce_to_3x3<float, 5, 5, uint16_t>(m);

python/src/module_config.hpp

@@ -0,0 +1,6 @@
+// SPDX-License-Identifier: MPL-2.0
+#pragma once
+#include <cstdint>
+
+// Configure module wide pedestal type for cluster finding
+using pd_type = double;

python/src/np_helper.hpp

@@ -41,31 +41,60 @@ template <typename T> py::array return_vector(std::vector<T> *vec) {
                           free_when_done); // numpy array references this parent
 }
 
-// todo rewrite generic
-template <class T, int Flags>
-auto get_shape_3d(const py::array_t<T, Flags> &arr) {
-    return aare::Shape<3>{arr.shape(0), arr.shape(1), arr.shape(2)};
+// Create a NDView of a numpy array. NDView assumes C-order strides so we
+// reject anything that is not Ndim dimensional and C-contiguous instead of
+// silently reading the data in the wrong order.
+template <ssize_t Ndim, class T, int Flags>
+auto get_checked_shape(const py::array_t<T, Flags> &arr) {
+    if (arr.ndim() != Ndim) {
+        throw py::value_error(
+            fmt::format("Expected {}D array, got {}D", Ndim, arr.ndim()));
+    }
+    if (!(arr.flags() & py::array::c_style)) {
+        throw py::value_error(
+            "Array is not C-contiguous. Use np.ascontiguousarray(arr)");
+    }
+    aare::Shape<Ndim> shape{};
+    for (ssize_t i = 0; i < Ndim; ++i) {
+        shape[i] = arr.shape(i);
+    }
+    return shape;
+}
+
+template <ssize_t Ndim, class T, int Flags>
+auto make_view(py::array_t<T, Flags> &arr) {
+    auto shape = get_checked_shape<Ndim>(arr);
+    return aare::NDView<T, Ndim>(arr.mutable_data(), shape);
+}
+
+// Read only view, also works for numpy arrays that are not writeable
+template <ssize_t Ndim, class T, int Flags>
+auto make_const_view(const py::array_t<T, Flags> &arr) {
+    auto shape = get_checked_shape<Ndim>(arr);
+    return aare::NDView<const T, Ndim>(arr.data(), shape);
 }
 
 template <class T, int Flags> auto make_view_3d(py::array_t<T, Flags> &arr) {
-    return aare::NDView<T, 3>(arr.mutable_data(), get_shape_3d<T, Flags>(arr));
+    return make_view<3>(arr);
 }
-
-template <class T, int Flags>
-auto get_shape_2d(const py::array_t<T, Flags> &arr) {
-    return aare::Shape<2>{arr.shape(0), arr.shape(1)};
-}
-
-template <class T, int Flags>
-auto get_shape_1d(const py::array_t<T, Flags> &arr) {
-    return aare::Shape<1>{arr.shape(0)};
-}
-
 template <class T, int Flags> auto make_view_2d(py::array_t<T, Flags> &arr) {
-    return aare::NDView<T, 2>(arr.mutable_data(), get_shape_2d<T, Flags>(arr));
+    return make_view<2>(arr);
 }
 template <class T, int Flags> auto make_view_1d(py::array_t<T, Flags> &arr) {
-    return aare::NDView<T, 1>(arr.mutable_data(), get_shape_1d<T, Flags>(arr));
+    return make_view<1>(arr);
+}
+
+template <class T, int Flags>
+auto make_const_view_3d(const py::array_t<T, Flags> &arr) {
+    return make_const_view<3>(arr);
+}
+template <class T, int Flags>
+auto make_const_view_2d(const py::array_t<T, Flags> &arr) {
+    return make_const_view<2>(arr);
+}
+template <class T, int Flags>
+auto make_const_view_1d(const py::array_t<T, Flags> &arr) {
+    return make_const_view<1>(arr);
 }
 
 template <typename ClusterType> struct fmt_format_trait; // forward declaration

python/src/pedestal.hpp

@@ -1,119 +0,0 @@
-// SPDX-License-Identifier: MPL-2.0
-
-#include "aare/Pedestal.hpp"
-#include "np_helper.hpp"
-
-#include <cstdint>
-#include <filesystem>
-#include <pybind11/numpy.h>
-#include <pybind11/pybind11.h>
-#include <pybind11/stl.h>
-
-namespace py = pybind11;
-
-template <typename SUM_TYPE>
-void define_pedestal_bindings(py::module &m, const std::string &name) {
-
-    py::class_<Pedestal<SUM_TYPE>>(m, name.c_str(), py::buffer_protocol())
-        .def(py::init<int, int, int>())
-        .def(py::init<int, int>())
-        .def("mean",
-             [](Pedestal<SUM_TYPE> &self) {
-                 auto mea = new NDArray<SUM_TYPE, 2>{};
-                 *mea = self.mean();
-                 return return_image_data(mea);
-             })
-        .def("view",
-             [](py::object self_py) {
-                 auto &self = self_py.cast<Pedestal<SUM_TYPE> &>();
-                 auto v = self.view();
-                 std::array<py::ssize_t, 2> shape{
-                     static_cast<py::ssize_t>(v.shape(0)),
-                     static_cast<py::ssize_t>(v.shape(1))};
-                 std::array<py::ssize_t, 2> byte_strides{
-                     static_cast<py::ssize_t>(v.strides()[0]) *
-                         static_cast<py::ssize_t>(sizeof(SUM_TYPE)),
-                     static_cast<py::ssize_t>(v.strides()[1]) *
-                         static_cast<py::ssize_t>(sizeof(SUM_TYPE))};
-                 auto arr = py::array_t<SUM_TYPE>(shape, byte_strides, v.data(),
-                                                  self_py);
-                 arr.attr("setflags")(py::arg("write") = false);
-                 return arr;
-             })
-        .def("variance",
-             [](Pedestal<SUM_TYPE> &self) {
-                 auto var = new NDArray<SUM_TYPE, 2>{};
-                 *var = self.variance();
-                 return return_image_data(var);
-             })
-        .def("std",
-             [](Pedestal<SUM_TYPE> &self) {
-                 auto std = new NDArray<SUM_TYPE, 2>{};
-                 *std = self.std();
-                 return return_image_data(std);
-             })
-        .def(
-            "__array_ufunc__",
-            [](py::object self, py::object ufunc, const std::string &method,
-               py::args inputs, py::kwargs kwargs) -> py::object {
-                if (method != "__call__" || inputs.size() != 2 ||
-                    inputs[1].ptr() != self.ptr() ||
-                    py::cast<std::string>(ufunc.attr("__name__")) !=
-                        "subtract") {
-                    return py::reinterpret_borrow<py::object>(
-                        Py_NotImplemented);
-                }
-
-                auto mean =
-                    py::module_::import("builtins").attr("memoryview")(self);
-                return ufunc(inputs[0], mean, **kwargs);
-            },
-            "Support subtracting a Pedestal from a NumPy array.")
-        .def("clear", py::overload_cast<>(&Pedestal<SUM_TYPE>::clear))
-        .def_property_readonly("rows", &Pedestal<SUM_TYPE>::rows)
-        .def_property_readonly("cols", &Pedestal<SUM_TYPE>::cols)
-        .def_property_readonly("n_samples", &Pedestal<SUM_TYPE>::n_samples)
-        .def_property_readonly("sum", &Pedestal<SUM_TYPE>::get_sum)
-        .def_property_readonly("sum2", &Pedestal<SUM_TYPE>::get_sum2)
-        .def("clone",
-             [&](Pedestal<SUM_TYPE> &pedestal) {
-                 return Pedestal<SUM_TYPE>(pedestal);
-             })
-        // TODO! add push for other data types
-        .def("push",
-             [](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> &f) {
-                 auto v = make_view_2d(f);
-                 pedestal.push(v);
-             })
-        .def(
-            "push_with_threshold",
-            [](Pedestal<SUM_TYPE> &pedestal,
-               py::array_t<uint16_t, py::array::c_style> &f,
-               py::array_t<SUM_TYPE, py::array::c_style> &threshold) {
-                auto frame_view = make_view_2d(f);
-                auto threshold_view = make_view_2d(threshold);
-                pedestal.push_with_threshold(frame_view, threshold_view);
-            },
-            py::arg("frame").noconvert(), py::arg("threshold").noconvert())
-        .def(
-            "push_no_update",
-            [](Pedestal<SUM_TYPE> &pedestal,
-               py::array_t<uint16_t, py::array::c_style> &f) {
-                auto v = make_view_2d(f);
-                pedestal.push_no_update(v);
-            },
-            py::arg().noconvert())
-        .def("update_mean", &Pedestal<SUM_TYPE>::update_mean)
-        .def_buffer([](Pedestal<SUM_TYPE> &self) {
-            auto mean = self.view();
-            return py::buffer_info(
-                const_cast<SUM_TYPE *>(mean.data()), sizeof(SUM_TYPE),
-                py::format_descriptor<SUM_TYPE>::format(), 2,
-                {static_cast<py::ssize_t>(mean.shape(0)),
-                 static_cast<py::ssize_t>(mean.shape(1))},
-                {static_cast<py::ssize_t>(mean.strides()[0] * sizeof(SUM_TYPE)),
-                 static_cast<py::ssize_t>(mean.strides()[1] *
-                                          sizeof(SUM_TYPE))},
-                true);
-        });
-}

python/src/raw_master_file.hpp

@@ -66,7 +66,8 @@ void define_raw_master_file_bindings(py::module &m) {
 
         .def_property_readonly("total_frames_expected",
                                &RawMasterFile::total_frames_expected)
-        .def_property_readonly("geometry", &RawMasterFile::geometry)
+        .def_property_readonly("detector_layout",
+                               &RawMasterFile::detector_layout)
         .def_property_readonly("udp_interfaces_per_module",
                                &RawMasterFile::udp_interfaces_per_module)
         .def_property_readonly("analog_samples", &RawMasterFile::analog_samples,
@@ -107,9 +108,46 @@ void define_raw_master_file_bindings(py::module &m) {
                     return std::nullopt;
                 }
             })
-        .def_property_readonly("period", [](RawMasterFile &self) {
-            double seconds =
-                std::chrono::duration<double>(self.period()).count();
-            return seconds;
-        });
+        .def_property_readonly("rois", &RawMasterFile::rois, R"(
+            Get the ROIs defined in the master file
+
+            Returns
+            ----------
+                List[ROI]
+                    List of ROIs (default complete ROI)
+            )")
+        .def_property_readonly("udp_port_types", &RawMasterFile::udp_port_types,
+                               R"(
+            Get the types of UDP ports
+
+            Returns
+            ----------
+                Optional[List[UDPPortPosition]]
+                    Optional vector of UDP port types as strings (only present for
+                    masterfile version >= 8.1)
+            )")
+        .def_property_readonly("disabled_udp_ports",
+                               &RawMasterFile::disabled_udp_ports, R"(
+            Get the indices of disabled UDP ports   
+
+            Returns
+            ----------
+                List[int]
+                    Vector of disabled UDP port indices relative to UDP port types (empty if none are disabled)
+            )")
+        .def_property_readonly(
+            "period",
+            [](RawMasterFile &self) {
+                double seconds =
+                    std::chrono::duration<double>(self.period()).count();
+                return seconds;
+            })
+
+        .def("__enter__",
+             [](RawMasterFile &self) -> RawMasterFile & { return self; })
+
+        .def("__exit__",
+             []([[maybe_unused]] RawMasterFile &self,
+                [[maybe_unused]] py::object, [[maybe_unused]] py::object,
+                [[maybe_unused]] py::object) { return; });
 }

python/src/var_cluster.hpp

@@ -21,40 +21,42 @@ void define_var_cluster_finder_bindings(py::module &m) {
                  auto *ptr = new NDArray<int, 2>(self.labeled());
                  return return_image_data(ptr);
              })
-        .def("set_noiseMap",
-             [](VarClusterFinder<double> &self,
-                py::array_t<double, py::array::c_style | py::array::forcecast>
-                    noise_map) {
-                 auto noise_map_span = make_view_2d(noise_map);
-                 self.set_noiseMap(noise_map_span);
-             })
+        .def(
+            "set_noiseMap", // conversions ok
+            [](VarClusterFinder<double> &self,
+               py::array_t<double, py::array::c_style | py::array::forcecast>
+                   noise_map) {
+                auto noise_map_span = make_view_2d(noise_map);
+                self.set_noiseMap(noise_map_span);
+            },
+            py::arg("noise_map"))
         .def("set_numberOfNeighbours",
              &VarClusterFinder<double>::set_numberOfNeighbours)
         .def("set_empty_surroundingPixels",
              &VarClusterFinder<double>::set_empty_surroundingPixels)
         .def("set_peripheralThresholdFactor",
              &VarClusterFinder<double>::set_peripheralThresholdFactor)
-        .def("find_clusters",
-             [](VarClusterFinder<double> &self,
-                py::array_t<double, py::array::c_style | py::array::forcecast>
-                    img) {
-                 auto view = make_view_2d(img);
-                 self.find_clusters(view);
-             })
-        .def("find_clusters_X",
-             [](VarClusterFinder<double> &self,
-                py::array_t<double, py::array::c_style | py::array::forcecast>
-                    img) {
-                 auto img_span = make_view_2d(img);
-                 self.find_clusters_X(img_span);
-             })
-        .def("single_pass",
-             [](VarClusterFinder<double> &self,
-                py::array_t<double, py::array::c_style | py::array::forcecast>
-                    img) {
-                 auto img_span = make_view_2d(img);
-                 self.single_pass(img_span);
-             })
+        .def(
+            "find_clusters",
+            [](VarClusterFinder<double> &self, py::array_t<double> img) {
+                auto view = make_view_2d(img);
+                self.find_clusters(view);
+            },
+            py::arg("img").noconvert())
+        .def(
+            "find_clusters_X",
+            [](VarClusterFinder<double> &self, py::array_t<double> img) {
+                auto img_span = make_view_2d(img);
+                self.find_clusters_X(img_span);
+            },
+            py::arg("img").noconvert())
+        .def(
+            "single_pass",
+            [](VarClusterFinder<double> &self, py::array_t<double> img) {
+                auto img_span = make_view_2d(img);
+                self.single_pass(img_span);
+            },
+            py::arg("img").noconvert())
         .def("hits",
              [](VarClusterFinder<double> &self) {
                  auto ptr = new std::vector<VarClusterFinder<double>::Hit>(

python/tests/test_Cluster.py

@@ -3,7 +3,7 @@ import pytest
 import numpy as np
 
 from aare import _aare #import the C++ module
-from aare import corner 
+from aare import ROI, corner 
 from conftest import test_data_path
 
 
@@ -91,6 +91,22 @@ def test_calculate_eta():
     assert eta2.sum == 4
 
 
+def test_cluster_representations():
+    """Test Cluster representation"""
+
+    cluster = _aare.Cluster2x2i(5, 6, np.array([1, 2, 3, 4], dtype=np.int32))
+    assert repr(cluster) == "Cluster2x2i(x=5, y=6, data=[[1, 2], [3, 4]])"
+
+def test_eta_representations():
+    """Test Eta representation"""
+    eta = _aare.Etai()
+    eta.x = 0.25
+    eta.y = 0.75
+    eta.c = corner.cBottomRight
+    eta.sum = 42
+    assert repr(eta) == "Etai(x=0.25, y=0.75, c=BottomRight, sum=42)"
+
+
 def test_max_sum(): 
     """Max 2x2 Sum"""
     cluster = _aare.Cluster3x3i(5,5,np.array([1, 1, 1, 2, 3, 1, 2, 2, 1], dtype=np.int32))
@@ -103,19 +119,59 @@ def test_max_sum():
 
 def test_cluster_finder(): 
     """Test ClusterFinder""" 
+    shape = [100,100]
+    cf = _aare.ClusterFinder_Cluster3x3i(shape)
 
-    clusterfinder = _aare.ClusterFinder_Cluster3x3i([100,100])
+    #Push 1000 frames to the pedestal
+    for i in range(1000):
+        frame = np.random.normal(loc = 100, scale = 5, size = shape).astype(np.uint16)
+        cf.push_pedestal_frame(frame)
+    cf.update_threshold()
+    frame = np.zeros(shape=shape, dtype=np.uint16)
+    cf.find_clusters(frame)
 
-    #frame = np.random.rand(100,100)
-    frame = np.zeros(shape=[100,100])
-
-    clusterfinder.find_clusters(frame)
-
-    clusters = clusterfinder.steal_clusters(False) #conversion does not work
+    clusters = cf.steal_clusters(False) #conversion does not work
 
     assert clusters.size == 0
 
 
+def test_cluster_finder_set_nSigma_before_pedestal_is_ready():
+    shape = [10, 10]
+    n_samples = 4
+    cf = _aare.ClusterFinder_Cluster3x3i(shape, min_pedestal_samples=n_samples)
+    assert not cf.pedestal_ready
+
+    cf.nSigma = 50
+    assert cf.nSigma == 50
+
+    # Pedestal of 100 with a noise of 2
+    for i in range(n_samples):
+        cf.push_pedestal_frame(np.full(shape, 98 + 4 * (i % 2), dtype=np.uint16))
+    assert cf.pedestal_ready
+
+    # 30 ADU is a cluster at the default nSigma of 5 but not at 50
+    frame = np.full(shape, 100, dtype=np.uint16)
+    frame[5, 5] = 130
+    cf.find_clusters(frame)
+    assert cf.steal_clusters().size == 0
+
+
+def test_cluster_finder_mt_set_nSigma_before_pedestal_is_ready():
+    shape = [10, 10]
+    n_samples = 4
+    cf = _aare.ClusterFinderMT_Cluster3x3i(
+        shape, n_threads=2, min_pedestal_samples=n_samples
+    )
+    assert not cf.pedestal_ready()
+
+    cf.set_nSigma(50)
+
+    for i in range(n_samples):
+        cf.push_pedestal_frame(np.full(shape, 98 + 4 * (i % 2), dtype=np.uint16))
+    assert cf.pedestal_ready()
+    cf.stop()
+
+
 def test_2x2_reduction(): 
     """Test 2x2 Reduction"""
     cluster = _aare.Cluster3x3i(5,5,np.array([1, 1, 1, 2, 3, 1, 2, 2, 1], dtype=np.int32))
@@ -144,3 +200,20 @@ def test_3x3_reduction():
 
 
 
+
+@pytest.mark.parametrize(
+    "cluster_type, n_pixels",
+    [(_aare.Cluster3x3i, 9), (_aare.Cluster2x2i, 4)],
+)
+def test_cluster_constructor_requires_exact_data_length(cluster_type, n_pixels):
+    """Data must hold exactly one value per pixel"""
+    for bad in (
+        np.ones(n_pixels - 1, dtype=np.int32),
+        np.ones(n_pixels + 1, dtype=np.int32),
+        np.ones((n_pixels, 1), dtype=np.int32),
+    ):
+        with pytest.raises(ValueError, match=f"{n_pixels} values"):
+            cluster_type(0, 0, bad)
+
+    cluster = cluster_type(0, 0, np.arange(n_pixels, dtype=np.int32))
+    assert (cluster.data.ravel() == np.arange(n_pixels)).all()

python/tests/test_ClusterFile.py

@@ -6,10 +6,279 @@ import boost_histogram as bh
 import time
 from pathlib import Path
 import pickle
+import struct
+import gc
+import weakref
 
-from aare import ClusterFile
+from aare import ClusterFile, ROI
 from conftest import test_data_path
 
+
+@pytest.fixture
+def iterator_file(tmp_path):
+    fname = tmp_path / "iteration.clust"
+    records = [
+        (-11, []),
+        (42, [(5, 6, 10), (7, 8, 20)]),
+        (103, []),
+        (104, [(5, 6, 30), (7, 8, 40)]),
+        (105, []),
+    ]
+    fname.write_bytes(b"".join(
+        struct.pack("@iI", number, len(clusters))
+        + b"".join(struct.pack("@HH9i", x, y, *([value] * 9))
+                   for x, y, value in clusters)
+        for number, clusters in records
+    ))
+    return fname
+
+
+def test_frames_preserve_empty_frames_and_frame_numbers(iterator_file):
+    with ClusterFile(iterator_file) as reader:
+        frames = reader.frames()
+        assert iter(frames) is frames
+        assert reader.tell() == 0
+        results = list(frames)
+        assert [frame.frame_number for frame in results] == [-11, 42, 103, 104, 105]
+        assert [frame.size for frame in results] == [0, 2, 0, 2, 0]
+        for _ in range(2):
+            with pytest.raises(StopIteration):
+                next(frames)
+        assert list(reader.frames()) == []
+
+
+@pytest.mark.parametrize("chunk_size, sizes", [(1, [1, 1, 1, 1]), (2, [2, 2]),
+                                             (3, [3, 1]), (10, [4])])
+@pytest.mark.parametrize("explicit", [False, True])
+def test_chunks_cross_frames(iterator_file, chunk_size, sizes, explicit):
+    with ClusterFile(iterator_file, chunk_size=chunk_size) as reader:
+        chunks = reader.chunks(chunk_size=chunk_size) if explicit else reader.chunks()
+        assert iter(chunks) is chunks
+        assert reader.tell() == 0
+        results = list(chunks)
+        assert [chunk.size for chunk in results] == sizes
+        data = np.concatenate([np.asarray(chunk) for chunk in results])
+        np.testing.assert_array_equal(data["x"], [5, 7, 5, 7])
+        np.testing.assert_array_equal(data["data"][:, 0, 0], [10, 20, 30, 40])
+        for _ in range(2):
+            with pytest.raises(StopIteration):
+                next(chunks)
+
+
+def test_chunk_override_preserves_default_iteration(iterator_file):
+    with ClusterFile(iterator_file, chunk_size=2) as reader:
+        assert iter(reader) is reader
+        assert next(reader.chunks(1)).size == 1
+        assert [chunk.size for chunk in reader] == [2, 1]
+
+
+@pytest.mark.parametrize("use_roi", [False, True])
+@pytest.mark.parametrize("use_noise", [False, True])
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+def test_iterators_apply_filters_and_gain(iterator_file, use_roi, use_noise, method):
+    with ClusterFile(iterator_file, chunk_size=3) as reader:
+        if use_roi:
+            reader.set_roi(ROI(0, 6, 0, 12))
+        if use_noise:
+            reader.set_noise_map(np.full((12, 12), 15, dtype=np.int32))
+        reader.set_gain_map(np.full((12, 12), 2.0))
+        results = list(getattr(reader, method)())
+        if method == "frames":
+            assert [frame.frame_number for frame in results] == [-11, 42, 103, 104, 105]
+        expected = [value / 2 for x, value in [(5, 10), (7, 20), (5, 30), (7, 40)]
+                    if (not use_roi or x == 5) and (not use_noise or value > 15)]
+        data = np.concatenate([np.asarray(result) for result in results])
+        np.testing.assert_array_equal(data["data"][:, 0, 0], expected)
+
+
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+def test_iterators_handle_all_rejected_clusters(iterator_file, method):
+    with ClusterFile(iterator_file) as reader:
+        reader.set_roi(ROI(0, 1, 0, 1))
+        results = list(getattr(reader, method)())
+        assert len(results) == (5 if method == "frames" else 0)
+        assert all(result.size == 0 for result in results)
+
+
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+def test_iterator_keeps_file_alive(iterator_file, method):
+    reader = ClusterFile(iterator_file)
+    owner = weakref.ref(reader)
+    iterator = getattr(reader, method)()
+    del reader
+    gc.collect()
+    assert owner() is not None
+    results = list(iterator)
+    assert sum(result.size for result in results) == 4
+    del iterator
+    gc.collect()
+    assert owner() is None
+
+
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+def test_retained_results_and_numpy_views_own_storage(iterator_file, method):
+    with ClusterFile(iterator_file, chunk_size=2) as reader:
+        reader.read_frame()  # Skip the initial empty frame.
+        iterator = getattr(reader, method)()
+        first = next(iterator)
+        array = np.asarray(first)
+        expected = array.copy()
+        results = list(iterator)
+        np.testing.assert_array_equal(np.asarray(first), expected)
+    del first, results, iterator, reader
+    gc.collect()
+    np.testing.assert_array_equal(array, expected)
+
+
+def test_frame_iteration_resumes_without_reading_ahead(iterator_file):
+    with ClusterFile(iterator_file) as reader:
+        reader.read_frame()
+        for frame in reader.frames():
+            assert frame.frame_number == 42
+            break
+        assert reader.read_frame().frame_number == 103
+        assert [frame.frame_number for frame in reader.frames()] == [104, 105]
+
+
+def test_partial_chunk_must_be_completed_before_frame_iteration(iterator_file):
+    with ClusterFile(iterator_file) as reader:
+        for chunk in reader.chunks(1):
+            assert chunk.size == 1
+            break
+        with pytest.raises(RuntimeError, match="clusters left"):
+            next(reader.frames())
+        assert reader.read_clusters(1).size == 1
+        assert [frame.frame_number for frame in reader.frames()] == [103, 104, 105]
+
+
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+def test_empty_iterators(tmp_path, method):
+    fname = tmp_path / "empty.clust"
+    fname.touch()
+    with ClusterFile(fname) as reader:
+        iterator = getattr(reader, method)()
+        for _ in range(2):
+            with pytest.raises(StopIteration):
+                next(iterator)
+
+
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+@pytest.mark.parametrize("started", [False, True])
+def test_iterators_reject_closed_files(iterator_file, method, started):
+    reader = ClusterFile(iterator_file, chunk_size=1)
+    iterator = getattr(reader, method)()
+    if started:
+        next(iterator)
+    reader.close()
+    with pytest.raises(RuntimeError, match="not opened for reading"):
+        next(iterator)
+
+
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+@pytest.mark.parametrize("mode", ["w", "a"])
+def test_iterators_reject_writing_modes(tmp_path, method, mode):
+    with ClusterFile(tmp_path / "output.clust", mode=mode) as reader:
+        with pytest.raises(RuntimeError, match="not opened for reading"):
+            next(getattr(reader, method)())
+
+
+def test_chunk_sizes_must_be_positive(iterator_file):
+    with ClusterFile(iterator_file, chunk_size=0) as reader:
+        with pytest.raises(ValueError, match="greater than zero"):
+            reader.chunks()
+        with pytest.raises(ValueError, match="greater than zero"):
+            reader.chunks(0)
+        with pytest.raises(ValueError, match="greater than zero"):
+            next(reader)
+        with pytest.raises(TypeError):
+            reader.chunks(-1)
+        assert reader.tell() == 0
+        assert reader.read_clusters(0).size == 0
+        assert len(list(reader.frames())) == 5
+
+
+@pytest.mark.parametrize("method", ["frames", "chunks"])
+def test_iterators_defer_read_errors_until_next_result(tmp_path, method):
+    fname = tmp_path / "truncated.clust"
+    cluster = struct.pack("@HH9i", 5, 6, *range(9))
+    frame = struct.pack("@iI", 42, 1) + cluster
+    fname.write_bytes(frame + frame[:-1])
+    with ClusterFile(fname, chunk_size=1) as reader:
+        iterator = getattr(reader, method)()
+        assert next(iterator).size == 1
+        assert reader.tell() == len(frame)
+        with pytest.raises(RuntimeError):
+            next(iterator)
+
+
+@pytest.mark.parametrize(
+    "shape, dtype",
+    [(shape, dtype)
+     for shape in [(2, 2), (3, 3), (5, 5), (7, 7), (9, 9)]
+     for dtype in [np.int32, np.float32, np.float64]]
+    + [((3, 3), np.int16)],
+)
+def test_iterators_are_bound_for_each_cluster_type(tmp_path, shape, dtype):
+    fname = tmp_path / "typed_frames.clust"
+    values = list(range(shape[0] * shape[1]))
+    record = struct.pack("@HH" + np.dtype(dtype).char * len(values), 5, 6, *values)
+    fname.write_bytes((struct.pack("@iI", -42, 1) + record) * 3)
+    with ClusterFile(fname, cluster_size=shape, dtype=dtype) as reader:
+        frame = next(reader.frames())
+        assert frame.frame_number == -42
+        np.testing.assert_array_equal(np.asarray(frame)["data"].reshape(-1), values)
+        chunks = list(reader.chunks())
+        assert len(chunks) == 1
+        assert chunks[0].size == 2
+        np.testing.assert_array_equal(np.asarray(chunks[0])["data"].reshape(-1), values * 2)
+
+
+def test_read_frame_returns_none_at_eof(tmp_path):
+    fname = tmp_path / "empty.clust"
+    fname.touch()
+
+    with ClusterFile(fname) as f:
+        assert f.read_frame() is None
+
+
+def test_read_frame_raises_for_malformed_file(tmp_path):
+    fname = tmp_path / "malformed.clust"
+    fname.write_bytes(b"\x00")
+
+    with ClusterFile(fname) as f, pytest.raises(RuntimeError):
+        f.read_frame()
+
+
+@pytest.mark.parametrize("use_roi", [False, True])
+@pytest.mark.parametrize("method", ["read_clusters", "default", "frames", "chunks"])
+@pytest.mark.parametrize(
+    "size",
+    [
+        pytest.param(3, id="partial-frame-number"),
+        pytest.param(4, id="missing-cluster-count"),
+        pytest.param(7, id="partial-cluster-count"),
+        pytest.param(48, id="missing-cluster-record"),
+        pytest.param(87, id="partial-cluster-record"),
+    ],
+)
+def test_chunk_reads_reject_incomplete_frames(tmp_path, use_roi, method, size):
+    fname = tmp_path / "incomplete.clust"
+    cluster = struct.pack("@HH9i", 5, 6, *range(9))
+    frame = struct.pack("@iI", 42, 2) + cluster * 2
+    fname.write_bytes(frame[:size])
+
+    with ClusterFile(fname) as reader:
+        if use_roi:
+            reader.set_roi(ROI(0, 10, 0, 10))
+        with pytest.raises(RuntimeError):
+            if method == "default":
+                next(reader)
+            elif method == "read_clusters":
+                reader.read_clusters(10)
+            else:
+                next(getattr(reader, method)())
+
+
 @pytest.mark.withdata
 def test_cluster_file(test_data_path): 
     """Test ClusterFile""" 

python/tests/test_ClusterVector.py

@@ -83,12 +83,15 @@ def test_make_a_hitmap_from_cluster_vector():
 
 def test_2x2_reduction(): 
     cv = ClusterVector((3,3))
+    cv.frame_number = -135
 
     cv.push_back(_aare.Cluster3x3i(5, 5, np.array([1, 1, 1, 2, 3, 1, 2, 2, 1], dtype=np.int32)))
     cv.push_back(_aare.Cluster3x3i(5, 5, np.array([2, 2, 1, 2, 3, 1, 1, 1, 1], dtype=np.int32)))
 
-    reduced_cv = np.array(_aare.reduce_to_2x2(cv), copy=False) 
+    reduced = _aare.reduce_to_2x2(cv)
+    reduced_cv = np.array(reduced, copy=False)
 
+    assert reduced.frame_number == cv.frame_number
     assert reduced_cv.size == 2
     assert reduced_cv[0]["x"] == 5
     assert reduced_cv[0]["y"] == 5
@@ -100,14 +103,17 @@ def test_2x2_reduction():
     
 def test_3x3_reduction(): 
     cv = _aare.ClusterVector_Cluster5x5d()
+    cv.frame_number = 246
     
     cv.push_back(_aare.Cluster5x5d(5,5,np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 2.0, 2.0, 3.0,
                                    1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], dtype=np.double)))
     cv.push_back(_aare.Cluster5x5d(5,5,np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 2.0, 2.0, 3.0,
                                    1.0, 1.0, 1.0, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0], dtype=np.double)))
     
-    reduced_cv = np.array(_aare.reduce_to_3x3(cv), copy=False)  
+    reduced = _aare.reduce_to_3x3(cv)
+    reduced_cv = np.array(reduced, copy=False)
 
+    assert reduced.frame_number == cv.frame_number
     assert reduced_cv.size == 2
     assert reduced_cv[0]["x"] == 5
     assert reduced_cv[0]["y"] == 5
@@ -130,3 +136,26 @@ def test_masking():
     assert cv_masked_array[0]["x"] == 1
     assert cv_masked_array[0]["y"] == 2
     assert (cv_masked_array[0]["data"] == np.ones((3,3),dtype=np.int32)).all()
+
+
+def test_masking_requires_c_contiguous_array():
+    cv = _aare.ClusterVector_Cluster3x3i()
+    cv.push_back(_aare.Cluster3x3i(1, 2, np.ones(9, dtype=np.int32)))
+    cv.push_back(_aare.Cluster3x3i(3, 4, np.ones(9, dtype=np.int32)))
+
+    mask = np.array([True, False, True, False], dtype=bool)[::2]
+    assert not mask.flags.c_contiguous
+
+    with pytest.raises(ValueError, match="C-contiguous"):
+        cv(mask)
+
+
+def test_masking_requires_one_dimension():
+    cv = _aare.ClusterVector_Cluster3x3i()
+    cv.push_back(_aare.Cluster3x3i(1, 2, np.ones(9, dtype=np.int32)))
+    cv.push_back(_aare.Cluster3x3i(3, 4, np.ones(9, dtype=np.int32)))
+
+    mask = np.array([[True, False]], dtype=bool)
+
+    with pytest.raises(ValueError):
+        cv(mask)

python/tests/test_CtbRawFile.py

@@ -0,0 +1,59 @@
+# SPDX-License-Identifier: MPL-2.0
+
+
+import pytest
+import numpy as np
+from aare import CtbRawFile
+
+
+@pytest.mark.withdata
+def test_basic_properites(test_data_path):
+    f = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+    assert f.total_frames == 47
+    assert f.image_size_in_bytes == 640
+    assert f.master.analog_samples == 10
+    assert f.master.digital_samples is None
+    assert f.master.transceiver_samples is None
+
+
+@pytest.mark.withdata
+def test_seek_and_read(test_data_path):
+    f = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+    f.seek(0)
+    header, raw_data = f.read_frame()
+    assert header['frameNumber'] == 48 #we know that the first frame in this file is 48
+
+    #Seek a multiple of frames in file
+    f.seek(9)
+    header, raw_data = f.read_frame()
+    assert header['frameNumber'] == 57 
+
+@pytest.mark.withdata
+def test_read_all_frames_direct_and_iterate(test_data_path):
+    """
+    We want to make sure that iterating and reading give the same result.
+    """
+    f0 = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+    f1 = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+
+    for h0, d0 in f0:
+        h1, d1 = f1.read_frame()
+        assert h0 == h1
+        assert np.all(d0 == d1)
+
+@pytest.mark.withdata
+def test_read_all_frames_with_index_and_iterate(test_data_path):
+    """
+    We want to make sure that iterating and reading give the same result.
+    """
+    f0 = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+    f1 = CtbRawFile(test_data_path/'raw/ctb/run_master_0.json')
+
+    frames_read = 0
+    for i,(h0, d0) in enumerate(f0):
+        h1, d1 = f1.read_frame(i)
+        assert h0 == h1
+        assert np.all(d0 == d1)
+        frames_read += 1
+
+    assert frames_read == f0.total_frames

python/tests/test_FastPedestal.py

@@ -0,0 +1,160 @@
+# SPDX-License-Identifier: MPL-2.0
+import numpy as np
+import pytest
+
+from aare import (
+    FastPedestal,
+    FastPedestal_d,
+    FastPedestal_f,
+    FastPedestal_i16,
+)
+
+
+@pytest.mark.parametrize(
+    ("dtype", "pedestal_type"),
+    [
+        (np.float64, FastPedestal_d),
+        (np.float32, FastPedestal_f),
+        (np.int16, FastPedestal_i16),
+    ],
+)
+def test_fast_pedestal_factory(dtype, pedestal_type):
+    pedestal = FastPedestal(2, 3, n_samples=4, dtype=dtype)
+
+    assert isinstance(pedestal, pedestal_type)
+    assert pedestal.rows == 2
+    assert pedestal.cols == 3
+    assert pedestal.n_samples == 4
+
+
+def test_fast_pedestal_factory_defaults_to_double():
+    assert isinstance(FastPedestal(2, 3), FastPedestal_d)
+
+
+def test_fast_pedestal_factory_rejects_unbound_dtype():
+    with pytest.raises(ValueError, match="Unsupported dtype for FastPedestal"):
+        FastPedestal(2, 3, dtype=np.int32)
+
+
+@pytest.mark.parametrize(
+    ("kwargs", "expected_n_samples"),
+    [
+        ({"rows": 2, "cols": 3}, 1000),
+        ({"rows": 2, "cols": 3, "n_samples": 4}, 4),
+    ],
+)
+def test_fast_pedestal_binding_accepts_constructor_keywords(
+    kwargs, expected_n_samples
+):
+    pedestal = FastPedestal_d(**kwargs)
+
+    assert pedestal.rows == 2
+    assert pedestal.cols == 3
+    assert pedestal.n_samples == expected_n_samples
+
+
+@pytest.mark.parametrize(
+    ("dtype", "pedestal_type", "expected_dtype"),
+    [
+        (np.float64, FastPedestal_d, np.float64),
+        (np.float32, FastPedestal_f, np.float32),
+        (np.int16, FastPedestal_i16, np.int16),
+    ],
+)
+def test_fast_pedestal_factory_from_file(
+    tmp_path, dtype, pedestal_type, expected_dtype
+):
+    frames = np.array(
+        [[[100, 100]], [[2, 4]], [[4, 6]], [[5, 7]]], dtype=np.uint16
+    )
+    filename = tmp_path / "frames.npy"
+    np.save(filename, frames)
+
+    pedestal = FastPedestal.from_file(
+        filename, n_samples=2, skip_first=1, dtype=dtype
+    )
+
+    assert isinstance(pedestal, pedestal_type)
+    assert pedestal.ready
+    assert pedestal.cur_samples == 2
+    assert pedestal.mean().dtype == expected_dtype
+    np.testing.assert_array_equal(pedestal.mean(), [[4, 6]])
+
+
+def test_fast_pedestal_factory_from_file_rejects_unbound_dtype():
+    with pytest.raises(ValueError, match="Unsupported dtype for FastPedestal"):
+        FastPedestal.from_file("unused.npy", dtype=np.int32)
+
+
+def test_fast_pedestal_from_file_rejects_skip_beyond_end(tmp_path):
+    filename = tmp_path / "frames.npy"
+    np.save(filename, np.zeros((1, 1, 1), dtype=np.uint16))
+
+    with pytest.raises(RuntimeError, match="less frames"):
+        FastPedestal.from_file(filename, n_samples=1, skip_first=2)
+
+
+@pytest.mark.parametrize(
+    ("pedestal_type", "expected_dtype"),
+    [(FastPedestal_d, np.float64), (FastPedestal_f, np.float32)],
+)
+def test_fast_pedestal_initialization(pedestal_type, expected_dtype):
+    pedestal = pedestal_type(2, 3, 2)
+    first = np.array([[2, 4, 6], [8, 10, 12]], dtype=np.uint16)
+    second = np.array([[4, 6, 8], [10, 12, 14]], dtype=np.uint16)
+
+    pedestal.add_init_frame(first)
+    pedestal.add_init_frame(second)
+
+
+    expected_mean = np.array(
+        [[3, 5, 7], [9, 11, 13]], dtype=expected_dtype
+    )
+    np.testing.assert_array_equal(pedestal.mean(), expected_mean)
+    np.testing.assert_array_equal(pedestal.std(), np.ones((2, 3)))
+
+
+def test_fast_pedestal_steady_state_push_ema():
+    pedestal = FastPedestal_d(1, 2, 2)
+    pedestal.add_init_frame(np.array([[2, 4]], dtype=np.uint16))
+    pedestal.add_init_frame(np.array([[4, 6]], dtype=np.uint16))
+
+
+    pedestal.push_ema(np.array([[6, 8]], dtype=np.uint16))
+
+    np.testing.assert_array_equal(pedestal.mean(), [[4.5, 6.5]])
+
+
+def test_fast_pedestal_exposes_read_only_buffer_and_subtraction():
+    pedestal = FastPedestal_d(1, 2, 1)
+    pedestal.add_init_frame(np.array([[2, 4]], dtype=np.uint16))
+
+
+    view = np.asarray(pedestal)
+    result = np.array([[12, 14]], dtype=np.uint16) - pedestal
+
+    np.testing.assert_array_equal(view, [[2, 4]])
+    np.testing.assert_array_equal(result, [[10, 10]])
+    assert np.shares_memory(view, pedestal.view())
+    assert not view.flags.writeable
+
+
+def test_fast_pedestal_rejects_wrong_shape():
+    pedestal = FastPedestal_d(2, 3)
+
+    with pytest.raises(RuntimeError, match="shape"):
+        pedestal.add_init_frame(np.zeros((2, 2), dtype=np.uint16))
+
+
+@pytest.mark.parametrize("dtype", [np.float64, np.float32, np.int16])
+def test_fast_pedestal_std_uses_double_variance(dtype):
+    pedestal = FastPedestal(1, 1, n_samples=2, dtype=dtype)
+    pedestal.add_init_frame(np.array([[0]], dtype=np.uint16))
+    pedestal.add_init_frame(np.array([[1000]], dtype=np.uint16))
+
+    assert pedestal.std().dtype == dtype
+    np.testing.assert_array_equal(pedestal.std(), [[500]])
+
+    pedestal.push_ema(np.array([[1000]], dtype=np.uint16))
+    expected = np.array([[np.sqrt(187500.0)]], dtype=dtype)
+    np.testing.assert_array_equal(pedestal.std(), expected)

python/tests/test_Fit.py

@@ -1,15 +1,401 @@
 # SPDX-License-Identifier: MPL-2.0
 import numpy as np
+import pytest
+
 import aare
 
+MINIMIZERS = [
+    aare.Minimizer.Migrad,
+    aare.Minimizer.Fumili,
+    aare.Minimizer.LevenbergMarquardt,
+    aare.Minimizer.VarPro,
+]
 
 
-def test_gaussian_model_evaluates_and_fits_data():
+def gaussian(x, A, mu, sigma):
+    return A * np.exp(-0.5 * ((x - mu) / sigma) ** 2)
+
+
+def uses_minuit2(minimizer):
+    return minimizer not in (
+        aare.Minimizer.LevenbergMarquardt,
+        aare.Minimizer.VarPro,
+    )
+
+
+@pytest.fixture(params=MINIMIZERS, ids=lambda m: m.name)
+def minimizer(request):
+    return request.param
+
+
+@pytest.fixture
+def gaussian_data():
+    x = np.linspace(-6.0, 6.0, 61)
+    y = gaussian(x, 120.0, 0.8, 1.3)
+    y_err = np.ones_like(y)
+    return x, y, y_err
+
+
+def test_gaussian_model_evaluates_and_fits_data(minimizer):
     x = np.linspace(-5.0, 5.0, 51)
     expected = np.array([20.0, 0.5, 1.2])
-    model = aare.Gaussian()
+    model = aare.Gaussian(minimizer=minimizer)
 
     y = model(x, expected)
     result = model.fit(x, y)
 
     np.testing.assert_allclose(result["par"], expected, atol=2e-3)
+
+
+def test_minimizer_can_be_selected():
+    assert aare.Gaussian().minimizer == aare.Minimizer.Migrad
+    assert set(aare.Minimizer.__members__) == {
+        "Migrad",
+        "Fumili",
+        "LevenbergMarquardt",
+        "VarPro",
+    }
+
+    for minimizer in MINIMIZERS:
+        assert aare.Gaussian(minimizer=minimizer).minimizer == minimizer
+
+    model = aare.Gaussian()
+    model.minimizer = aare.Minimizer.LevenbergMarquardt
+    assert model.minimizer == aare.Minimizer.LevenbergMarquardt
+
+
+def test_fit_returns_errors_and_chi2(gaussian_data, minimizer):
+    x, y, y_err = gaussian_data
+    model = aare.Gaussian(compute_errors=True, minimizer=minimizer)
+    result = model.fit(x, y, y_err)
+
+    assert set(result) == {"par", "par_err", "chi2"}
+    assert result["par"].shape == (3,)
+    assert result["par_err"].shape == (3,)
+    assert result["chi2"].shape == (1,)
+    np.testing.assert_allclose(result["par"], [120.0, 0.8, 1.3], rtol=1e-4)
+    np.testing.assert_allclose(
+        result["par_err"], [0.361, 0.00451, 0.00451], rtol=2e-2
+    )
+    assert result["chi2"][0] == pytest.approx(0.0, abs=1e-3)
+
+    without_errors = aare.fit(aare.Gaussian(minimizer=minimizer), x, y)
+    assert set(without_errors) == {"par", "chi2"}
+
+
+def test_pol1_weighted_errors_match_analytic_formula(minimizer):
+    x = np.arange(20, dtype=float)
+    y = 2.0 + 0.5 * x
+    y_err = 1.0 + 0.1 * x
+    w = 1.0 / y_err**2
+    S, Sx, Sxx = w.sum(), (w * x).sum(), (w * x * x).sum()
+    delta = S * Sxx - Sx * Sx
+
+    result = aare.Pol1(compute_errors=True, minimizer=minimizer).fit(x, y, y_err)
+
+    np.testing.assert_allclose(result["par"], [2.0, 0.5], rtol=1e-6)
+    np.testing.assert_allclose(
+        result["par_err"],
+        [np.sqrt(Sxx / delta), np.sqrt(S / delta)],
+        rtol=1e-6,
+    )
+
+
+def test_fixed_parameter_reports_zero_error(gaussian_data, minimizer):
+    x, y, y_err = gaussian_data
+    model = aare.Gaussian(compute_errors=True, minimizer=minimizer)
+    model.FixParameter("sigma", 1.5)
+    result = model.fit(x, y, y_err)
+
+    assert result["par"][2] == 1.5
+    assert result["par_err"][2] == 0.0
+    assert result["par"][0] == pytest.approx(111.15, rel=1e-3)
+    assert result["chi2"][0] == pytest.approx(1684.25, rel=1e-3)
+
+    model.ReleaseParameter("sigma")
+    result = model.fit(x, y, y_err)
+    assert result["par"][2] == pytest.approx(1.3, rel=1e-4)
+
+
+def test_active_limit_ends_on_the_limit(gaussian_data, minimizer):
+    x, y, y_err = gaussian_data
+    model = aare.Gaussian(compute_errors=True, minimizer=minimizer)
+    model.SetParLimits("mu", -5.0, 0.5)
+    result = model.fit(x, y, y_err)
+
+    margin = 1e-6 if uses_minuit2(minimizer) else 1e-9
+    assert result["par"][1] == pytest.approx(0.5, abs=margin)
+    assert result["par_err"][1] == 0.0
+    assert result["par"][0] == pytest.approx(116.88, rel=1e-3)
+    assert result["chi2"][0] == pytest.approx(4301.95, rel=1e-3)
+
+    inactive = aare.Gaussian(minimizer=minimizer)
+    inactive.SetParLimits(1, -5.0, 5.0)
+    np.testing.assert_allclose(
+        inactive.fit(x, y)["par"], [120.0, 0.8, 1.3], rtol=1e-4
+    )
+
+
+def test_user_start_value(gaussian_data, minimizer):
+    x, y, _ = gaussian_data
+    model = aare.Gaussian(minimizer=minimizer)
+    model.SetParameter("mu", 0.0)
+    np.testing.assert_allclose(model.fit(x, y)["par"], [120.0, 0.8, 1.3], rtol=1e-4)
+
+
+def test_max_calls_failure_returns_zeros(gaussian_data, minimizer):
+    x, y, y_err = gaussian_data
+    model = aare.Gaussian(max_calls=3, compute_errors=True, minimizer=minimizer)
+    # start far from the minimum so that no minimizer converges in 3 calls
+    model.SetParameter("mu", -3.0)
+    model.SetParameter("sigma", 3.0)
+    result = model.fit(x, y, y_err)
+    assert np.all(result["par"] == 0.0)
+    assert np.all(result["par_err"] == 0.0)
+    assert np.all(result["chi2"] == 0.0)
+
+
+def test_strategy_is_ignored_by_levenberg_marquardt(gaussian_data):
+    x, y, _ = gaussian_data
+    lm = aare.Minimizer.LevenbergMarquardt
+    a = aare.Gaussian(strategy=0, minimizer=lm).fit(x, y)
+    b = aare.Gaussian(strategy=1, minimizer=lm).fit(x, y)
+    np.testing.assert_array_equal(a["par"], b["par"])
+
+
+def test_invalid_parameter_access_raises():
+    model = aare.Gaussian()
+    assert model.par_names == ["A", "mu", "sigma"]
+    assert model.n_par == 3
+    with pytest.raises(RuntimeError):
+        model.SetParameter("nope", 1.0)
+    with pytest.raises(IndexError):
+        model.SetParameter(3, 1.0)
+    with pytest.raises(RuntimeError):
+        model.SetParLimits("mu", 1.0, 1.0)
+
+
+def test_mismatched_array_sizes_raise(gaussian_data, minimizer):
+    x, y, y_err = gaussian_data
+    model = aare.Gaussian(minimizer=minimizer)
+    with pytest.raises(RuntimeError):
+        model.fit(x, y[:-1])
+    with pytest.raises(RuntimeError):
+        model.fit(x, y, y_err[:-1])
+
+
+@pytest.mark.parametrize("cls", [aare.RisingScurve, aare.FallingScurve])
+def test_scurve_fit(cls, minimizer):
+    truth = np.array([10.0, 0.1, 50.0, 3.0, 1000.0, 2.0])
+    x = np.linspace(0.0, 100.0, 100)
+    model = cls(max_calls=500, compute_errors=True, minimizer=minimizer)
+    y = model(x, truth)
+    result = model.fit(x, y, np.ones_like(y))
+
+    np.testing.assert_allclose(result["par"], truth, rtol=1e-3)
+    assert np.all(result["par_err"] > 0.0)
+
+
+@pytest.mark.parametrize("with_errors", [False, True])
+@pytest.mark.parametrize("weighted", [False, True])
+def test_fit_3d_matches_fit_1d(with_errors, weighted, minimizer):
+    x = np.linspace(-6.0, 6.0, 61)
+    rows, cols = 4, 3
+    y = np.empty((rows, cols, x.size))
+    for r in range(rows):
+        for c in range(cols):
+            y[r, c] = gaussian(x, 80.0 + 10.0 * r, -0.6 + 0.3 * c, 0.9 + 0.1 * r)
+    y_err = np.ones_like(y) if weighted else None
+
+    model = aare.Gaussian(compute_errors=with_errors, minimizer=minimizer)
+    cube = model.fit(x, y, y_err, n_threads=3)
+
+    expected_keys = {"par", "par_err", "chi2"} if with_errors else {"par", "chi2"}
+    assert set(cube) == expected_keys
+    assert cube["par"].shape == (rows, cols, 3)
+    assert cube["chi2"].shape == (rows, cols)
+    if with_errors:
+        assert cube["par_err"].shape == (rows, cols, 3)
+    for r in range(rows):
+        for c in range(cols):
+            single = model.fit(x, y[r, c], None if y_err is None else y_err[r, c])
+            np.testing.assert_allclose(cube["par"][r, c], single["par"], rtol=1e-10)
+            assert cube["chi2"][r, c] == pytest.approx(single["chi2"][0], abs=1e-10)
+            if with_errors:
+                np.testing.assert_allclose(
+                    cube["par_err"][r, c], single["par_err"], rtol=1e-10
+                )
+    np.testing.assert_allclose(cube["par"][3, 2], [110.0, 0.0, 1.2], rtol=1e-4, atol=1e-6)
+
+
+def test_all_minimizers_agree_on_3d_data_with_errors():
+    x = np.linspace(-5.0, 5.0, 51)
+    expected = np.array([20.0, 0.5, 1.2])
+    single = aare.Gaussian()(x, expected)
+    y = np.tile(single, (2, 3, 1))
+    y_err = np.ones_like(y)
+
+    results = {
+        minimizer: aare.Gaussian(compute_errors=True, minimizer=minimizer).fit(
+            x, y, y_err, n_threads=2
+        )
+        for minimizer in MINIMIZERS
+    }
+    reference = results[aare.Minimizer.Migrad]
+    for minimizer, result in results.items():
+        assert result["par"].shape == (2, 3, 3)
+        assert result["par_err"].shape == (2, 3, 3)
+        np.testing.assert_allclose(result["par"], np.broadcast_to(expected, (2, 3, 3)), atol=2e-3)
+        np.testing.assert_allclose(result["chi2"], 0.0, atol=1e-3)
+        np.testing.assert_allclose(result["par_err"], reference["par_err"], rtol=2e-2)
+
+
+# Regressions of the built-in solvers: fits that used to be reported as valid
+# although they had stopped away from the minimum or returned garbage.
+
+
+def test_gaussian_fit_to_noise_fails_or_stays_bounded(minimizer):
+    # A width collapsing between two scan points let the amplitude diverge
+    # (up to 1e159) in a fit that was still reported as valid.
+    x = np.linspace(-6.0, 6.0, 61)
+    y = np.random.default_rng(0).normal(0.0, 1.0, (20, 20, x.size))
+    result = aare.Gaussian(minimizer=minimizer).fit(x, y, n_threads=2)
+
+    par = result["par"]
+    failed = np.all(par == 0.0, axis=-1)
+    bounded = np.all(np.isfinite(par) & (np.abs(par) < 1e6), axis=-1)
+    assert np.all(failed | bounded)
+    assert np.all(np.isfinite(result["chi2"]))
+
+
+def test_step_reflected_onto_start_is_not_convergence(minimizer):
+    # The step of +2 reflects at the limit 1 back onto the start 0; the
+    # truncated step must still be tried.
+    model = aare.Pol1(minimizer=minimizer)
+    model.SetParameter("p0", 0.0)
+    model.FixParameter("p1", 0.0)
+    model.SetParLimits("p0", -1.0, 1.0)
+    result = model.fit(np.arange(4.0), np.full(4, 2.2))
+
+    assert result["par"][0] == pytest.approx(1.0, abs=1e-6)
+    assert result["chi2"][0] == pytest.approx(5.76, rel=1e-6)
+
+
+@pytest.mark.parametrize(
+    "amplitude, offset", [(1e8, 0.0), (1e10, 0.0), (120.0, 1e9)]
+)
+def test_large_parameter_does_not_stop_the_fit(minimizer, amplitude, offset):
+    # One relative step norm over all parameters let a large amplitude or
+    # position stop the fit at the start estimate of the width.
+    x = np.linspace(-6.0, 6.0, 61)
+    y = gaussian(x, amplitude, 0.8, 1.3)
+    y_err = np.sqrt(y) + 1.0
+    result = aare.Gaussian(minimizer=minimizer, max_calls=500).fit(
+        x + offset, y, y_err
+    )
+
+    par = result["par"]
+    assert par[0] == pytest.approx(amplitude, rel=1e-3)
+    assert par[1] - offset == pytest.approx(0.8, abs=1e-3)
+    assert par[2] == pytest.approx(1.3, rel=1e-3)
+
+
+def test_rejected_trial_on_limit_does_not_freeze_interior_parameter(minimizer):
+    # A rejected trial truncated onto the lower limit froze mu at its start.
+    x = np.linspace(-0.1, 0.1, 11)
+    mu = np.sqrt(2.0 * np.log(2.0))
+    y = gaussian(x, 10.0, mu, 1.0)
+    model = aare.Gaussian(minimizer=minimizer)
+    model.FixParameter("A", 10.0)
+    model.FixParameter("sigma", 1.0)
+    model.SetParameter("mu", 3.0)
+    model.SetParLimits("mu", 0.0, 5.0)
+    result = model.fit(x, y)
+
+    assert result["par"][1] == pytest.approx(mu, abs=1e-3)
+    assert result["chi2"][0] < 1e-4
+
+
+# Samples with a zero error are ignored, also when they hold NaN (a masked
+# dead frame): a zero weight does not remove a NaN from a sum.
+
+
+@pytest.mark.parametrize("masked", [4, 0, 9])
+def test_masked_nan_sample_is_ignored(minimizer, masked):
+    # The start estimate of Pol1 uses the first and the last sample.
+    x = np.arange(10.0)
+    y = 2.0 + 0.5 * x
+    y_err = np.ones_like(y)
+    y[masked] = np.nan
+    y_err[masked] = 0.0
+    result = aare.Pol1(minimizer=minimizer, compute_errors=True).fit(x, y, y_err)
+
+    np.testing.assert_allclose(result["par"], [2.0, 0.5], atol=1e-5)
+    assert np.all(result["par_err"] > 0.0)
+    assert result["chi2"][0] == pytest.approx(0.0, abs=1e-6)
+
+
+def test_too_few_unmasked_samples_fail(minimizer):
+    x = np.arange(10.0)
+    y = np.full_like(x, np.nan)
+    y_err = np.zeros_like(x)
+    y[0], y_err[0] = 2.0, 1.0  # one sample left for two free parameters
+    result = aare.Pol1(minimizer=minimizer).fit(x, y, y_err)
+    assert np.all(result["par"] == 0.0)
+    assert np.all(result["chi2"] == 0.0)
+
+
+def test_masked_samples_in_cube_match_1d_fits(minimizer):
+    x = np.linspace(-6.0, 6.0, 61)
+    y = np.stack([gaussian(x, 80.0 + 10.0 * c, 0.8, 1.2) for c in range(4)])
+    y = y.reshape(1, 4, x.size)
+    y_err = np.ones_like(y)
+    for col, masked in [(1, 10), (2, 0), (3, 34)]:
+        y[0, col, masked] = np.nan
+        y_err[0, col, masked] = 0.0
+    model = aare.Gaussian(minimizer=minimizer, compute_errors=True)
+    cube = model.fit(x, y, y_err, n_threads=1)
+
+    for col in range(4):
+        single = model.fit(x, y[0, col], y_err[0, col])
+        np.testing.assert_allclose(cube["par"][0, col], single["par"], atol=1e-12)
+        np.testing.assert_allclose(cube["par_err"][0, col], single["par_err"], atol=1e-12)
+        np.testing.assert_allclose(
+            cube["par"][0, col], [80.0 + 10.0 * col, 0.8, 1.2], rtol=1e-3
+        )
+
+
+def test_parameter_whose_minimum_lies_on_its_limit_reports_no_error(minimizer):
+    # At a minimum that lies exactly on the limit the gradient is zero.
+    x = np.array([-1.0, -1.0, 1.0, 1.0])
+    model = aare.Pol1(minimizer=minimizer, compute_errors=True)
+    model.SetParLimits("p0", 0.0, 1.0)
+    result = model.fit(x, 1.0 + 2.0 * x)
+
+    assert result["par"][0] == pytest.approx(1.0, abs=1e-6)
+    assert result["par_err"][0] == 0.0
+    assert result["par_err"][1] == pytest.approx(0.5, rel=1e-6)
+
+
+@pytest.mark.parametrize("amplitude", [0.1, 1e-4])
+def test_fit_does_not_stop_at_start_without_a_step(minimizer, amplitude):
+    # The EDM tolerance is absolute in chi2 units: an unweighted fit of small
+    # data passes it at the start estimate (sigma 1.3617) already.
+    x = np.linspace(-6.0, 6.0, 61)
+    y = gaussian(x, amplitude, 0.8, 1.3)
+    result = aare.Gaussian(minimizer=minimizer).fit(x, y)
+    assert abs(result["par"][2] - 1.3) < 0.03  # Migrad ends at 1.3205
+
+    one_call = aare.Gaussian(minimizer=minimizer, max_calls=1).fit(x, y)
+    assert np.all(one_call["par"] == 0.0)
+
+
+@pytest.mark.parametrize("lo, span", [(8000, 1000), (20000, 500), (60000, 200)])
+def test_polynomial_far_from_zero_keeps_its_precision(minimizer, lo, span):
+    # Nearly collinear columns 1, x, x^2: VarPro's normal equations lost six
+    # digits at x = 20000.
+    x = np.linspace(lo, lo + span, 40)
+    c = np.array([1.5 - 2e-3 * lo + 3e-7 * lo**2, 2e-3 - 6e-7 * lo, 3e-7])
+    y = c[0] + c[1] * x + c[2] * x**2
+    result = aare.Pol2(minimizer=minimizer).fit(x, y)
+    np.testing.assert_allclose(result["par"], c, rtol=1e-9)

python/tests/test_InterpolationAPI.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 import pytest 
 
 from aare import Interpolator, ClusterVector, Etai, Cluster

python/tests/test_Pedestal.py

@@ -1,7 +1,107 @@
+# SPDX-License-Identifier: MPL-2.0
 import numpy as np
 import pytest
 
-from aare import Pedestal_d, Pedestal_f
+from aare import Pedestal, Pedestal_d, Pedestal_f, Pedestal_i16
+
+
+@pytest.mark.parametrize(
+    ("dtype", "pedestal_type"),
+    [(np.float64, Pedestal_d), (np.float32, Pedestal_f), (np.int16, Pedestal_i16)],
+)
+def test_pedestal_factory(dtype, pedestal_type):
+    pedestal = Pedestal(rows=2, cols=3, n_samples=4, dtype=dtype)
+
+    assert isinstance(pedestal, pedestal_type)
+    assert pedestal.rows == 2
+    assert pedestal.cols == 3
+    assert pedestal.n_samples == 4
+
+
+def test_pedestal_factory_defaults_to_double():
+    pedestal = Pedestal(2, 3)
+
+    assert isinstance(pedestal, Pedestal_d)
+    assert pedestal.n_samples == 1000
+
+
+def test_pedestal_factory_rejects_unbound_dtype():
+    with pytest.raises(ValueError, match="Unsupported dtype for Pedestal"):
+        Pedestal(2, 3, dtype=np.int32)
+
+
+@pytest.mark.parametrize("dtype", [np.float64, np.float32, np.int16])
+@pytest.mark.parametrize("name", ["rows", "cols", "n_samples"])
+def test_pedestal_factory_rejects_negative_parameters(dtype, name):
+    parameters = dict(rows=1, cols=1, n_samples=10)
+    parameters[name] = -1
+    with pytest.raises(TypeError):
+        Pedestal(**parameters, dtype=dtype)
+
+
+@pytest.mark.parametrize("pedestal_type", [Pedestal_d, Pedestal_f, Pedestal_i16])
+@pytest.mark.parametrize(
+    "args", [(-1, 1), (1, -1), (-1, 1, 10), (1, -1, 10), (1, 1, -1)]
+)
+def test_pedestal_constructor_rejects_negative_parameters(pedestal_type, args):
+    with pytest.raises(TypeError):
+        pedestal_type(*args)
+
+
+@pytest.mark.parametrize("dtype", [np.float64, np.float32, np.int16])
+def test_pedestal_std_stays_finite_after_settling(dtype):
+    pedestal = Pedestal(1, 1, n_samples=10, dtype=dtype)
+    pedestal.push(np.array([[16382]], dtype=np.uint16))
+    frame = np.array([[16383]], dtype=np.uint16)
+    for _ in range(201):
+        pedestal.push(frame)
+
+    noise = pedestal.std().item()
+    assert np.isfinite(noise)
+    assert 0 <= noise < 1e-3
+
+
+@pytest.mark.parametrize(
+    ("pedestal_type", "expected_dtype"),
+    [(Pedestal_d, np.float64), (Pedestal_f, np.float32), (Pedestal_i16, np.int16)],
+)
+def test_double_precision_moments(pedestal_type, expected_dtype):
+    pedestal = pedestal_type(1, 1, 2)
+    for value in [30000, 30003, 30002]:
+        pedestal.push(np.array([[value]], dtype=np.uint16))
+
+    assert pedestal.mean().dtype == expected_dtype
+    assert pedestal.std().dtype == expected_dtype
+    np.testing.assert_array_equal(
+        pedestal.mean(), np.array([[30001.75]], dtype=expected_dtype)
+    )
+    np.testing.assert_allclose(
+        pedestal.std(),
+        np.array([[np.sqrt(1.1875)]], dtype=expected_dtype),
+        rtol=1e-6,
+    )
+
+
+@pytest.mark.parametrize("dtype", [np.float64, np.float32, np.int16])
+@pytest.mark.parametrize("method", ["push", "push_with_threshold"])
+@pytest.mark.parametrize("shape", [(2, 2), (3, 2)])
+def test_pedestal_rejects_mismatched_frame_shapes(dtype, method, shape):
+    pedestal = Pedestal(2, 3, dtype=dtype)
+    initial = np.full((2, 3), 7, dtype=np.uint16)
+    pedestal.push(initial)
+    frame = np.zeros(shape, dtype=np.uint16)
+    args = (
+        (np.full((2, 3), 10, dtype=dtype),)
+        if method == "push_with_threshold"
+        else ()
+    )
+
+    with pytest.raises(
+        RuntimeError, match="Frame shape does not match pedestal shape"
+    ):
+        getattr(pedestal, method)(frame, *args)
+
+    np.testing.assert_array_equal(pedestal.mean(), initial)
 
 
 @pytest.mark.parametrize(
@@ -38,3 +138,67 @@ def test_pedestal_exposes_mean_as_read_only_buffer():
     np.testing.assert_array_equal(mean, pedestal.view())
     assert np.shares_memory(mean, pedestal.view())
     assert not mean.flags.writeable
+
+
+@pytest.mark.parametrize("dtype", [np.float64, np.float32, np.int16])
+@pytest.mark.parametrize("method", ["push", "push_with_threshold"])
+def test_pedestal_push_accepts_contiguous_frames(dtype, method):
+    pedestal = Pedestal(2, 3, dtype=dtype)
+    frame = np.arange(6, dtype=np.uint16).reshape(2, 3)
+    args = (
+        (np.full((2, 3), 10, dtype=dtype),)
+        if method == "push_with_threshold"
+        else ()
+    )
+
+    getattr(pedestal, method)(frame, *args)
+
+    np.testing.assert_array_equal(pedestal.mean(), frame)
+
+
+@pytest.mark.parametrize("dtype", [np.float64, np.float32, np.int16])
+@pytest.mark.parametrize("input_name", ["push_frame", "threshold_frame", "threshold"])
+@pytest.mark.parametrize("layout", ["transpose", "slice", "reverse"])
+def test_pedestal_rejects_noncontiguous_inputs(dtype, input_name, layout):
+    pedestal = Pedestal(2, 3, dtype=dtype)
+    frame = np.zeros((2, 3), dtype=np.uint16)
+    threshold = np.full((2, 3), 10, dtype=dtype)
+    input_dtype = dtype if input_name == "threshold" else np.uint16
+    if layout == "transpose":
+        invalid = np.ones((3, 2), dtype=input_dtype).T
+    elif layout == "slice":
+        invalid = np.ones((2, 6), dtype=input_dtype)[:, ::2]
+    else:
+        invalid = np.ones((4, 6), dtype=input_dtype)[:2, :3][:, ::-1]
+    assert not invalid.flags.c_contiguous
+
+    with pytest.raises(ValueError, match="C-contiguous"):
+        if input_name == "push_frame":
+            pedestal.push(invalid)
+        elif input_name == "threshold_frame":
+            pedestal.push_with_threshold(invalid, threshold)
+        else:
+            pedestal.push_with_threshold(frame, invalid)
+
+    np.testing.assert_array_equal(pedestal.mean(), np.zeros((2, 3)))
+
+
+@pytest.mark.parametrize("dtype", [np.float64, np.float32, np.int16])
+@pytest.mark.parametrize("input_name", ["push_frame", "threshold_frame", "threshold"])
+@pytest.mark.parametrize("shape", [(), (6,), (2, 3, 2)])
+def test_pedestal_rejects_inputs_with_wrong_ndim(dtype, input_name, shape):
+    pedestal = Pedestal(2, 3, dtype=dtype)
+    frame = np.ones((2, 3), dtype=np.uint16)
+    threshold = np.full((2, 3), 10, dtype=dtype)
+    input_dtype = dtype if input_name == "threshold" else np.uint16
+    invalid = np.ones(shape, dtype=input_dtype)
+
+    with pytest.raises(ValueError):
+        if input_name == "push_frame":
+            pedestal.push(invalid)
+        elif input_name == "threshold_frame":
+            pedestal.push_with_threshold(invalid, threshold)
+        else:
+            pedestal.push_with_threshold(frame, invalid)
+
+    np.testing.assert_array_equal(pedestal.mean(), np.zeros((2, 3)))

python/tests/test_ROI.py

@@ -0,0 +1,22 @@
+# SPDX-License-Identifier: MPL-2.0
+import pytest 
+
+import numpy as np
+
+from aare import ROI 
+
+
+def test_roi_slice():
+    """test the ROI slice method"""
+    roi = ROI(0, 4, 2, 4)
+
+    array = np.array([[1, 2, 3, 4],
+                      [5, 6, 7, 8],
+                      [9, 10, 11, 12],
+                      [13, 14, 15, 16]])
+    
+    sliced_array = array[roi.slice()]
+
+    assert sliced_array.shape == (2, 4)
+    assert np.array_equal(sliced_array, np.array([[9, 10, 11, 12],
+                                                   [13, 14, 15, 16]]))

python/tests/test_RawFile.py

@@ -1,8 +1,299 @@
 # SPDX-License-Identifier: MPL-2.0
 import pytest
-from aare import RawFile
+import json
+from aare import File, RawFile, RawSubFile, DetectorType, ROI, UDPPortPosition
 import numpy as np
 
+
+@pytest.fixture
+def small_raw_file(tmp_path):
+    master_path = tmp_path / "run_master_0.json"
+    master_path.write_text(json.dumps({
+        "Version": 7.2,
+        "Detector Type": "Jungfrau",
+        "Timing Mode": "auto",
+        "Geometry": {"x": 1, "y": 1},
+        "Image Size in bytes": 12,
+        "Pixels": {"x": 3, "y": 2},
+        "Max Frames Per File": 1,
+        "Total Frames": 2,
+        "Frames in File": 2,
+        "Frame Padding": 1,
+        "Frame Discard Policy": "nodiscard",
+    }))
+    for index in range(2):
+        data_path = tmp_path / f"run_d0_f{index}_0.raw"
+        # A detector header followed by six uint16 pixels.
+        data_path.write_bytes(bytes(112) + np.arange(6, dtype=np.uint16).tobytes())
+    return master_path
+
+
+@pytest.mark.parametrize("reader_type", [RawFile, File])
+@pytest.mark.parametrize("legacy_master", [False, True])
+@pytest.mark.parametrize("padding, policy, supported", [
+    (0, "nodiscard", False),
+    (0, "discard", False),
+    (0, "discardpartial", True),
+    (1, "nodiscard", True),
+    (1, "discard", True),
+    (1, "discardpartial", True),
+])
+def test_raw_frame_policy(small_raw_file, reader_type, legacy_master,
+                          padding, policy, supported):
+    master_path = small_raw_file
+    if legacy_master:
+        master_path = master_path.with_suffix(".raw")
+        master_path.write_text(
+            "Version : 6.4\n"
+            "Detector Type : Jungfrau\n"
+            "Timing Mode : auto\n"
+            "Geometry : [1, 1]\n"
+            "Image Size : 12\n"
+            "Pixels : [3, 2]\n"
+            "Dynamic Range : 16\n"
+            "Max Frames Per File : 1\n"
+            "Total Frames : 2\n"
+            "Frames in File : 2\n"
+            f"Frame Padding : {padding}\n"
+            f"Frame Discard Policy : {policy}\n"
+        )
+    else:
+        metadata = json.loads(master_path.read_text())
+        metadata["Frame Padding"] = padding
+        metadata["Frame Discard Policy"] = policy
+        master_path.write_text(json.dumps(metadata))
+
+    if supported:
+        reader = reader_type(master_path)
+        assert reader.total_frames == 2
+        frame = reader.read_frame()
+        if reader_type is RawFile:
+            _, frame = frame
+        np.testing.assert_array_equal(frame, np.arange(6).reshape(2, 3))
+    else:
+        message = "requires frame padding or discardpartial"
+        with pytest.raises(RuntimeError, match=message) as error:
+            reader_type(master_path)
+        assert str(master_path) in str(error.value)
+
+        for index in range(2):
+            (master_path.parent / f"run_d0_f{index}_0.raw").unlink()
+        with pytest.raises(RuntimeError, match=message) as error:
+            reader_type(master_path)
+        assert str(master_path) in str(error.value)
+
+
+@pytest.mark.parametrize("method, args, kwargs", [
+    ("read_frame", (), {}),
+    ("read_n", (1,), {}),
+    ("read_roi", (0,), {}),
+    ("read_rois", (), {}),
+    ("read_n_with_roi", (1,), {"roi_index": 0}),
+    ("frame_number", (2,), {}),
+])
+def test_raw_read_error_context(small_raw_file, method, args, kwargs):
+    reader = RawFile(small_raw_file)
+    reader.seek(2)
+    with pytest.raises(RuntimeError) as error:
+        getattr(reader, method)(*args, **kwargs)
+    assert "frame index 2" in str(error.value)
+    assert str(small_raw_file) in str(error.value)
+    assert ".raw" not in str(error.value)
+    if method == "read_frame":
+        assert "file contains 2 frames (indices are zero-based)" in str(error.value)
+
+
+@pytest.mark.parametrize("frame_index", [2, 17])
+def test_generic_raw_read_error_context(small_raw_file, frame_index):
+    reader = File(small_raw_file)
+    with pytest.raises(RuntimeError) as error:
+        reader.read_frame(frame_index)
+    assert f"frame index {frame_index}" in str(error.value)
+    assert str(small_raw_file) in str(error.value)
+
+
+@pytest.mark.parametrize("raw_subfile", [False, True])
+def test_truncated_raw_read_error_context(small_raw_file, raw_subfile):
+    data_path = small_raw_file.parent / "run_d0_f1_0.raw"
+    if raw_subfile:
+        reader = RawSubFile(
+            small_raw_file.parent / "run_d0_f0_0.raw",
+            DetectorType.Jungfrau, 2, 3, 16,
+        )
+    else:
+        reader = RawFile(small_raw_file)
+    data_path.write_bytes(bytes(112))
+    reader.seek(1)
+    with pytest.raises(RuntimeError) as error:
+        reader.read_frame()
+    assert "frame index 1" in str(error.value)
+    assert str(data_path) in str(error.value)
+    assert "End of file" in str(error.value)
+    assert "master" not in str(error.value)
+    assert "RawSubFile.cpp" not in str(error.value)
+
+
+@pytest.mark.parametrize("method, args", [
+    ("read_frame", ()),
+    ("read_n", (1,)),
+    ("read_roi", (0,)),
+    ("frame_number", (1,)),
+])
+def test_short_raw_subfile_sets_read_bounds(small_raw_file, method, args):
+    (small_raw_file.parent / "run_d0_f1_0.raw").unlink()
+    reader = RawFile(small_raw_file)
+    reader.seek(1)
+    with pytest.raises(RuntimeError) as error:
+        getattr(reader, method)(*args)
+    assert reader.total_frames == 1
+    assert "frame index 1" in str(error.value)
+    assert str(small_raw_file) in str(error.value)
+    assert ".raw" not in str(error.value)
+
+
+@pytest.mark.parametrize("counts, master_count, warn", [
+    ([2], 2, False),
+    ([2, 2], 2, False),
+    ([2, 2], 0, True),
+    ([2, 2], 1, True),
+    ([2, 2], 5, True),
+    ([1, 2], 2, True),
+    ([2, 1], 2, True),
+    ([1, 2], 1, True),
+    ([0, 2], 2, True),
+    ([0, 0], 0, False),
+])
+def test_raw_subfile_frame_counts(small_raw_file, capfd, counts, master_count, warn):
+    metadata = json.loads(small_raw_file.read_text())
+    metadata["Geometry"]["y"] = len(counts)
+    metadata["Frames in File"] = master_count
+    metadata["Total Frames"] = 1000
+    small_raw_file.write_text(json.dumps(metadata))
+    data = (small_raw_file.parent / "run_d0_f0_0.raw").read_bytes()
+    for module, count in enumerate(counts):
+        for index in range(2):
+            path = small_raw_file.parent / f"run_d{module}_f{index}_0.raw"
+            if index < count:
+                path.write_bytes(data)
+            elif index == 0:
+                path.write_bytes(b"")
+            else:
+                path.unlink(missing_ok=True)
+
+    reader = RawFile(small_raw_file)
+    warning = capfd.readouterr().err
+    if warn:
+        assert warning.count("WARNING") == 1
+        assert str(small_raw_file) in warning
+    else:
+        assert "WARNING" not in warning
+
+    assert reader.total_frames == min(counts)
+    assert len(reader) == min(counts)
+    assert reader.master.frames_in_file == master_count
+    if min(counts):
+        _, frames = reader.read()
+        assert frames.shape == (min(counts), 2 * len(counts), 3)
+        reader.seek(0)
+        _, frames = reader.read_n(10)
+        assert len(frames) == min(counts)
+        assert reader.frame_number(min(counts) - 1) == 0
+    with pytest.raises(RuntimeError):
+        reader.frame_number(min(counts))
+    reader.seek(min(counts))
+    with pytest.raises(RuntimeError):
+        reader.read_frame()
+    assert "WARNING" not in capfd.readouterr().err
+
+
+@pytest.mark.parametrize("short_module", [0, 1])
+def test_raw_frame_count_across_rois(small_raw_file, short_module):
+    metadata = json.loads(small_raw_file.read_text())
+    metadata.update({
+        "Version": 8.0,
+        "Geometry": {"x": 1, "y": 2},
+        "Image Size": 12,
+        "Receiver Rois": [
+            {"xmin": 0, "xmax": 2, "ymin": 0, "ymax": 1},
+            {"xmin": 0, "xmax": 2, "ymin": 2, "ymax": 3},
+        ],
+    })
+    small_raw_file.write_text(json.dumps(metadata))
+    for index in range(2):
+        data = (small_raw_file.parent / f"run_d0_f{index}_0.raw").read_bytes()
+        (small_raw_file.parent / f"run_d1_f{index}_0.raw").write_bytes(data)
+    (small_raw_file.parent / f"run_d{short_module}_f1_0.raw").unlink()
+
+    reader = RawFile(small_raw_file)
+    assert reader.total_frames == 1
+    _, rois = reader.read_rois()
+    assert len(rois) == 2
+    assert all(roi.shape == (2, 3) for roi in rois)
+    for roi_index in range(2):
+        reader.seek(0)
+        _, frames = reader.read_n_with_roi(10, roi_index=roi_index)
+        assert frames.shape == (1, 2, 3)
+        with pytest.raises(RuntimeError):
+            reader.read_roi(roi_index)
+
+
+@pytest.mark.parametrize("disabled_port", [0, 1])
+def test_raw_frame_count_with_disabled_udp_port(small_raw_file, disabled_port):
+    metadata = json.loads(small_raw_file.read_text())
+    metadata.update({
+        "Geometry": {"x": 1, "y": 2},
+        "Frames in File": 0,
+        "Number of UDP Interfaces": 2,
+        "UDP Ports Type": ["bottom", "top"],
+        "UDP Ports Disabled": [disabled_port],
+    })
+    small_raw_file.write_text(json.dumps(metadata))
+    data = (small_raw_file.parent / "run_d0_f0_0.raw").read_bytes()
+    for index in range(2):
+        (small_raw_file.parent / f"run_d0_f{index}_0.raw").unlink()
+        path = small_raw_file.parent / f"run_d{1 - disabled_port}_f{index}_0.raw"
+        path.write_bytes(data)
+
+    reader = RawFile(small_raw_file)
+    assert reader.total_frames == 2
+    assert reader.master.frames_in_file == 0
+    assert reader.frame_number(1) == 0
+    _, frames = reader.read()
+    assert frames.shape == (2, 2, 3)
+
+
+def test_raw_frame_count_without_selected_subfiles(small_raw_file):
+    metadata = json.loads(small_raw_file.read_text())
+    metadata.update({"Version": 8.0, "Image Size": 12, "Receiver Rois": []})
+    small_raw_file.write_text(json.dumps(metadata))
+    with pytest.raises(RuntimeError, match="No raw subfiles selected"):
+        RawFile(small_raw_file)
+
+
+@pytest.mark.parametrize("lagging_module", [0, 1])
+def test_raw_synchronization_error_identifies_subfile(
+    small_raw_file, lagging_module
+):
+    metadata = json.loads(small_raw_file.read_text())
+    metadata["Geometry"]["y"] = 2
+    small_raw_file.write_text(json.dumps(metadata))
+    for index in range(2):
+        data = (small_raw_file.parent / f"run_d0_f{index}_0.raw").read_bytes()
+        (small_raw_file.parent / f"run_d1_f{index}_0.raw").write_bytes(data)
+
+    ahead_path = small_raw_file.parent / f"run_d{1 - lagging_module}_f1_0.raw"
+    with ahead_path.open("r+b") as output:
+        output.write(np.array([1], dtype=np.uint64).tobytes())
+
+    reader = RawFile(small_raw_file)
+    reader.seek(1)
+    with pytest.raises(RuntimeError) as error:
+        reader.read_frame()
+    assert "frame index 1" in str(error.value)
+    assert "synchroniz" in str(error.value)
+    lagging_path = small_raw_file.parent / f"run_d{lagging_module}_f1_0.raw"
+    assert str(lagging_path) in str(error.value)
+
 @pytest.mark.withdata
 def test_read_rawfile_with_roi_spanning_over_one_module(test_data_path):
 
@@ -110,3 +401,41 @@ def test_read_rawfile_eiger_and_compare_to_numpy(test_data_path):
         header, image1 = f.read_frame()
 
     assert (image == image1).all()
+
+
+@pytest.mark.withdata 
+def test_read_eiger_udp_port_disabled(test_data_path):
+    with RawFile(test_data_path / "raw/eiger/one_udp_port_disabled_master_0.json") as f:
+        _, frame = f.read_rois() 
+
+        assert(len(frame) == 2)
+        assert frame[0].shape == (256, 512)
+        assert frame[1].shape == (256, 1024)
+        assert f.master.udp_port_types == [UDPPortPosition.LEFT, UDPPortPosition.RIGHT]
+        rois = f.master.rois
+        assert len(rois) == 2
+        assert rois[0] == ROI(512, 1024, 0, 256)
+        assert rois[1] == ROI(0, 1024, 256, 512)
+
+    with RawFile(test_data_path / "raw/eiger/quad_eiger_disabled_bottom_port_master_0.json") as f:
+        _, frame = f.read_frame() 
+
+        assert frame.shape == (256, 512)
+        assert(f.master.disabled_udp_ports == [1])
+        assert f.master.udp_port_types == [UDPPortPosition.TOP, UDPPortPosition.BOTTOM]
+        rois = f.master.rois
+        assert len(rois) == 1
+        assert rois[0] == ROI(0, 512, 256, 512)
+
+    with RawFile(test_data_path / "raw/eiger/2_modules_eiger_disabled_udp_port_master_0.json") as f:
+        _, frame = f.read_rois() 
+
+        assert(len(frame) == 2)
+        assert frame[0].shape == (512, 512)
+        assert frame[1].shape == (512, 512)
+        assert (f.master.disabled_udp_ports == [1, 3, 5, 7])
+        assert f.master.udp_port_types == [UDPPortPosition.LEFT, UDPPortPosition.RIGHT]
+        rois = f.master.rois
+        assert len(rois) == 2
+        assert rois[0] == ROI(0, 512, 0, 512)
+        assert rois[1] == ROI(1024, 1536, 0, 512)

python/tests/test_RawMasterFile.py

@@ -1,7 +1,10 @@
+# SPDX-License-Identifier: MPL-2.0
 
 
+import json
+
 import pytest
-from aare import RawMasterFile, ReadoutMode, DetectorType
+from aare import RawMasterFile, ReadoutMode, DetectorType, FrameDiscardPolicy
 
 
 @pytest.mark.withdata
@@ -11,4 +14,37 @@ def test_read_rawfile_quad_eiger_and_compare_to_numpy(test_data_path):
    
     f = RawMasterFile(file_name)
     assert(f.reading_mode == ReadoutMode.UNKNOWN)
-    assert(f.detector_type == DetectorType.Jungfrau)
+    assert(f.detector_type == DetectorType.Jungfrau)
+
+
+@pytest.mark.parametrize(
+    "policy, expected_policy",
+    [
+        ("nodiscard", FrameDiscardPolicy.NoDiscard),
+        ("discard", FrameDiscardPolicy.Discard),
+        ("discardpartial", FrameDiscardPolicy.DiscardPartial),
+    ],
+)
+def test_raw_master_file_context_manager(tmp_path, policy, expected_policy):
+    file_name = tmp_path / "run_master_0.json"
+    file_name.write_text(json.dumps({
+        "Version": 7.2,
+        "Detector Type": "Jungfrau",
+        "Timing Mode": "auto",
+        "Geometry": {"x": 1, "y": 1},
+        "Image Size in bytes": 12,
+        "Pixels": {"x": 3, "y": 2},
+        "Max Frames Per File": 1,
+        "Total Frames": 2,
+        "Frames in File": 2,
+        "Frame Padding": 1,
+        "Frame Discard Policy": policy,
+    }))
+
+    
+    with RawMasterFile(file_name) as context_file:
+        assert context_file.reading_mode == ReadoutMode.UNKNOWN
+        assert context_file.detector_type == DetectorType.Jungfrau
+        frame_discard_policy = context_file.frame_discard_policy
+        assert isinstance(frame_discard_policy, FrameDiscardPolicy)
+        assert frame_discard_policy == expected_policy

python/tests/test_Transformation.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 import pytest 
 import numpy as np
 from aare import CtbRawFile, transform

python/tests/test_interpolation_simulated.py

@@ -1,3 +1,4 @@
+# SPDX-License-Identifier: MPL-2.0
 
 #test script to test interpolation on simulated data 
 

python/tests/test_make_view.py

@@ -0,0 +1,74 @@
+# SPDX-License-Identifier: MPL-2.0
+"""Probe make_view in np_helper.hpp through the bindings in _aare.testing"""
+import numpy as np
+import pytest
+
+from aare._aare import testing
+
+
+
+
+def test_make_view_roundtrip():
+    """Send an array to C++ and check that we get the same array back"""
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    out = testing.view_roundtrip_2d(arr)
+    np.testing.assert_array_equal(out, arr)
+
+    #Sanity check to make sure we got a copy an not just the same thing back
+    arr[3,2] = 31655
+    assert arr[3,2] != out[3,2]
+
+def test_make_view_const_roundtrip():
+    """Const works if the array is not writeable"""
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    arr.flags.writeable = False
+    out = testing.const_view_roundtrip_2d(arr)
+    np.testing.assert_array_equal(out, arr)
+
+def test_make_view_const_roundtrip_writable():
+    """Const works also if the array is  writeable"""
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    out = testing.const_view_roundtrip_2d(arr)
+    np.testing.assert_array_equal(out, arr)
+
+def test_make_view_fill():
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    testing.view_fill_2d(arr, 7)
+    assert (arr == 7).all()
+
+def test_make_view_fill_rejects_nonwriteable():
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    arr.flags.writeable = False
+    with pytest.raises(ValueError, match="not writeable"):
+        testing.view_fill_2d(arr, 7)
+
+def test_make_view_2d_rejects_1d():
+    arr = np.arange(30, dtype=np.uint16)
+    with pytest.raises(ValueError, match="Expected 2D array, got 1D"):
+        testing.view_roundtrip_2d(arr)
+
+def test_make_view_2d_rejects_3d():
+    arr = np.arange(30, dtype=np.uint16).reshape(3, 5, 2)
+    with pytest.raises(ValueError, match="Expected 2D array, got 3D"):
+        testing.view_roundtrip_2d(arr)
+
+def test_make_view_rejects_noncontiguous():
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    arr = arr[::2, ::3]  # non-contiguous
+    assert not arr.flags.c_contiguous
+    with pytest.raises(ValueError, match="not C-contiguous"):
+        testing.view_roundtrip_2d(arr)
+
+def test_make_view_rejects_noncontiguous_const():
+    arr = np.arange(30, dtype=np.float64).reshape(5, 6)
+    arr = arr[::2, ::3]  # non-contiguous
+    assert not arr.flags.c_contiguous
+    with pytest.raises(ValueError, match="not C-contiguous"):
+        testing.const_view_roundtrip_2d(arr)
+
+def test_make_view_rejects_fortran_order():
+    arr = np.asfortranarray(np.arange(30, dtype=np.float64).reshape(5, 6))
+    assert not arr.flags.c_contiguous
+    with pytest.raises(ValueError, match="not C-contiguous"):
+        testing.view_roundtrip_2d(arr)
+