solved merge conflict

.gitea/workflows/cmake_build.yml

@@ -0,0 +1,58 @@
+name: Build the package using cmake then documentation
+
+on:
+  workflow_dispatch:
+  push:
+    
+
+
+permissions:
+  contents: read
+  pages: write
+  id-token: write
+
+jobs:
+  build:
+    strategy:
+      fail-fast: false
+      matrix:
+        platform: [ubuntu-latest, ] # macos-12, windows-2019]
+        python-version: ["3.12",]
+
+    runs-on: ${{ matrix.platform }}
+
+    # The setup-miniconda action needs this to activate miniconda
+    defaults:
+      run:
+        shell: "bash -l {0}"
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Setup dev env
+        run: |
+            sudo apt-get update
+            sudo apt-get -y install cmake gcc g++
+
+      - name: Get conda
+        uses: conda-incubator/setup-miniconda@v3.0.4
+        with:
+          python-version: ${{ matrix.python-version }}
+          channels: conda-forge
+
+      - name: Prepare
+        run: conda install doxygen sphinx=7.1.2 breathe pybind11 sphinx_rtd_theme furo nlohmann_json zeromq fmt numpy
+
+      - name: Build library 
+        run: |
+            mkdir build 
+            cd build 
+            cmake .. -DAARE_SYSTEM_LIBRARIES=ON -DAARE_DOCS=ON 
+            make -j 2 
+            make docs
+
+
+
+
+
+

CMakeLists.txt

@@ -81,15 +81,29 @@ set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 
 if(AARE_FETCH_LMFIT)
     #TODO! Should we fetch lmfit from the web or inlcude a tar.gz in the repo?
-    set(lmfit_patch git apply ${CMAKE_CURRENT_SOURCE_DIR}/patches/lmfit.patch)
+    set(LMFIT_PATCH_COMMAND git apply ${CMAKE_CURRENT_SOURCE_DIR}/patches/lmfit.patch)
-    FetchContent_Declare(
+
-        lmfit
+    # For cmake < 3.28 we can't supply EXCLUDE_FROM_ALL to FetchContent_Declare
-        GIT_REPOSITORY https://jugit.fz-juelich.de/mlz/lmfit.git
+    # so we need this workaround
-        GIT_TAG        main
+    if (${CMAKE_VERSION} VERSION_LESS "3.28")
-        PATCH_COMMAND ${lmfit_patch}
+        FetchContent_Declare(
-        UPDATE_DISCONNECTED 1
+            lmfit
-        #EXCLUDE_FROM_ALL 1
+            GIT_REPOSITORY https://jugit.fz-juelich.de/mlz/lmfit.git
-    )
+            GIT_TAG        main
+            PATCH_COMMAND ${LMFIT_PATCH_COMMAND}
+            UPDATE_DISCONNECTED 1
+        )
+    else()
+        FetchContent_Declare(
+            lmfit
+            GIT_REPOSITORY https://jugit.fz-juelich.de/mlz/lmfit.git
+            GIT_TAG        main
+            PATCH_COMMAND ${LMFIT_PATCH_COMMAND}
+            UPDATE_DISCONNECTED 1
+            EXCLUDE_FROM_ALL 1
+        )
+    endif()
+    
     #Disable what we don't need from lmfit
     set(BUILD_TESTING OFF CACHE BOOL "")
     set(LMFIT_CPPTEST OFF CACHE BOOL "")
@@ -97,8 +111,15 @@ if(AARE_FETCH_LMFIT)
     set(LMFIT_CPPTEST OFF CACHE BOOL "")
     set(BUILD_SHARED_LIBS OFF CACHE BOOL "")
 
+    if (${CMAKE_VERSION} VERSION_LESS "3.28")
+        if(NOT lmfit_POPULATED)
+            FetchContent_Populate(lmfit)
+            add_subdirectory(${lmfit_SOURCE_DIR} ${lmfit_BINARY_DIR} EXCLUDE_FROM_ALL)
+        endif()
+    else()
+        FetchContent_MakeAvailable(lmfit)
+    endif()
     
-    FetchContent_MakeAvailable(lmfit)
     set_property(TARGET lmfit PROPERTY POSITION_INDEPENDENT_CODE ON)
 else()
     find_package(lmfit REQUIRED)
