clear pedestal

slsdetectorgroup · Jan 10, 2025 · 7ce0200 · 7ce0200
1 parent 7550a2c
commit 7ce0200
Show file tree

Hide file tree

Showing 6 changed files with 46 additions and 21 deletions.
diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp
@@ -33,6 +33,12 @@ typedef enum {
     pTopRight = 8
 } pixel;
 
+struct Eta2 {
+    double x;
+    double y;
+    corner c;
+};
+
 struct ClusterAnalysis {
     uint32_t c;
     int32_t tot;
@@ -74,7 +80,8 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
 int analyze_cluster(Cluster3x3& cl, int32_t *t2, int32_t *t3, char *quad,
                     double *eta2x, double *eta2y, double *eta3x, double *eta3y);
 
+
 NDArray<double, 2> calculate_eta2( ClusterVector<int>& clusters);
-std::array<double,2> calculate_eta2( Cluster3x3& cl);
+Eta2 calculate_eta2( Cluster3x3& cl);
 
 } // namespace aare
diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp
@@ -47,6 +47,7 @@ class ClusterFinder {
 
     NDArray<PEDESTAL_TYPE, 2> pedestal() { return m_pedestal.mean(); }
     NDArray<PEDESTAL_TYPE, 2> noise() { return m_pedestal.std(); }
+    void clear_pedestal() { m_pedestal.clear(); }
 
     /**
      * @brief Move the clusters from the ClusterVector in the ClusterFinder to a

diff --git a/include/aare/ClusterFinderMT.hpp b/include/aare/ClusterFinderMT.hpp
@@ -215,6 +215,15 @@ class ClusterFinderMT {
         m_current_thread++;
     }
 
+    void clear_pedestal() {
+        if (!m_processing_threads_stopped) {
+            throw std::runtime_error("ClusterFinderMT is still running");
+        }
+        for (auto &cf : m_cluster_finders) {
+            cf->clear_pedestal();
+        }
+    }
+
     /**
      * @brief Return the pedestal currently used by the cluster finder
      * @param thread_index index of the thread

diff --git a/include/aare/Pedestal.hpp b/include/aare/Pedestal.hpp
@@ -89,6 +89,7 @@ template <typename SUM_TYPE = double> class Pedestal {
         m_sum = 0;
         m_sum2 = 0;
         m_cur_samples = 0;
+        m_mean = 0;
     }
 
 
@@ -97,6 +98,7 @@ template <typename SUM_TYPE = double> class Pedestal {
         m_sum(row, col) = 0;
         m_sum2(row, col) = 0;
         m_cur_samples(row, col) = 0;
+        m_mean(row, col) = 0;
     }
 
 
@@ -119,7 +121,7 @@ template <typename SUM_TYPE = double> class Pedestal {
 
     /**
      * Push but don't update the cached mean. Speeds up the process
-     * when intitializing the pedestal.
+     * when initializing the pedestal.
      * 
      */
     template <typename T> void push_no_update(NDView<T, 2> frame) {

diff --git a/python/src/cluster.hpp b/python/src/cluster.hpp
@@ -72,6 +72,7 @@ void define_cluster_finder_mt_bindings(py::module &m) {
                 return;
             },
             py::arg(), py::arg("frame_number") = 0)
+        .def("clear_pedestal", &ClusterFinderMT<uint16_t, pd_type>::clear_pedestal)
         .def("sync", &ClusterFinderMT<uint16_t, pd_type>::sync)
         .def("stop", &ClusterFinderMT<uint16_t, pd_type>::stop)
         .def("start", &ClusterFinderMT<uint16_t, pd_type>::start)
@@ -121,6 +122,7 @@ void define_cluster_finder_bindings(py::module &m) {
                  auto view = make_view_2d(frame);
                  self.push_pedestal_frame(view);
              })
+        .def("clear_pedestal", &ClusterFinder<uint16_t, pd_type>::clear_pedestal)
         .def_property_readonly("pedestal",
                                [](ClusterFinder<uint16_t, pd_type> &self) {
                                    auto pd = new NDArray<pd_type, 2>{};

diff --git a/src/ClusterFile.cpp b/src/ClusterFile.cpp
@@ -262,19 +262,19 @@ std::vector<Cluster3x3> ClusterFile::read_cluster_with_cut(size_t n_clusters,
 NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters) {
     NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
     for (size_t i = 0; i < clusters.size(); i++) {
-        // int32_t t2;
-        // auto* ptr = reinterpret_cast<int32_t*> (clusters.element_ptr(i) + 2 *
-        // sizeof(int16_t)); analyze_cluster(clusters.at<Cluster3x3>(i), &t2,
-        // nullptr, nullptr, &eta2(i,0), &eta2(i,1) , nullptr, nullptr);
-        auto [x, y] = calculate_eta2(clusters.at<Cluster3x3>(i));
-        eta2(i, 0) = x;
-        eta2(i, 1) = y;
+        auto e = calculate_eta2(clusters.at<Cluster3x3>(i));
+        eta2(i, 0) = e.x;
+        eta2(i, 1) = e.y;
     }
     return eta2;
 }
 
-std::array<double, 2> calculate_eta2(Cluster3x3 &cl) {
-    std::array<double, 2> eta2{};
+/** 
+ * @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2 struct
+ * containing etay, etax and the corner of the cluster. 
+*/
+Eta2 calculate_eta2(Cluster3x3 &cl) {
+    Eta2 eta{};
 
     std::array<int32_t, 4> tot2;
     tot2[0] = cl.data[0] + cl.data[1] + cl.data[3] + cl.data[4];
@@ -287,39 +287,43 @@ std::array<double, 2> calculate_eta2(Cluster3x3 &cl) {
     switch (c) {
     case cBottomLeft:
         if ((cl.data[3] + cl.data[4]) != 0)
-            eta2[0] =
+            eta.x =
                 static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
         if ((cl.data[1] + cl.data[4]) != 0)
-            eta2[1] =
+            eta.y =
                 static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
+        eta.c = cBottomLeft;
         break;
     case cBottomRight:
         if ((cl.data[2] + cl.data[5]) != 0)
-            eta2[0] =
+            eta.x =
                 static_cast<double>(cl.data[5]) / (cl.data[4] + cl.data[5]);
         if ((cl.data[1] + cl.data[4]) != 0)
-            eta2[1] =
+            eta.y =
                 static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
+        eta.c = cBottomRight;
         break;
     case cTopLeft:
         if ((cl.data[7] + cl.data[4]) != 0)
-            eta2[0] =
+            eta.x =
                 static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
         if ((cl.data[7] + cl.data[4]) != 0)
-            eta2[1] =
+            eta.y =
                 static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
+        eta.c = cTopLeft;
         break;
     case cTopRight:
         if ((cl.data[5] + cl.data[4]) != 0)
-            eta2[0] =
+            eta.x =
                 static_cast<double>(cl.data[5]) / (cl.data[5] + cl.data[4]);
         if ((cl.data[7] + cl.data[4]) != 0)
-            eta2[1] =
+            eta.y =
                 static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
+        eta.c = cTopRight;
         break;
-    // default:;
+    // no default to allow compiler to warn about missing cases
     }
-    return eta2;
+    return eta;
 }
 
 int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad,