From 2efdd17d81ba8e0e3b38d0df95ce90c5b7b400b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Erik=20Fr=C3=B6jdh?= Date: Wed, 30 Sep 2026 14:54:59 +0200 Subject: [PATCH] check ready before updating pd (#391) - For set_nSigma only call update_pedestal if we have accumulated enough samples otherwise just write it - For update threshold check if pd is ready and if not throw closes #390 --- RELEASE.md | 7 +++++ include/aare/ClusterFinder.hpp | 29 ++++++++++++++++---- include/aare/ClusterFinderMT.hpp | 30 ++++++++++++++++----- python/src/bind_ClusterFinder.hpp | 4 +++ python/src/bind_ClusterFinderMT.hpp | 4 +++ python/tests/test_Cluster.py | 37 ++++++++++++++++++++++++++ src/ClusterFinder.test.cpp | 37 ++++++++++++++++++++++++++ src/ClusterFinderMT.test.cpp | 41 +++++++++++++++++++++++++++++ 8 files changed, 177 insertions(+), 12 deletions(-) diff --git a/RELEASE.md b/RELEASE.md index 6cf4c316..7fa0910e 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -4,6 +4,9 @@ ### New Features: +- Added ``pedestal_ready`` to ``ClusterFinder`` and ``ClusterFinderMT`` in C++ + and Python to check whether enough pedestal frames have been pushed to find + clusters. - Added the Python ``FrameDiscardPolicy`` enum with ``NoDiscard``, ``Discard``, and ``DiscardPartial``, enabling access to ``RawMasterFile.frame_discard_policy``. - Added the Python ``Pedestal`` factory with ``dtype`` selection, matching @@ -90,6 +93,10 @@ ### Bugfixes: +- Setting nSigma on ``ClusterFinder`` and ``ClusterFinderMT`` before the + pedestal is ready no longer raises an error. The new value is used when the + threshold is initialized. Previously ``ClusterFinderMT.set_nSigma()`` also + left the worker threads with different nSigma values. - The Python ``Cluster`` constructors validate that the data array holds exactly one value per pixel and raise ``ValueError`` otherwise. Previously a longer array wrote past the end of the cluster data. diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp index 53ac4f73..ca78ad06 100644 --- a/include/aare/ClusterFinder.hpp +++ b/include/aare/ClusterFinder.hpp @@ -68,13 +68,17 @@ class ClusterFinder { } /** - * @brief Set the noise multiplier used for threshold calculation and - * recompute the threshold. + * @brief Set the noise multiplier used for threshold calculation. + * + * The threshold is recomputed immediately if the pedestal is ready, + * otherwise when the pedestal becomes ready. * @param nSigma New per-pixel noise multiplier. */ void set_nSigma(PEDESTAL_TYPE nSigma) { m_nSigma = nSigma; - update_threshold(); + if (m_pedestal.ready()) { + update_threshold(); + } } /** @brief Return the current noise multiplier used for threshold @@ -110,8 +114,23 @@ class ClusterFinder { /** @brief Clear the pedestal and mark it as not ready. */ void clear_pedestal() { m_pedestal.clear(); } - /** @brief Recompute the threshold as noise multiplied by nSigma. */ - void update_threshold() { m_threshold = m_pedestal.std() * m_nSigma; } + /** + * @brief Return whether enough pedestal frames have been pushed to find + * clusters. + */ + bool pedestal_ready() const { return m_pedestal.ready(); } + + /** + * @brief Recompute the threshold as noise multiplied by nSigma. + * @throws std::runtime_error if the pedestal is not ready. + */ + void update_threshold() { + if (!m_pedestal.ready()) { + throw std::runtime_error( + "Pedestal is not ready, cannot update threshold"); + } + m_threshold = m_pedestal.std() * m_nSigma; + } /** * @brief Move out all accumulated clusters and reset the internal vector. diff --git a/include/aare/ClusterFinderMT.hpp b/include/aare/ClusterFinderMT.hpp index b4f4fd22..f57fec88 100644 --- a/include/aare/ClusterFinderMT.hpp +++ b/include/aare/ClusterFinderMT.hpp @@ -1,5 +1,6 @@ // SPDX-License-Identifier: MPL-2.0 #pragma once +#include #include #include #include @@ -131,6 +132,14 @@ class ClusterFinderMT { } } + void wait_for_input_queues() const { + for (auto &q : m_input_queues) { + while (!q->isEmpty()) { + std::this_thread::sleep_for(m_default_wait); + } + } + } + bool output_queues_are_empty() const { for (auto &q : m_output_queues) { if (!q->isEmpty()) { @@ -392,21 +401,28 @@ class ClusterFinderMT { // } /** - * @brief Set the nSigma value for all the cluster finders. + * @brief Set the nSigma value for all the cluster finders. Waits for the + * queued frames to be processed first. * @param nSigma number of sigma above the pedestal to consider a photon * during cluster finding. */ void set_nSigma(const PEDESTAL_TYPE nSigma) { - // Wait for all queues to be empty before changing the sigma - for (auto &q : m_input_queues) { - while (!q->isEmpty()) { - std::this_thread::sleep_for(m_default_wait); - } - } + wait_for_input_queues(); for (auto &cf : m_cluster_finders) { cf->set_nSigma(nSigma); } } + + /** + * @brief Return whether all cluster finders have received enough pedestal + * frames to find clusters. Waits for the queued frames to be processed + * first. + */ + bool pedestal_ready() { + wait_for_input_queues(); + return std::all_of(m_cluster_finders.begin(), m_cluster_finders.end(), + [](const auto &cf) { return cf->pedestal_ready(); }); + } }; } // namespace aare \ No newline at end of file diff --git a/python/src/bind_ClusterFinder.hpp b/python/src/bind_ClusterFinder.hpp index 661b9d75..034b1b88 100644 --- a/python/src/bind_ClusterFinder.hpp +++ b/python/src/bind_ClusterFinder.hpp @@ -55,6 +55,10 @@ void define_ClusterFinder(py::module &m, const std::string &typestr) { &ClusterFinder::clear_pedestal) .def("update_threshold", &ClusterFinder::update_threshold) + .def_property_readonly( + "pedestal_ready", + &ClusterFinder::pedestal_ready, + R"(True once enough pedestal frames have been pushed to find clusters.)") .def_property_readonly( "pedestal", [](ClusterFinder &self) { diff --git a/python/src/bind_ClusterFinderMT.hpp b/python/src/bind_ClusterFinderMT.hpp index 2c166a1a..fdd4b580 100644 --- a/python/src/bind_ClusterFinderMT.hpp +++ b/python/src/bind_ClusterFinderMT.hpp @@ -62,6 +62,10 @@ void define_ClusterFinderMT(py::module &m, const std::string &typestr) { &ClusterFinderMT::clear_pedestal) .def("update_threshold", &ClusterFinderMT::update_threshold) + .def( + "pedestal_ready", + &ClusterFinderMT::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::sync) .def("stop", &ClusterFinderMT::stop) .def("start", &ClusterFinderMT::start) diff --git a/python/tests/test_Cluster.py b/python/tests/test_Cluster.py index 3978ab0c..28b90e9d 100644 --- a/python/tests/test_Cluster.py +++ b/python/tests/test_Cluster.py @@ -135,6 +135,43 @@ def test_cluster_finder(): 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)) diff --git a/src/ClusterFinder.test.cpp b/src/ClusterFinder.test.cpp index fd7eeb37..5ce214fd 100644 --- a/src/ClusterFinder.test.cpp +++ b/src/ClusterFinder.test.cpp @@ -2,6 +2,7 @@ #include "aare/ClusterFinder.hpp" #include "aare/Pedestal.hpp" #include +#include #include #include #include @@ -45,6 +46,42 @@ TEST_CASE("Construct a cluster finder") { // REQUIRE(clusterFinder.get_nSigma() == 1); } +TEST_CASE("nSigma set before the pedestal is ready is used for the threshold") { + const Shape<2> image_size{10, 10}; + const size_t min_pedestal_samples = 4; + + // Pedestal of 100 with a noise of 2 and a 30 ADU hit: a cluster at the + // default nSigma of 5, below threshold at nSigma 50. + NDArray low(image_size, 98); + NDArray high(image_size, 102); + NDArray hit(image_size, 100); + hit(5, 5) = 130; + + const double nSigma = GENERATE(5.0, 50.0); + const size_t expected_clusters = nSigma == 5.0 ? 1 : 0; + + ClusterFinder> cf(image_size, 5.0, 100, + min_pedestal_samples); + REQUIRE_FALSE(cf.pedestal_ready()); + REQUIRE_NOTHROW(cf.set_nSigma(nSigma)); + REQUIRE(cf.get_nSigma() == nSigma); + + for (size_t i = 0; i < min_pedestal_samples / 2; ++i) { + cf.push_pedestal_frame(low.view()); + REQUIRE_FALSE(cf.pedestal_ready()); + cf.push_pedestal_frame(high.view()); + } + REQUIRE(cf.pedestal_ready()); + + cf.find_clusters(hit.view()); + REQUIRE(cf.steal_clusters().size() == expected_clusters); +} + +TEST_CASE("update_threshold throws until the pedestal is ready") { + ClusterFinder> cf({10, 10}); + REQUIRE_THROWS_AS(cf.update_threshold(), std::runtime_error); +} + // TEST_CASE("test cluster finder") { // aare::Pedestal pedestal(10, 10, 5); // NDArray frame({10, 10}); diff --git a/src/ClusterFinderMT.test.cpp b/src/ClusterFinderMT.test.cpp index ca630de4..6973d0df 100644 --- a/src/ClusterFinderMT.test.cpp +++ b/src/ClusterFinderMT.test.cpp @@ -144,6 +144,47 @@ TEST_CASE("frame buffers are recycled when pushing more frames than the pool " CHECK(cf.m_input_queues_are_empty() == true); } +TEST_CASE("nSigma set before the pedestal is ready reaches all threads") { + using ClusterType = Cluster; + + const size_t n_threads = 2; + const size_t min_pedestal_samples = 4; + const Shape<2> image_size{10, 10}; + + ClusterFinderMTWrapper cf(image_size, 5, 200, n_threads, 16, + min_pedestal_samples); + + // Pedestal of 100 with a noise of 2 and a 30 ADU hit: a cluster at the + // default nSigma of 5, below threshold at nSigma 50. + NDArray low(image_size, 98); + NDArray high(image_size, 102); + NDArray hit(image_size, 100); + hit(5, 5) = 130; + + CHECK_FALSE(cf.pedestal_ready()); + REQUIRE_NOTHROW(cf.set_nSigma(50)); + + for (size_t i = 0; i < min_pedestal_samples / 2; ++i) { + cf.push_pedestal_frame(low.view()); + cf.push_pedestal_frame(high.view()); + } + CHECK(cf.pedestal_ready()); + + for (size_t i = 0; i < n_threads; ++i) { + cf.find_clusters(hit.view(), i); + } + cf.stop(); + + ClusterCollector collector(&cf); + collector.stop(); + + auto clusters = collector.steal_clusters(); + REQUIRE(clusters.size() == n_threads); + for (const auto &frame_clusters : clusters) { + CHECK(frame_clusters.size() == 0); + } +} + TEST_CASE("cluster collector accepts finders with a matching cluster type") { using ClusterType = Cluster; using Finder = ClusterFinderMTWrapper;