@@ -111,10 +132,13 @@ if(AARE_FETCH_ZMQ)
     if (${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.30")
         cmake_policy(SET CMP0169 OLD)
     endif()
+    set(ZMQ_PATCH_COMMAND git apply ${CMAKE_CURRENT_SOURCE_DIR}/patches/libzmq_cmake_version.patch)
     FetchContent_Declare(
         libzmq
         GIT_REPOSITORY https://github.com/zeromq/libzmq.git
         GIT_TAG        v4.3.4
+        PATCH_COMMAND ${ZMQ_PATCH_COMMAND}
+        UPDATE_DISCONNECTED 1
     )
     # Disable unwanted options from libzmq
     set(BUILD_TESTS OFF CACHE BOOL "Switch off libzmq test build")
@@ -395,6 +419,7 @@ if(AARE_TESTS)
             ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinder.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterVector.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/Cluster.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFile.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp

conda-recipe/meta.yaml

@@ -1,6 +1,7 @@
 package:
   name: aare
-  version: 2025.2.18 #TODO! how to not duplicate this?
+  version: 2025.4.1 #TODO! how to not duplicate this?
+
 
 
 

include/aare/Cluster.hpp

@@ -8,7 +8,9 @@
 
 #pragma once
 
+#include <algorithm>
 #include <cstdint>
+#include <numeric>
 #include <type_traits>
 
 namespace aare {
@@ -28,6 +30,61 @@ struct Cluster {
     CoordType x;
     CoordType y;
     T data[ClusterSizeX * ClusterSizeY];
+
+    T sum() const {
+        return std::accumulate(data, data + ClusterSizeX * ClusterSizeY, 0);
+    }
+
+    T max_sum_2x2() const {
+
+        constexpr size_t num_2x2_subclusters =
+            (ClusterSizeX - 1) * (ClusterSizeY - 1);
+        std::array<T, num_2x2_subclusters> sum_2x2_subcluster;
+        for (size_t i = 0; i < ClusterSizeY - 1; ++i) {
+            for (size_t j = 0; j < ClusterSizeX - 1; ++j)
+                sum_2x2_subcluster[i * (ClusterSizeX - 1) + j] =
+                    data[i * ClusterSizeX + j] +
+                    data[i * ClusterSizeX + j + 1] +
+                    data[(i + 1) * ClusterSizeX + j] +
+                    data[(i + 1) * ClusterSizeX + j + 1];
+        }
+
+        return *std::max_element(sum_2x2_subcluster.begin(),
+                                 sum_2x2_subcluster.end());
+    }
+};
+
+// Specialization for 2x2 clusters (only one sum exists)
+template <typename T> struct Cluster<T, 2, 2, int16_t> {
+    int16_t x;
+    int16_t y;
+    T data[4];
+
+    T sum() const { return std::accumulate(data, data + 4, 0); }
+
+    T max_sum_2x2() const {
+        return data[0] + data[1] + data[2] +
+               data[3]; // Only one possible 2x2 sum
+    }
+};
+
+// Specialization for 3x3 clusters
+template <typename T> struct Cluster<T, 3, 3, int16_t> {
+    int16_t x;
+    int16_t y;
+    T data[9];
+
+    T sum() const { return std::accumulate(data, data + 9, 0); }
+
+    T max_sum_2x2() const {
+        std::array<T, 4> sum_2x2_subclusters;
+        sum_2x2_subclusters[0] = data[0] + data[1] + data[3] + data[4];
+        sum_2x2_subclusters[1] = data[1] + data[2] + data[4] + data[5];
+        sum_2x2_subclusters[2] = data[3] + data[4] + data[6] + data[7];
+        sum_2x2_subclusters[3] = data[4] + data[5] + data[7] + data[8];
+        return *std::max_element(sum_2x2_subclusters.begin(),
+                                 sum_2x2_subclusters.end());
+    }
 };
 
 // Type Traits for is_cluster_type

include/aare/ClusterFile.hpp

@@ -6,6 +6,7 @@
 #include "aare/defs.hpp"
 #include <filesystem>
 #include <fstream>
+#include <optional>
 
 namespace aare {
 
@@ -28,7 +29,7 @@ uint32_t number_of_clusters
  *
  *       int32_t frame_number
  *       uint32_t number_of_clusters
- *       int16_t x, int16_t y, int32_t data[9] x number_of_clusters
+ *       int16_t x, int16_t y, int32_t data[9] * number_of_clusters
  *       int32_t frame_number
  *       uint32_t number_of_clusters
  *       etc.
@@ -37,9 +38,15 @@ template <typename ClusterType,
           typename Enable = std::enable_if_t<is_cluster_v<ClusterType>, bool>>
 class ClusterFile {
     FILE *fp{};
-    uint32_t m_num_left{};
+    uint32_t m_num_left{};    /*Number of photons left in frame*/
-    size_t m_chunk_size{};
+    size_t m_chunk_size{};    /*Number of clusters to read at a time*/
-    const std::string m_mode;
+    const std::string m_mode; /*Mode to open the file in*/
+    std::optional<ROI> m_roi; /*Region of interest, will be applied if set*/
+    std::optional<NDArray<int32_t, 2>>
+        m_noise_map; /*Noise map to cut photons, will be applied if set*/
+    std::optional<NDArray<double, 2>>
+        m_gain_map; /*Gain map to apply to the clusters, will be applied if
+                       set*/
 
   public:
     /**
@@ -75,21 +82,44 @@ class ClusterFile {
 
     void write_frame(const ClusterVector<ClusterType> &clusters);
 
-    // Need to be migrated to support NDArray and return a ClusterVector
-    // std::vector<Cluster3x3>
-    // read_cluster_with_cut(size_t n_clusters, double *noise_map, int nx, int
-    // ny);
-
     /**
      * @brief Return the chunk size
      */
     size_t chunk_size() const { return m_chunk_size; }
 
+    /**
+     * @brief Set the region of interest to use when reading clusters. If set
+     * only clusters within the ROI will be read.
+     */
+    void set_roi(ROI roi);
+
+    /**
+     * @brief Set the noise map to use when reading clusters. If set clusters
+     * below the noise level will be discarded. Selection criteria one of:
+     * Central pixel above noise, highest 2x2 sum above 2 * noise, total sum
+     * above 3 * noise.
+     */
+    void set_noise_map(const NDView<int32_t, 2> noise_map);
+
+    /**
+     * @brief Set the gain map to use when reading clusters. If set the gain map
+     * will be applied to the clusters that pass ROI and noise_map selection.
+     */
+    void set_gain_map(const NDView<double, 2> gain_map);
+
     /**
      * @brief Close the file. If not closed the file will be closed in the
      * destructor
      */
     void close();
+
+  private:
+    ClusterVector<ClusterType> read_clusters_with_cut(size_t n_clusters);
+    ClusterVector<ClusterType> read_clusters_without_cut(size_t n_clusters);
+    ClusterVector<ClusterType> read_frame_with_cut();
+    ClusterVector<ClusterType> read_frame_without_cut();
+    bool is_selected(ClusterType &cl);
+    ClusterType read_one_cluster();
 };
 
 template <typename ClusterType, typename Enable>
@@ -133,6 +163,20 @@ void ClusterFile<ClusterType, Enable>::close() {
         fp = nullptr;
     }
 }
+template <typename ClusterType, typename Enable>
+void ClusterFile<ClusterType, Enable>::set_roi(ROI roi) {
+    m_roi = roi;
+}
+template <typename ClusterType, typename Enable>
+void ClusterFile<ClusterType, Enable>::set_noise_map(
+    const NDView<int32_t, 2> noise_map) {
+    m_noise_map = NDArray<int32_t, 2>(noise_map);
+}
+template <typename ClusterType, typename Enable>
+void ClusterFile<ClusterType, Enable>::set_gain_map(
+    const NDView<double, 2> gain_map) {
+    m_gain_map = NDArray<double, 2>(gain_map);
+}
 
 // TODO generally supported for all clsuter types
 template <typename ClusterType, typename Enable>
@@ -158,6 +202,19 @@ ClusterFile<ClusterType, Enable>::read_clusters(size_t n_clusters) {
     if (m_mode != "r") {
         throw std::runtime_error("File not opened for reading");
     }
+    if (m_noise_map || m_roi) {
+        return read_clusters_with_cut(n_clusters);
+    } else {
+        return read_clusters_without_cut(n_clusters);
+    }
+}
+
+template <typename ClusterType, typename Enable>
+ClusterVector<ClusterType>
+ClusterFile<ClusterType, Enable>::read_clusters_without_cut(size_t n_clusters) {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
 
     ClusterVector<ClusterType> clusters(n_clusters);
 
@@ -185,6 +242,7 @@ ClusterFile<ClusterType, Enable>::read_clusters(size_t n_clusters) {
     if (nph_read < n_clusters) {
         // keep on reading frames and photons until reaching n_clusters
         while (fread(&iframe, sizeof(iframe), 1, fp)) {
+            clusters.set_frame_number(iframe);
             // read number of clusters in frame
             if (fread(&nph, sizeof(nph), 1, fp)) {
                 if (nph > (n_clusters - nph_read))
@@ -204,93 +262,121 @@ ClusterFile<ClusterType, Enable>::read_clusters(size_t n_clusters) {
     // Resize the vector to the number of clusters.
     // No new allocation, only change bounds.
     clusters.resize(nph_read);
+    if (m_gain_map)
+        clusters.apply_gain_map(m_gain_map->view());
     return clusters;
 }
 
 template <typename ClusterType, typename Enable>
 ClusterVector<ClusterType>
-ClusterFile<ClusterType, Enable>::read_clusters(size_t n_clusters, ROI roi) {
+ClusterFile<ClusterType, Enable>::read_clusters_with_cut(size_t n_clusters) {
-    if (m_mode != "r") {
-        throw std::runtime_error("File not opened for reading");
-    }
-
     ClusterVector<ClusterType> clusters;
     clusters.reserve(n_clusters);
 
-    int32_t iframe = 0; // frame number needs to be 4 bytes!
-    size_t nph_read = 0;
-    uint32_t nn = m_num_left;
-    uint32_t nph = m_num_left; // number of clusters in frame needs to be 4
-
-    // auto buf = reinterpret_cast<Cluster3x3 *>(clusters.data());
-    // auto buf = clusters.data();
-
-    ClusterType tmp; // this would break if the cluster size changes
-
     // if there are photons left from previous frame read them first
-    if (nph) {
+    if (m_num_left) {
-        if (nph > n_clusters) {
+        while (m_num_left && clusters.size() < n_clusters) {
-            // if we have more photons left in the frame then photons to read we
+            ClusterType c = read_one_cluster();
-            // read directly the requested number
+            if (is_selected(c)) {
-            nn = n_clusters;
+                clusters.push_back(c);
-        } else {
-            nn = nph;
-        }
-        // Read one cluster, in the ROI push back
-        //  nph_read += fread((buf + nph_read*clusters.item_size()),
-        //                    clusters.item_size(), nn, fp);
-        for (size_t i = 0; i < nn; i++) {
-            fread(&tmp, sizeof(tmp), 1, fp);
-            if (tmp.x >= roi.xmin && tmp.x <= roi.xmax && tmp.y >= roi.ymin &&
-                tmp.y <= roi.ymax) {
-                // clusters.push_back(tmp.x, tmp.y,
-                // reinterpret_cast<std::byte *>(tmp.data));
-                clusters.push_back(tmp);
-                nph_read++;
             }
         }
-
-        m_num_left = nph - nn; // write back the number of photons left
     }
 
-    if (nph_read < n_clusters) {
+    // we did not have enough clusters left in the previous frame
-        // keep on reading frames and photons until reaching n_clusters
+    // keep on reading frames until reaching n_clusters
-        while (fread(&iframe, sizeof(iframe), 1, fp)) {
+    if (clusters.size() < n_clusters) {
-            // read number of clusters in frame
+        // sanity check
-            if (fread(&nph, sizeof(nph), 1, fp)) {
+        if (m_num_left) {
-                if (nph > (n_clusters - nph_read))
+            throw std::runtime_error(
-                    nn = n_clusters - nph_read;
+                LOCATION + "Entered second loop with clusters left\n");
-                else
+        }
-                    nn = nph;
 
-                // nph_read += fread((buf + nph_read*clusters.item_size()),
+        int32_t frame_number = 0; // frame number needs to be 4 bytes!
-                //                   clusters.item_size(), nn, fp);
+        while (fread(&frame_number, sizeof(frame_number), 1, fp)) {
-                for (size_t i = 0; i < nn; i++) {
+            if (fread(&m_num_left, sizeof(m_num_left), 1, fp)) {
-                    fread(&tmp, sizeof(tmp), 1, fp);
+                clusters.set_frame_number(
-                    if (tmp.x >= roi.xmin && tmp.x <= roi.xmax &&
+                    frame_number); // cluster vector will hold the last frame
-                        tmp.y >= roi.ymin && tmp.y <= roi.ymax) {
+                                   // number
-                        // clusters.push_back(
+                while (m_num_left && clusters.size() < n_clusters) {
-                        // tmp.x, tmp.y,
+                    ClusterType c = read_one_cluster();
-                        // reinterpret_cast<std::byte *>(tmp.data));
+                    if (is_selected(c)) {
-                        clusters.push_back(tmp);
+                        clusters.push_back(c);
-                        nph_read++;
                     }
                 }
-                m_num_left = nph - nn;
             }
-            if (nph_read >= n_clusters)
+
+            // we have enough clusters, break out of the outer while loop
+            if (clusters.size() >= n_clusters)
                 break;
         }
     }
+    if (m_gain_map)
+        clusters.apply_gain_map(m_gain_map->view());
 
-    // Resize the vector to the number of clusters.
-    // No new allocation, only change bounds.
-    clusters.resize(nph_read);
     return clusters;
 }
 
+template <typename ClusterType, typename Enable>
+ClusterType ClusterFile<ClusterType, Enable>::read_one_cluster() {
+    ClusterType c;
+    auto rc = fread(&c, sizeof(c), 1, fp);
+    if (rc != 1) {
+        throw std::runtime_error(LOCATION + "Could not read cluster");
+    }
+    --m_num_left;
+    return c;
+}
+
 template <typename ClusterType, typename Enable>
 ClusterVector<ClusterType> ClusterFile<ClusterType, Enable>::read_frame() {
+    if (m_mode != "r") {
+        throw std::runtime_error(LOCATION + "File not opened for reading");
+    }
+    if (m_noise_map || m_roi) {
+        return read_frame_with_cut();
+    } else {
+        return read_frame_without_cut();
+    }
+}
+
+template <typename ClusterType, typename Enable>
+ClusterVector<ClusterType>
+ClusterFile<ClusterType, Enable>::read_frame_without_cut() {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+    if (m_num_left) {
+        throw std::runtime_error(
+            "There are still photons left in the last frame");
+    }
+    int32_t frame_number;
+    if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) {
+        throw std::runtime_error(LOCATION + "Could not read frame number");
+    }
+
+    int32_t n_clusters; // Saved as 32bit integer in the cluster file
+    if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
+        throw std::runtime_error(LOCATION +
+                                 "Could not read number of clusters");
+    }
+
+    ClusterVector<ClusterType> clusters(n_clusters);
+    clusters.set_frame_number(frame_number);
+
+    if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
+        static_cast<size_t>(n_clusters)) {
+        throw std::runtime_error(LOCATION + "Could not read clusters");
+    }
+    clusters.resize(n_clusters);
+    if (m_gain_map)
+        clusters.apply_gain_map(m_gain_map->view());
+    return clusters;
+}
+
+template <typename ClusterType, typename Enable>
+ClusterVector<ClusterType>
+ClusterFile<ClusterType, Enable>::read_frame_with_cut() {
     if (m_mode != "r") {
         throw std::runtime_error("File not opened for reading");
     }
@@ -303,20 +389,45 @@ ClusterVector<ClusterType> ClusterFile<ClusterType, Enable>::read_frame() {
         throw std::runtime_error("Could not read frame number");
     }
 
-    int32_t n_clusters; // Saved as 32bit integer in the cluster file
+    if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) {
-    if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
         throw std::runtime_error("Could not read number of clusters");
     }
-    // std::vector<Cluster3x3> clusters(n_clusters);
-    ClusterVector<ClusterType> clusters(n_clusters);
-    clusters.set_frame_number(frame_number);
 
-    if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
+    ClusterVector<ClusterType> clusters;
-        static_cast<size_t>(n_clusters)) {
+    clusters.reserve(m_num_left);
-        throw std::runtime_error("Could not read clusters");
+    clusters.set_frame_number(frame_number);
+    while (m_num_left) {
+        ClusterType c = read_one_cluster();
+        if (is_selected(c)) {
+            clusters.push_back(c);
+        }
     }
-    clusters.resize(n_clusters);
+    if (m_gain_map)
+        clusters.apply_gain_map(m_gain_map->view());
     return clusters;
 }
 
+template <typename ClusterType, typename Enable>
+bool ClusterFile<ClusterType, Enable>::is_selected(ClusterType &cl) {
+    // Should fail fast
+    if (m_roi) {
+        if (!(m_roi->contains(cl.x, cl.y))) {
+            return false;
+        }
+    }
+    if (m_noise_map) {
+        int32_t sum_1x1 = cl.data[4];       // central pixel
+        int32_t sum_2x2 = cl.max_sum_2x2(); // highest sum of 2x2 subclusters
+        int32_t sum_3x3 = cl.sum();         // sum of all pixels
+
+        auto noise =
+            (*m_noise_map)(cl.y, cl.x); // TODO! check if this is correct
+        if (sum_1x1 <= noise || sum_2x2 <= 2 * noise || sum_3x3 <= 3 * noise) {
+            return false;
+        }
+    }
+    // we passed all checks
+    return true;
+}
+
 } // namespace aare

include/aare/ClusterVector.hpp

@@ -9,6 +9,9 @@
 
 #include <fmt/core.h>
 
+#include "aare/Cluster.hpp"
+#include "aare/NDView.hpp"
+
 namespace aare {
 
 template <typename ClusterType,
@@ -235,6 +238,10 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
         return *reinterpret_cast<const ClusterType *>(element_ptr(i));
     }
 
+    template <typename V> const V &at(size_t i) const {
+        return *reinterpret_cast<const V *>(element_ptr(i));
+    }
+
     const std::string_view fmt_base() const {
         // TODO! how do we match on coord_t?
         return m_fmt_base;
@@ -265,6 +272,28 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
         m_size = new_size;
     }
 
+    void apply_gain_map(const NDView<double> gain_map){
+        //in principle we need to know the size of the image for this lookup
+        //TODO! check orientations
+        std::array<int64_t, 9> xcorr = {-1, 0, 1, -1, 0, 1, -1, 0, 1};
+        std::array<int64_t, 9> ycorr = {-1, -1, -1, 0, 0, 0, 1, 1, 1};
+        for (size_t i=0; i<m_size; i++){
+            auto& cl = at<Cluster3x3>(i);
+
+            if (cl.x > 0 && cl.y > 0 && cl.x < gain_map.shape(1)-1 && cl.y < gain_map.shape(0)-1){
+                for (size_t j=0; j<9; j++){
+                    size_t x = cl.x + xcorr[j];
+                    size_t y = cl.y + ycorr[j];
+                    cl.data[j] = static_cast<T>(cl.data[j] * gain_map(y, x));
+                }
+            }else{
+                memset(cl.data, 0, 9*sizeof(T)); //clear edge clusters
+            }
+
+            
+        }
+    }
+
   private:
     void allocate_buffer(size_t new_capacity) {
         size_t num_bytes = item_size() * new_capacity;

include/aare/defs.hpp

@@ -1,11 +1,9 @@
 #pragma once
 
 #include "aare/Dtype.hpp"
-// #include "aare/utils/logger.hpp"
 
 #include <array>
 #include <stdexcept>
-
 #include <cassert>
 #include <cstdint>
 #include <cstring>
@@ -43,6 +41,7 @@ inline constexpr size_t bits_per_byte = 8;
 void assert_failed(const std::string &msg);
 
 
+
 class DynamicCluster {
   public:
     int cluster_sizeX;
@@ -215,6 +214,9 @@ struct ROI{
   
     int64_t height() const { return ymax - ymin; }
     int64_t width() const { return xmax - xmin; }
+    bool contains(int64_t x, int64_t y) const {
+        return x >= xmin && x < xmax && y >= ymin && y < ymax;
+    }
   };
 
 

patches/libzmq_cmake_version.patch

@@ -0,0 +1,18 @@
+diff --git a/CMakeLists.txt b/CMakeLists.txt
+index dd3d8eb9..c0187747 100644
+--- a/CMakeLists.txt
++++ b/CMakeLists.txt
+@@ -1,11 +1,8 @@
+ # CMake build script for ZeroMQ
+ project(ZeroMQ)
+ 
+-if(${CMAKE_SYSTEM_NAME} STREQUAL Darwin)
+-  cmake_minimum_required(VERSION 3.0.2)
+-else()
+-  cmake_minimum_required(VERSION 2.8.12)
+-endif()
++cmake_minimum_required(VERSION 3.15)
++message(STATUS "Patched cmake version")
+ 
+ include(CheckIncludeFiles)
+ include(CheckCCompilerFlag)

pyproject.toml

@@ -4,7 +4,8 @@ build-backend = "scikit_build_core.build"
 
 [project]
 name = "aare"
-version = "2025.2.18"
+version = "2025.4.1"
+
 
 
 [tool.scikit-build]

python/aare/__init__.py

@@ -12,7 +12,7 @@ from ._aare import ROI
 from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i
 
 from ._aare import fit_gaus, fit_pol1
-
+from ._aare import Interpolator
 from .CtbRawFile import CtbRawFile
 from .RawFile import RawFile
 from .ScanParameters import ScanParameters

python/examples/play.py

@@ -12,6 +12,7 @@ import math
 import time
 
 
+
 def gaussian_2d(mx, my, sigma = 1, res=100, grid_size = 2):
     """
     Generate a 2D gaussian as position mx, my, with sigma=sigma. 
@@ -54,6 +55,7 @@ p = Interpolator(eta, etabinsx[0:-1], etabinsy[0:-1], ebins[0:-1])
 
 
 
+
 #Generate the hit
 
 

python/src/cluster_file.hpp

@@ -41,16 +41,17 @@ void define_cluster_file_io_bindings(py::module &m) {
                 auto v = new ClusterVector<int32_t>(self.read_frame());
                 return v;
              })
+        .def("set_roi", &ClusterFile::set_roi)
+        .def("set_noise_map", [](ClusterFile &self, py::array_t<int32_t> noise_map) {
+            auto view = make_view_2d(noise_map);
+            self.set_noise_map(view);
+        })
+        .def("set_gain_map", [](ClusterFile &self, py::array_t<double> gain_map) {
+            auto view = make_view_2d(gain_map);
+            self.set_gain_map(view);
+        })
+        .def("close", &ClusterFile::close)
         .def("write_frame", &ClusterFile::write_frame)
-        // .def("read_cluster_with_cut",
-        //      [](ClusterFile &self, size_t n_clusters,
-        //         py::array_t<double> noise_map, int nx, int ny) {
-        //          auto view = make_view_2d(noise_map);
-        //          auto *vec =
-        //              new std::vector<Cluster3x3>(self.read_cluster_with_cut(
-        //                  n_clusters, view.data(), nx, ny));
-        //          return return_vector(vec);
-        //      })
         .def("__enter__", [](ClusterFile &self) { return &self; })
         .def("__exit__",
              [](ClusterFile &self,

src/ClusterFile.cpp

@@ -34,9 +34,25 @@ ClusterFile<ClusterType>::ClusterFile(const std::filesystem::path &fname,
     }
 }
 
+<<<<<<< HEAD
 template <typename ClusterType> ClusterFile<ClusterType>::~ClusterFile() {
     close();
 }
+=======
+void ClusterFile::set_roi(ROI roi){
+    m_roi = roi;
+}
+
+void ClusterFile::set_noise_map(const NDView<int32_t, 2> noise_map){
+    m_noise_map = NDArray<int32_t, 2>(noise_map);
+}
+
+void ClusterFile::set_gain_map(const NDView<double, 2> gain_map){
+    m_gain_map = NDArray<double, 2>(gain_map);
+}
+
+ClusterFile::~ClusterFile() { close(); }
+>>>>>>> developer
 
 template <typename ClusterType> void ClusterFile<ClusterType>::close() {
     if (fp) {
@@ -56,16 +72,43 @@ void ClusterFile<ClusterType>::write_frame(
         !(clusters.cluster_size_y() == 3)) {
         throw std::runtime_error("Only 3x3 clusters are supported");
     }
+    //First write the frame number - 4 bytes
     int32_t frame_number = clusters.frame_number();
-    fwrite(&frame_number, sizeof(frame_number), 1, fp);
+    if(fwrite(&frame_number, sizeof(frame_number), 1, fp)!=1){
+        throw std::runtime_error(LOCATION + "Could not write frame number");
+    }
+
+    //Then write the number of clusters - 4 bytes
     uint32_t n_clusters = clusters.size();
-    fwrite(&n_clusters, sizeof(n_clusters), 1, fp);
+    if(fwrite(&n_clusters, sizeof(n_clusters), 1, fp)!=1){
-    fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp);
+        throw std::runtime_error(LOCATION + "Could not write number of clusters");
+    }
+
+    //Now write the clusters in the frame
+    if(fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp)!=clusters.size()){
+        throw std::runtime_error(LOCATION + "Could not write clusters");
+    }
 }
 
+<<<<<<< HEAD
 template <typename ClusterType>
 ClusterVector<ClusterType>
 ClusterFile<ClusterType>::read_clusters(size_t n_clusters) {
+=======
+
+ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters){
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+    if (m_noise_map || m_roi){
+        return read_clusters_with_cut(n_clusters);
+    }else{
+        return read_clusters_without_cut(n_clusters);
+    }
+}
+
+ClusterVector<int32_t> ClusterFile::read_clusters_without_cut(size_t n_clusters) {
+>>>>>>> developer
     if (m_mode != "r") {
         throw std::runtime_error("File not opened for reading");
     }
@@ -96,6 +139,7 @@ ClusterFile<ClusterType>::read_clusters(size_t n_clusters) {
     if (nph_read < n_clusters) {
         // keep on reading frames and photons until reaching n_clusters
         while (fread(&iframe, sizeof(iframe), 1, fp)) {
+            clusters.set_frame_number(iframe);
             // read number of clusters in frame
             if (fread(&nph, sizeof(nph), 1, fp)) {
                 if (nph > (n_clusters - nph_read))
@@ -115,9 +159,12 @@ ClusterFile<ClusterType>::read_clusters(size_t n_clusters) {
     // Resize the vector to the number of clusters.
     // No new allocation, only change bounds.
     clusters.resize(nph_read);
+    if(m_gain_map)
+        clusters.apply_gain_map(m_gain_map->view());
     return clusters;
 }
 
+<<<<<<< HEAD
 template <typename ClusterType>
 ClusterVector<ClusterType>
 ClusterFile<ClusterType>::read_clusters(size_t n_clusters, ROI roi) {
@@ -200,6 +247,109 @@ ClusterFile<ClusterType>::read_clusters(size_t n_clusters, ROI roi) {
 
 template <typename ClusterType>
 ClusterVector<ClusterType> ClusterFile<ClusterType>::read_frame() {
+=======
+
+
+ClusterVector<int32_t> ClusterFile::read_clusters_with_cut(size_t n_clusters) {
+    ClusterVector<int32_t> clusters(3,3);
+    clusters.reserve(n_clusters);
+
+    // if there are photons left from previous frame read them first
+    if (m_num_left) {
+        while(m_num_left && clusters.size() < n_clusters){
+            Cluster3x3 c = read_one_cluster();
+            if(is_selected(c)){
+                clusters.push_back(c.x, c.y, reinterpret_cast<std::byte*>(c.data));
+            }
+        }
+    }
+
+    // we did not have enough clusters left in the previous frame
+    // keep on reading frames until reaching n_clusters
+    if (clusters.size() < n_clusters) {
+        // sanity check
+        if (m_num_left) {
+            throw std::runtime_error(LOCATION + "Entered second loop with clusters left\n");
+        }
+        
+        int32_t frame_number = 0; // frame number needs to be 4 bytes!
+        while (fread(&frame_number, sizeof(frame_number), 1, fp)) {
+            if (fread(&m_num_left, sizeof(m_num_left), 1, fp)) {
+                clusters.set_frame_number(frame_number); //cluster vector will hold the last frame number
+                while(m_num_left && clusters.size() < n_clusters){
+                    Cluster3x3 c = read_one_cluster();
+                    if(is_selected(c)){
+                        clusters.push_back(c.x, c.y, reinterpret_cast<std::byte*>(c.data));
+                    }
+                }
+            }
+
+            // we have enough clusters, break out of the outer while loop
+            if (clusters.size() >= n_clusters)
+                break;
+        }
+
+    }
+    if(m_gain_map)
+        clusters.apply_gain_map(m_gain_map->view());
+
+    return clusters;
+}
+
+Cluster3x3 ClusterFile::read_one_cluster(){
+    Cluster3x3 c;
+    auto rc = fread(&c, sizeof(c), 1, fp);
+    if (rc != 1) {
+        throw std::runtime_error(LOCATION + "Could not read cluster");
+    }
+    --m_num_left;
+    return c;
+}
+
+ClusterVector<int32_t> ClusterFile::read_frame(){
+    if (m_mode != "r") {
+        throw std::runtime_error(LOCATION + "File not opened for reading");
+    }
+    if (m_noise_map || m_roi){
+        return read_frame_with_cut();
+    }else{
+        return read_frame_without_cut();
+    }
+}
+
+ClusterVector<int32_t> ClusterFile::read_frame_without_cut() {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+    if (m_num_left) {
+        throw std::runtime_error(
+            "There are still photons left in the last frame");
+    }
+    int32_t frame_number;
+    if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) {
+        throw std::runtime_error(LOCATION + "Could not read frame number");
+    }
+
+    int32_t n_clusters; // Saved as 32bit integer in the cluster file
+    if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
+        throw std::runtime_error(LOCATION + "Could not read number of clusters");
+    }
+
+    ClusterVector<int32_t> clusters(3, 3, n_clusters);
+    clusters.set_frame_number(frame_number);
+
+    if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
+        static_cast<size_t>(n_clusters)) {
+        throw std::runtime_error(LOCATION + "Could not read clusters");
+    }
+    clusters.resize(n_clusters);
+    if (m_gain_map)
+        clusters.apply_gain_map(m_gain_map->view());
+    return clusters;
+}
+
+ClusterVector<int32_t> ClusterFile::read_frame_with_cut() {
+>>>>>>> developer
     if (m_mode != "r") {
         throw std::runtime_error("File not opened for reading");
     }
@@ -212,22 +362,31 @@ ClusterVector<ClusterType> ClusterFile<ClusterType>::read_frame() {
         throw std::runtime_error("Could not read frame number");
     }
 
-    int32_t n_clusters; // Saved as 32bit integer in the cluster file
+
-    if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
+    if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) {
         throw std::runtime_error("Could not read number of clusters");
     }
+<<<<<<< HEAD
     // std::vector<Cluster3x3> clusters(n_clusters);
     ClusterVector<ClusterType> clusters(n_clusters);
+=======
+    
+    ClusterVector<int32_t> clusters(3, 3);
+    clusters.reserve(m_num_left);
+>>>>>>> developer
     clusters.set_frame_number(frame_number);
-
+    while(m_num_left){
-    if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
+        Cluster3x3 c = read_one_cluster();
-        static_cast<size_t>(n_clusters)) {
+        if(is_selected(c)){
-        throw std::runtime_error("Could not read clusters");
+            clusters.push_back(c.x, c.y, reinterpret_cast<std::byte*>(c.data));
+        }
     }
-    clusters.resize(n_clusters);
+    if (m_gain_map)
+        clusters.apply_gain_map(m_gain_map->view());
     return clusters;
 }
 
+<<<<<<< HEAD
 // std::vector<Cluster3x3> ClusterFile::read_cluster_with_cut(size_t n_clusters,
 //                                                            double *noise_map,
 //                                                            int nx, int ny) {
@@ -245,17 +404,23 @@ ClusterVector<ClusterType> ClusterFile<ClusterType>::read_frame() {
 //     //     uint32_t nn = *n_left;
 //     uint32_t nn = m_num_left;
 //     size_t nph_read = 0;
+=======
 
-//     int32_t t2max, tot1;
+>>>>>>> developer
-//     int32_t tot3;
-//     // Cluster *ptr = buf;
-//     Cluster3x3 *ptr = clusters.data();
-//     int good = 1;
-//     double noise;
-//     // read photons left from previous frame
-//     if (noise_map)
-//         printf("Using noise map\n");
 
+bool ClusterFile::is_selected(Cluster3x3 &cl) {
+    //Should fail fast
+    if (m_roi) {
+        if (!(m_roi->contains(cl.x, cl.y))) {
+            return false;
+        }
+    }
+    if (m_noise_map){
+        int32_t sum_1x1 = cl.data[4]; // central pixel
+        int32_t sum_2x2 = cl.sum_2x2(); // highest sum of 2x2 subclusters
+        int32_t sum_3x3 = cl.sum(); // sum of all pixels
+
+<<<<<<< HEAD
 //     if (nph) {
 //         if (nph > n_clusters) {
 //             // if we have more photons left in the frame then photons to
@@ -359,11 +524,22 @@ ClusterVector<ClusterType> ClusterFile<ClusterType>::read_frame() {
 //     clusters.resize(nph_read);
 //     return clusters;
 // }
+=======
+        auto noise = (*m_noise_map)(cl.y, cl.x); //TODO! check if this is correct
+        if (sum_1x1 <= noise || sum_2x2 <= 2 * noise || sum_3x3 <= 3 * noise) {
+            return false;
+        }
+    }
+    //we passed all checks
+    return true;
+}
+>>>>>>> developer
 
 template <typename ClusterType>
 NDArray<double, 2> calculate_eta2(ClusterVector<ClusterType> &clusters) {
     // TOTO! make work with 2x2 clusters
     NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
+<<<<<<< HEAD
 
     for (size_t i = 0; i < clusters.size(); i++) {
         auto e = calculate_eta2<ClusterType>(clusters.at(i));
@@ -371,6 +547,25 @@ NDArray<double, 2> calculate_eta2(ClusterVector<ClusterType> &clusters) {
         eta2(i, 1) = e.y;
     }
 
+=======
+    
+    if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
+        for (size_t i = 0; i < clusters.size(); i++) {
+            auto e = calculate_eta2(clusters.at<Cluster3x3>(i));
+            eta2(i, 0) = e.x;
+            eta2(i, 1) = e.y;
+        }
+    }else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 2){
+        for (size_t i = 0; i < clusters.size(); i++) {
+            auto e = calculate_eta2(clusters.at<Cluster2x2>(i));
+            eta2(i, 0) = e.x;
+            eta2(i, 1) = e.y;
+        }
+    }else{
+        throw std::runtime_error("Only 3x3 and 2x2 clusters are supported");
+    }
+    
+>>>>>>> developer
     return eta2;
 }
 
@@ -476,6 +671,7 @@ template <typename T> Eta2 calculate_eta2(Cluster<T, 3, 3> &cl) {
     return eta;
 }
 
+<<<<<<< HEAD
 template <typename T> Eta2 calculate_eta2(Cluster<T, 2, 2> &cl) {
     Eta2 eta{};
 
@@ -490,107 +686,19 @@ template <typename T> Eta2 calculate_eta2(Cluster<T, 2, 2> &cl) {
 int analyze_cluster(Cluster<int32_t, 3, 3> &cl, int32_t *t2, int32_t *t3,
                     char *quad, double *eta2x, double *eta2y, double *eta3x,
                     double *eta3y) {
+=======
+>>>>>>> developer
 
-    return analyze_data(cl.data, t2, t3, quad, eta2x, eta2y, eta3x, eta3y);
+Eta2 calculate_eta2(Cluster2x2 &cl) {
+    Eta2 eta{};
+    if ((cl.data[0] + cl.data[1]) != 0)
+        eta.x = static_cast<double>(cl.data[1]) / (cl.data[0] + cl.data[1]);
+    if ((cl.data[0] + cl.data[2]) != 0)
+        eta.y = static_cast<double>(cl.data[2]) / (cl.data[0] + cl.data[2]);
+    eta.sum = cl.data[0] + cl.data[1] + cl.data[2]+ cl.data[3];
+    eta.c = cBottomLeft; //TODO! This is not correct, but need to put something
+    return eta;
 }
 
-int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
-                 double *eta2x, double *eta2y, double *eta3x, double *eta3y) {
-
-    int ok = 1;
-
-    int32_t tot2[4];
-    int32_t t2max = 0;
-    char c = 0;
-    int32_t val, tot3;
-
-    tot3 = 0;
-    for (int i = 0; i < 4; i++)
-        tot2[i] = 0;
-
-    for (int ix = 0; ix < 3; ix++) {
-        for (int iy = 0; iy < 3; iy++) {
-            val = data[iy * 3 + ix];
-            //	printf ("%d ",data[iy * 3 + ix]);
-            tot3 += val;
-            if (ix <= 1 && iy <= 1)
-                tot2[cBottomLeft] += val;
-            if (ix >= 1 && iy <= 1)
-                tot2[cBottomRight] += val;
-            if (ix <= 1 && iy >= 1)
-                tot2[cTopLeft] += val;
-            if (ix >= 1 && iy >= 1)
-                tot2[cTopRight] += val;
-        }
-        //	printf ("\n");
-    }
-    // printf ("\n");
-
-    if (t2 || quad) {
-
-        t2max = tot2[0];
-        c = cBottomLeft;
-        for (int i = 1; i < 4; i++) {
-            if (tot2[i] > t2max) {
-                t2max = tot2[i];
-                c = i;
-            }
-        }
-        // printf("*** %d %d %d %d --
-        // %d\n",tot2[0],tot2[1],tot2[2],tot2[3],t2max);
-        if (quad)
-            *quad = c;
-        if (t2)
-            *t2 = t2max;
-    }
-
-    if (t3)
-        *t3 = tot3;
-
-    if (eta2x || eta2y) {
-        if (eta2x)
-            *eta2x = 0;
-        if (eta2y)
-            *eta2y = 0;
-        switch (c) {
-        case cBottomLeft:
-            if (eta2x && (data[3] + data[4]) != 0)
-                *eta2x = static_cast<double>(data[4]) / (data[3] + data[4]);
-            if (eta2y && (data[1] + data[4]) != 0)
-                *eta2y = static_cast<double>(data[4]) / (data[1] + data[4]);
-            break;
-        case cBottomRight:
-            if (eta2x && (data[2] + data[5]) != 0)
-                *eta2x = static_cast<double>(data[5]) / (data[4] + data[5]);
-            if (eta2y && (data[1] + data[4]) != 0)
-                *eta2y = static_cast<double>(data[4]) / (data[1] + data[4]);
-            break;
-        case cTopLeft:
-            if (eta2x && (data[7] + data[4]) != 0)
-                *eta2x = static_cast<double>(data[4]) / (data[3] + data[4]);
-            if (eta2y && (data[7] + data[4]) != 0)
-                *eta2y = static_cast<double>(data[7]) / (data[7] + data[4]);
-            break;
-        case cTopRight:
-            if (eta2x && t2max != 0)
-                *eta2x = static_cast<double>(data[5]) / (data[5] + data[4]);
-            if (eta2y && t2max != 0)
-                *eta2y = static_cast<double>(data[7]) / (data[7] + data[4]);
-            break;
-        default:;
-        }
-    }
-
-    if (eta3x || eta3y) {
-        if (eta3x && (data[3] + data[4] + data[5]) != 0)
-            *eta3x = static_cast<double>(-data[3] + data[3 + 2]) /
-                     (data[3] + data[4] + data[5]);
-        if (eta3y && (data[1] + data[4] + data[7]) != 0)
-            *eta3y = static_cast<double>(-data[1] + data[2 * 3 + 1]) /
-                     (data[1] + data[4] + data[7]);
-    }
-
-    return ok;
-}
 
 } // namespace aare

src/ClusterFile.test.cpp

@@ -0,0 +1,80 @@
+#include "aare/ClusterFile.hpp"
+#include "test_config.hpp"
+
+
+#include "aare/defs.hpp"
+#include <catch2/catch_test_macros.hpp>
+#include <filesystem>
+
+
+
+
+using aare::ClusterFile;
+
+TEST_CASE("Read one frame from a a cluster file", "[.integration]") {
+    //We know that the frame has 97 clusters
+    auto fpath = test_data_path() / "clusters" / "single_frame_97_clustrers.clust";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    ClusterFile f(fpath);
+    auto clusters = f.read_frame();
+    REQUIRE(clusters.size() == 97);
+    REQUIRE(clusters.frame_number() == 135);
+}
+
+TEST_CASE("Read one frame using ROI", "[.integration]") {
+    //We know that the frame has 97 clusters
+    auto fpath = test_data_path() / "clusters" / "single_frame_97_clustrers.clust";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    ClusterFile f(fpath);
+    aare::ROI roi;
+    roi.xmin = 0;
+    roi.xmax = 50;
+    roi.ymin = 200;
+    roi.ymax = 249;
+    f.set_roi(roi);
+    auto clusters = f.read_frame();
+    REQUIRE(clusters.size() == 49);
+    REQUIRE(clusters.frame_number() == 135);
+
+    //Check that all clusters are within the ROI
+    for (size_t i = 0; i < clusters.size(); i++) {
+        auto c = clusters.at<aare::Cluster3x3>(i);
+        REQUIRE(c.x >= roi.xmin);
+        REQUIRE(c.x <= roi.xmax);
+        REQUIRE(c.y >= roi.ymin);
+        REQUIRE(c.y <= roi.ymax);
+    }
+
+}
+
+
+TEST_CASE("Read clusters from single frame file", "[.integration]") {
+
+    auto fpath = test_data_path() / "clusters" / "single_frame_97_clustrers.clust";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    SECTION("Read fewer clusters than available") {
+        ClusterFile f(fpath);
+        auto clusters = f.read_clusters(50);
+        REQUIRE(clusters.size() == 50);
+         REQUIRE(clusters.frame_number() == 135);
+    }
+    SECTION("Read more clusters than available") {
+        ClusterFile f(fpath);
+        // 100 is the maximum number of clusters read
+        auto clusters = f.read_clusters(100);
+        REQUIRE(clusters.size() == 97);
+        REQUIRE(clusters.frame_number() == 135);
+    }
+    SECTION("Read all clusters") {
+        ClusterFile f(fpath);
+        auto clusters = f.read_clusters(97);
+        REQUIRE(clusters.size() == 97);
+        REQUIRE(clusters.frame_number() == 135);
+    }
+
+
+    
+}

src/Interpolator.cpp

@@ -14,41 +14,41 @@ Interpolator::Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
             "The shape of the etacube does not match the shape of the bins");
     }
 
-    // Cumulative sum in the x direction, can maybe be combined with a copy?
+    // Cumulative sum in the x direction
-    for (ssize_t k = 0; k < m_ietax.shape(2); k++) {
+    for (ssize_t i = 1; i < m_ietax.shape(0); i++) {
         for (ssize_t j = 0; j < m_ietax.shape(1); j++) {
-            for (ssize_t i = 1; i < m_ietax.shape(0); i++) {
+            for (ssize_t k = 0; k < m_ietax.shape(2); k++) {
                 m_ietax(i, j, k) += m_ietax(i - 1, j, k);
             }
         }
     }
 
     // Normalize by the highest row, if norm less than 1 don't do anything
-    for (ssize_t k = 0; k < m_ietax.shape(2); k++) {
+    for (ssize_t i = 0; i < m_ietax.shape(0); i++) {
         for (ssize_t j = 0; j < m_ietax.shape(1); j++) {
-            auto val = m_ietax(m_ietax.shape(0) - 1, j, k);
+            for (ssize_t k = 0; k < m_ietax.shape(2); k++) {
-            double norm = val < 1 ? 1 : val;
+                auto val = m_ietax(m_ietax.shape(0) - 1, j, k);
-            for (ssize_t i = 0; i < m_ietax.shape(0); i++) {
+                double norm = val < 1 ? 1 : val;
                 m_ietax(i, j, k) /= norm;
             }
         }
     }
 
     // Cumulative sum in the y direction
-    for (ssize_t k = 0; k < m_ietay.shape(2); k++) {
+    for (ssize_t i = 0; i < m_ietay.shape(0); i++) {
-        for (ssize_t i = 0; i < m_ietay.shape(0); i++) {
+        for (ssize_t j = 1; j < m_ietay.shape(1); j++) {
-            for (ssize_t j = 1; j < m_ietay.shape(1); j++) {
+            for (ssize_t k = 0; k < m_ietay.shape(2); k++) {
                 m_ietay(i, j, k) += m_ietay(i, j - 1, k);
             }
         }
     }
 
     // Normalize by the highest column, if norm less than 1 don't do anything
-    for (ssize_t k = 0; k < m_ietay.shape(2); k++) {
+    for (ssize_t i = 0; i < m_ietay.shape(0); i++) {
-        for (ssize_t i = 0; i < m_ietay.shape(0); i++) {
+        for (ssize_t j = 0; j < m_ietay.shape(1); j++) {
-            auto val = m_ietay(i, m_ietay.shape(1) - 1, k);
+            for (ssize_t k = 0; k < m_ietay.shape(2); k++) {
-            double norm = val < 1 ? 1 : val;
+                auto val = m_ietay(i, m_ietay.shape(1) - 1, k);
-            for (ssize_t j = 0; j < m_ietay.shape(1); j++) {
+                double norm = val < 1 ? 1 : val;
                 m_ietay(i, j, k) /= norm;
             }
         }
@@ -110,8 +110,8 @@ Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
                 dY = -1.;
                 break;
             }
-            photon.x += m_ietax(ix, iy, 0) * 2 + dX;
+            photon.x += m_ietax(ix, iy, ie) * 2 + dX;
-            photon.y += m_ietay(ix, iy, 0) * 2 + dY;
+            photon.y += m_ietay(ix, iy, ie) * 2 + dY;
             photons.push_back(photon);
         }
     } else if (clusters.cluster_size_x() == 2 ||
@@ -136,10 +136,10 @@ Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
             auto ix = last_smaller(m_etabinsx, eta.x);
             auto iy = last_smaller(m_etabinsy, eta.y);
 
-            photon.x +=
+            photon.x += m_ietax(ix, iy, ie) *
-                m_ietax(ix, iy, 0) * 2; // eta goes between 0 and 1 but we could
+                        2; // eta goes between 0 and 1 but we could move the hit
-                                        // move the hit anywhere in the 2x2
+                           // anywhere in the 2x2
-            photon.y += m_ietay(ix, iy, 0) * 2;
+            photon.y += m_ietay(ix, iy, ie) * 2;
             photons.push_back(photon);
         }
 

src/algorithm.test.cpp

@@ -1,8 +1,7 @@
 
 
-#include <catch2/catch_test_macros.hpp>
 #include <aare/algorithm.hpp>
-
+#include <catch2/catch_test_macros.hpp>
 
 TEST_CASE("Find the closed index in a 1D array", "[algorithm]") {
     aare::NDArray<double, 1> arr({5});
@@ -17,7 +16,7 @@ TEST_CASE("Find the closed index in a 1D array", "[algorithm]") {
     REQUIRE(aare::nearest_index(arr, -1.0) == 0);
 }
 
-TEST_CASE("Passing integers to nearest_index works"){
+TEST_CASE("Passing integers to nearest_index works", "[algorithm]") {
     aare::NDArray<int, 1> arr({5});
     for (size_t i = 0; i < arr.size(); i++) {
         arr[i] = i;
@@ -30,8 +29,7 @@ TEST_CASE("Passing integers to nearest_index works"){
     REQUIRE(aare::nearest_index(arr, -1) == 0);
 }
 
-
+TEST_CASE("nearest_index works with std::vector", "[algorithm]") {
-TEST_CASE("nearest_index works with std::vector"){
     std::vector<double> vec = {0, 1, 2, 3, 4};
     REQUIRE(aare::nearest_index(vec, 2.123) == 2);
     REQUIRE(aare::nearest_index(vec, 2.66) == 3);
@@ -40,7 +38,7 @@ TEST_CASE("nearest_index works with std::vector"){
     REQUIRE(aare::nearest_index(vec, -10.0) == 0);
 }
 
-TEST_CASE("nearest index works with std::array"){
+TEST_CASE("nearest index works with std::array", "[algorithm]") {
     std::array<double, 5> arr = {0, 1, 2, 3, 4};
     REQUIRE(aare::nearest_index(arr, 2.123) == 2);
     REQUIRE(aare::nearest_index(arr, 2.501) == 3);
@@ -49,8 +47,7 @@ TEST_CASE("nearest index works with std::array"){
     REQUIRE(aare::nearest_index(arr, -10.0) == 0);
 }
 
-
+TEST_CASE("last smaller", "[algorithm]") {
-TEST_CASE("last smaller"){
     aare::NDArray<double, 1> arr({5});
     for (size_t i = 0; i < arr.size(); i++) {
         arr[i] = i;
@@ -60,4 +57,13 @@ TEST_CASE("last smaller"){
     REQUIRE(aare::last_smaller(arr, 0.0) == 0);
     REQUIRE(aare::last_smaller(arr, 2.3) == 2);
     REQUIRE(aare::last_smaller(arr, 253.) == 4);
+}
+
+TEST_CASE("returns last bin strictly smaller", "[algorithm]") {
+    aare::NDArray<double, 1> arr({5});
+    for (size_t i = 0; i < arr.size(); i++) {
+        arr[i] = i;
+    }
+    // arr 0, 1, 2, 3, 4
+    REQUIRE(aare::last_smaller(arr, 2.0) == 2);
 }

tests/test_config.hpp.in

@@ -7,6 +7,6 @@ inline auto test_data_path(){
     if(const char* env_p = std::getenv("AARE_TEST_DATA")){
         return std::filesystem::path(env_p);
     }else{
-        throw std::runtime_error("AARE_TEST_DATA_PATH not set");
+        throw std::runtime_error("Path to test data: $AARE_TEST_DATA not set");
     }
 }