Developer (#138)

- Fully functioning variable size cluster finder - Added interpolation - Bit reordering for ADC SAR 05 --------- Co-authored-by: Patrick <patrick.sieberer@psi.ch> Co-authored-by: JulianHeymes <julian.heymes@psi.ch> Co-authored-by: Dhanya Thattil <dhanya.thattil@psi.ch> Co-authored-by: xiangyu.xie <xiangyu.xie@psi.ch>
2025-04-19 13:20:03 +02:00 · 2025-03-20 12:52:04 +01:00 · 2025-03-20 12:52:04 +01:00 · 5d8ad27b21
commit 5d8ad27b21
parent b7a47576a1
30 changed files with 743 additions and 61 deletions
--- a/.clang-tidy
+++ b/.clang-tidy
@ -0,0 +1,42 @@
+
+---
+Checks:          '*, 
+                    -altera-*,
+                    -android-cloexec-fopen,
+                    -cppcoreguidelines-pro-bounds-array-to-pointer-decay, 
+                    -cppcoreguidelines-pro-bounds-pointer-arithmetic, 
+                    -fuchsia*,
+                    -readability-else-after-return,
+                    -readability-avoid-const-params-in-decls,
+                    -readability-identifier-length,
+                    -cppcoreguidelines-pro-bounds-constant-array-index,
+                    -cppcoreguidelines-pro-type-reinterpret-cast,
+                    -llvm-header-guard,
+                    -modernize-use-nodiscard,
+                    -misc-non-private-member-variables-in-classes,
+                    -readability-static-accessed-through-instance,
+                    -readability-braces-around-statements,
+                    -readability-isolate-declaration,
+                    -readability-implicit-bool-conversion,
+                    -readability-identifier-length,
+                    -readability-identifier-naming,
+                    -hicpp-signed-bitwise,
+                    -hicpp-no-array-decay, 
+                    -hicpp-braces-around-statements,
+                    -google-runtime-references, 
+                    -google-readability-todo, 
+                    -google-readability-braces-around-statements,
+                    -modernize-use-trailing-return-type,
+                    -llvmlibc-*'
+
+HeaderFilterRegex: \.hpp
+FormatStyle:     none
+CheckOptions:
+  - { key: readability-identifier-naming.NamespaceCase,       value: lower_case }
+ # - { key: readability-identifier-naming.FunctionCase,        value: lower_case }
+  - { key: readability-identifier-naming.ClassCase,           value: CamelCase  }
+ # - { key: readability-identifier-naming.MethodCase,          value: CamelCase  }
+ # - { key: readability-identifier-naming.StructCase,          value: CamelCase  }
+ # - { key: readability-identifier-naming.VariableCase,        value: lower_case }
+  - { key: readability-identifier-naming.GlobalConstantCase,  value: UPPER_CASE }
+...
--- a/.gitea/workflows/cmake_build.yml
+++ b/.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
+
+
+
+
+
+
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@ -60,6 +60,8 @@ if(AARE_SYSTEM_LIBRARIES)
    set(AARE_FETCH_CATCH OFF CACHE BOOL "Disabled FetchContent for catch2" FORCE)
    set(AARE_FETCH_JSON OFF CACHE BOOL "Disabled FetchContent for nlohmann::json" FORCE)
    set(AARE_FETCH_ZMQ OFF CACHE BOOL "Disabled FetchContent for libzmq" FORCE)
+    # Still fetch lmfit when setting AARE_SYSTEM_LIBRARIES since this is not available
+    # on conda-forge
 endif()

 if(AARE_VERBOSE)
@ -78,6 +80,7 @@ endif()
 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)
    FetchContent_Declare(
        lmfit
@ -343,6 +346,7 @@ set(SourceFiles
    ${CMAKE_CURRENT_SOURCE_DIR}/src/geo_helpers.cpp
    ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.cpp
    ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolator.cpp
    ${CMAKE_CURRENT_SOURCE_DIR}/src/PixelMap.cpp
    ${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.cpp
    ${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.cpp
@ -382,6 +386,7 @@ endif()

 if(AARE_TESTS)
    set(TestSources
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.test.cpp
            ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.test.cpp
            ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp
            ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp
--- a/conda-recipe/meta.yaml
+++ b/conda-recipe/meta.yaml
@ -1,6 +1,7 @@
 package:
  name: aare
-  version: 2025.2.18 #TODO! how to not duplicate this?
+  version: 2025.3.18 #TODO! how to not duplicate this?
+



--- a/include/aare/ClusterFile.hpp
+++ b/include/aare/ClusterFile.hpp
@ -8,11 +8,17 @@

 namespace aare {

+//TODO! Template this? 
 struct Cluster3x3 {
    int16_t x;
    int16_t y;
    int32_t data[9];
 };
+struct Cluster2x2 {
+    int16_t x;
+    int16_t y;
+    int32_t data[4];
+};

 typedef enum {
    cBottomLeft = 0,
@ -37,6 +43,7 @@ struct Eta2 {
    double x;
    double y;
    corner c;
+    int32_t sum;
 };

 struct ClusterAnalysis {
@ -97,6 +104,8 @@ class ClusterFile {
     */
    ClusterVector<int32_t> read_clusters(size_t n_clusters);

+    ClusterVector<int32_t> read_clusters(size_t n_clusters, ROI roi);
+
    /**
     * @brief Read a single frame from the file and return the clusters. The
     * cluster vector will have the frame number set.
@ -131,5 +140,6 @@ int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad,

 NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters);
 Eta2 calculate_eta2(Cluster3x3 &cl);
+Eta2 calculate_eta2(Cluster2x2 &cl);

 } // namespace aare
--- a/include/aare/ClusterVector.hpp
+++ b/include/aare/ClusterVector.hpp
@ -231,6 +231,10 @@ template <typename T, typename CoordType = int16_t> class ClusterVector {
        return *reinterpret_cast<V *>(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;
--- a/include/aare/Interpolator.hpp
+++ b/include/aare/Interpolator.hpp
@ -0,0 +1,29 @@
+#pragma once
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/ClusterFile.hpp" //Cluster_3x3
+namespace aare{
+
+struct Photon{
+    double x;
+    double y;
+    double energy;
+};
+
+class Interpolator{
+    NDArray<double, 3> m_ietax;
+    NDArray<double, 3> m_ietay;
+
+    NDArray<double, 1> m_etabinsx;
+    NDArray<double, 1> m_etabinsy;
+    NDArray<double, 1> m_energy_bins;
+    public:
+        Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins, NDView<double, 1> ybins, NDView<double, 1> ebins);
+        NDArray<double, 3> get_ietax(){return m_ietax;}
+        NDArray<double, 3> get_ietay(){return m_ietay;}
+
+        std::vector<Photon> interpolate(const ClusterVector<int32_t>& clusters);
+};
+
+} // namespace aare
--- a/include/aare/NDArray.hpp
+++ b/include/aare/NDArray.hpp
@ -102,6 +102,9 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
    auto begin() { return data_; }
    auto end() { return data_ + size_; }

+    auto begin() const { return data_; }
+    auto end() const { return data_ + size_; }
+
    using value_type = T;

    NDArray &operator=(NDArray &&other) noexcept; // Move assign
@ -388,12 +391,12 @@ NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const T &value) {
    result *= value;
    return result;
 }
-template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print() {
-    if (shape_[0] < 20 && shape_[1] < 20)
-        Print_all();
-    else
-        Print_some();
-}
+// template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print() {
+//     if (shape_[0] < 20 && shape_[1] < 20)
+//         Print_all();
+//     else
+//         Print_some();
+// }

 template <typename T, int64_t Ndim>
 std::ostream &operator<<(std::ostream &os, const NDArray<T, Ndim> &arr) {
--- a/include/aare/RawSubFile.hpp
+++ b/include/aare/RawSubFile.hpp
@ -64,7 +64,7 @@ class RawSubFile {

    size_t bytes_per_frame() const { return m_bytes_per_frame; }
    size_t pixels_per_frame() const { return m_rows * m_cols; }
-    size_t bytes_per_pixel() const { return m_bitdepth / 8; }
+    size_t bytes_per_pixel() const { return m_bitdepth / bits_per_byte; }

 private:
  template <typename T>
--- a/include/aare/VarClusterFinder.hpp
+++ b/include/aare/VarClusterFinder.hpp
@ -7,7 +7,7 @@

 #include "aare/NDArray.hpp"

-const int MAX_CLUSTER_SIZE = 200;
+const int MAX_CLUSTER_SIZE = 50;
 namespace aare {

 template <typename T> class VarClusterFinder {
--- a/include/aare/algorithm.hpp
+++ b/include/aare/algorithm.hpp
@ -0,0 +1,55 @@
+
+#pragma once
+#include <algorithm>
+#include <array>
+#include <vector>
+#include <aare/NDArray.hpp>
+
+namespace aare {
+/**
+ * @brief Find the index of the last element smaller than val
+ * assume a sorted array
+ */
+template <typename T>
+size_t last_smaller(const T* first, const T* last, T val) {
+    for (auto iter = first+1; iter != last; ++iter) {
+        if (*iter > val) {
+            return std::distance(first, iter-1);
+        }
+    }
+    return std::distance(first, last-1);
+}
+
+template <typename T>
+size_t last_smaller(const NDArray<T, 1>& arr, T val) {
+    return last_smaller(arr.begin(), arr.end(), val);
+}
+
+
+template <typename T>
+size_t nearest_index(const T* first, const T* last, T val) {
+    auto iter = std::min_element(first, last,
+    [val](T a, T b) {
+        return std::abs(a - val) < std::abs(b - val);
+    });
+    return std::distance(first, iter);
+}
+
+template <typename T>
+size_t nearest_index(const NDArray<T, 1>& arr, T val) {
+    return nearest_index(arr.begin(), arr.end(), val);
+}
+
+template <typename T>
+size_t nearest_index(const std::vector<T>& vec, T val) {
+    return nearest_index(vec.data(), vec.data()+vec.size(), val);
+}
+
+template <typename T, size_t N>
+size_t nearest_index(const std::array<T,N>& arr, T val) {
+    return nearest_index(arr.data(), arr.data()+arr.size(), val);
+}
+
+
+
+} // namespace aare
--- a/include/aare/defs.hpp
+++ b/include/aare/defs.hpp
@ -38,6 +38,8 @@

 namespace aare {

+inline constexpr size_t bits_per_byte = 8;
+
 void assert_failed(const std::string &msg);


--- a/pyproject.toml
+++ b/pyproject.toml
@ -4,7 +4,8 @@ build-backend = "scikit_build_core.build"

 [project]
 name = "aare"
-version = "2025.2.18"
+version = "2025.3.18"
+


 [tool.scikit-build]
--- a/python/aare/init.py
+++ b/python/aare/init.py
@ -7,11 +7,12 @@ from ._aare import Pedestal_d, Pedestal_f, ClusterFinder, VarClusterFinder
 from ._aare import DetectorType
 from ._aare import ClusterFile
 from ._aare import hitmap
+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
--- a/python/examples/play.py
+++ b/python/examples/play.py
@ -1,50 +1,79 @@
 import sys
 sys.path.append('/home/l_msdetect/erik/aare/build')

-#Our normal python imports
-from pathlib import Path
-import matplotlib.pyplot as plt
+from aare._aare import ClusterVector_i, Interpolator
+
+import pickle 
 import numpy as np
+import matplotlib.pyplot as plt
 import boost_histogram as bh
+import torch
+import math
 import time


-import aare

-data = np.random.normal(10, 1, 1000)
+def gaussian_2d(mx, my, sigma = 1, res=100, grid_size = 2):
+    """
+    Generate a 2D gaussian as position mx, my, with sigma=sigma. 
+    The gaussian is placed on a 2x2 pixel matrix with resolution 
+    res in one dimesion.
+    """
+    x = torch.linspace(0, pixel_size*grid_size, res)
+    x,y = torch.meshgrid(x,x, indexing="ij")
+    return 1 / (2*math.pi*sigma**2) * \
+      torch.exp(-((x - my)**2 / (2*sigma**2) + (y - mx)**2 / (2*sigma**2)))

-hist = bh.Histogram(bh.axis.Regular(10, 0, 20))
-hist.fill(data)
+scale = 1000 #Scale factor when converting to integer
+pixel_size = 25 #um
+grid = 2
+resolution = 100
+sigma_um = 10
+xa = np.linspace(0,grid*pixel_size,resolution)
+ticks = [0, 25, 50]
+
+hit = np.array((20,20))
+etahist_fname = "/home/l_msdetect/erik/tmp/test_hist.pkl"
+
+local_resolution = 99
+grid_size = 3
+xaxis = np.linspace(0,grid_size*pixel_size, local_resolution)
+t = gaussian_2d(hit[0],hit[1], grid_size = grid_size, sigma = 10, res = local_resolution)
+pixels = t.reshape(grid_size, t.shape[0] // grid_size, grid_size, t.shape[1] // grid_size).sum(axis = 3).sum(axis = 1)
+pixels = pixels.numpy()
+pixels = (pixels*scale).astype(np.int32)
+v = ClusterVector_i(3,3)
+v.push_back(1,1, pixels)
+
+with open(etahist_fname, "rb") as f:
+        hist = pickle.load(f)
+eta = hist.view().copy()
+etabinsx = np.array(hist.axes.edges.T[0].flat)
+etabinsy = np.array(hist.axes.edges.T[1].flat)
+ebins = np.array(hist.axes.edges.T[2].flat)
+p = Interpolator(eta, etabinsx[0:-1], etabinsy[0:-1], ebins[0:-1])


-x = hist.axes[0].centers
-y = hist.values()
-y_err = np.sqrt(y)+1
-res = aare.fit_gaus(x, y, y_err, chi2 = True)


-        
-t_elapsed = time.perf_counter()-t0
-print(f'Histogram filling took: {t_elapsed:.3f}s {total_clusters/t_elapsed/1e6:.3f}M clusters/s')
+#Generate the hit

-histogram_data = hist3d.counts()
-x = hist3d.axes[2].edges[:-1]

-y = histogram_data[100,100,:]
-xx = np.linspace(x[0], x[-1])
-# fig, ax = plt.subplots()
-# ax.step(x, y, where = 'post')

-y_err = np.sqrt(y)
-y_err = np.zeros(y.size)
-y_err += 1

-# par = fit_gaus2(y,x, y_err)
-# ax.plot(xx, gaus(xx,par))
-# print(par)
+tmp = p.interpolate(v)
+print(f'tmp:{tmp}')
+pos = np.array((tmp['x'], tmp['y']))*25

-res = fit_gaus(y,x)
-res2 = fit_gaus(y,x, y_err)
-print(res)
-print(res2)

+print(pixels)
+fig, ax = plt.subplots(figsize = (7,7))
+ax.pcolormesh(xaxis, xaxis, t)
+ax.plot(*pos, 'o')
+ax.set_xticks([0,25,50,75])
+ax.set_yticks([0,25,50,75])
+ax.set_xlim(0,75)
+ax.set_ylim(0,75)
+ax.grid()
+print(f'{hit=}')
+print(f'{pos=}')
--- a/python/src/cluster.hpp
+++ b/python/src/cluster.hpp
@ -20,7 +20,13 @@ template <typename T>
 void define_cluster_vector(py::module &m, const std::string &typestr) {
    auto class_name = fmt::format("ClusterVector_{}", typestr);
    py::class_<ClusterVector<T>>(m, class_name.c_str(), py::buffer_protocol())
-        .def(py::init<int, int>())
+        .def(py::init<int, int>(),
+    py::arg("cluster_size_x") = 3, py::arg("cluster_size_y") = 3)
+        .def("push_back",
+             [](ClusterVector<T> &self, int x, int y, py::array_t<T> data) {
+                //  auto view = make_view_2d(data);
+                 self.push_back(x, y, reinterpret_cast<const std::byte*>(data.data()));
+             })
        .def_property_readonly("size", &ClusterVector<T>::size)
        .def("item_size", &ClusterVector<T>::item_size)
        .def_property_readonly("fmt",
@ -38,6 +44,8 @@ void define_cluster_vector(py::module &m, const std::string &typestr) {
            auto *vec = new std::vector<T>(self.sum_2x2());
            return return_vector(vec);
        })
+        .def_property_readonly("cluster_size_x", &ClusterVector<T>::cluster_size_x)
+        .def_property_readonly("cluster_size_y", &ClusterVector<T>::cluster_size_y)
        .def_property_readonly("capacity", &ClusterVector<T>::capacity)
        .def_property("frame_number", &ClusterVector<T>::frame_number,
                      &ClusterVector<T>::set_frame_number)
--- a/python/src/cluster_file.hpp
+++ b/python/src/cluster_file.hpp
@ -31,6 +31,11 @@ void define_cluster_file_io_bindings(py::module &m) {
                auto v = new ClusterVector<int32_t>(self.read_clusters(n_clusters));
                return v;
             },py::return_value_policy::take_ownership)
+        .def("read_clusters",
+                [](ClusterFile &self, size_t n_clusters, ROI roi) {
+                   auto v = new ClusterVector<int32_t>(self.read_clusters(n_clusters, roi));
+                   return v;
+                },py::return_value_policy::take_ownership)
        .def("read_frame",
             [](ClusterFile &self) {
                auto v = new ClusterVector<int32_t>(self.read_frame());
--- a/python/src/ctb_raw_file.hpp
+++ b/python/src/ctb_raw_file.hpp
@ -32,7 +32,7 @@ m.def("adc_sar_05_decode64to16", [](py::array_t<uint8_t> input) {
    }

    //Create a 2D output array with the same shape as the input
-    std::vector<ssize_t> shape{input.shape(0), input.shape(1)/8};
+    std::vector<ssize_t> shape{input.shape(0), input.shape(1)/static_cast<int64_t>(bits_per_byte)};
    py::array_t<uint16_t> output(shape);

    //Create a view of the input and output arrays
@ -53,7 +53,7 @@ m.def("adc_sar_04_decode64to16", [](py::array_t<uint8_t> input) {
    }

    //Create a 2D output array with the same shape as the input
-    std::vector<ssize_t> shape{input.shape(0), input.shape(1)/8};
+    std::vector<ssize_t> shape{input.shape(0), input.shape(1)/static_cast<int64_t>(bits_per_byte)};
    py::array_t<uint16_t> output(shape);

    //Create a view of the input and output arrays
--- a/python/src/file.hpp
+++ b/python/src/file.hpp
@ -195,6 +195,8 @@ void define_file_io_bindings(py::module &m) {

    py::class_<ROI>(m, "ROI")
        .def(py::init<>())
+        .def(py::init<int64_t, int64_t, int64_t, int64_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)
--- a/python/src/interpolation.hpp
+++ b/python/src/interpolation.hpp
@ -0,0 +1,58 @@
+#include "aare/Interpolator.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include "np_helper.hpp"
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+void define_interpolation_bindings(py::module &m) {
+
+    PYBIND11_NUMPY_DTYPE(aare::Photon, x,y,energy);
+
+    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> xbins,
+                         py::array_t<double> ybins, py::array_t<double> ebins) {
+            return Interpolator(make_view_3d(etacube), make_view_1d(xbins),
+                                make_view_1d(ybins), make_view_1d(ebins));
+        }))
+        .def("get_ietax", [](Interpolator& self){
+            auto*ptr = new NDArray<double,3>{};
+            *ptr = self.get_ietax();
+            return return_image_data(ptr);
+        })
+        .def("get_ietay", [](Interpolator& self){
+            auto*ptr = new NDArray<double,3>{};
+            *ptr = self.get_ietay();
+            return return_image_data(ptr);
+        })
+        .def("interpolate", [](Interpolator& self, const ClusterVector<int32_t>& clusters){
+            auto photons = self.interpolate(clusters);
+            auto* ptr = new std::vector<Photon>{photons};
+            return return_vector(ptr);
+        });
+
+    // TODO! Evaluate without converting to double
+    m.def(
+        "hej",
+        []() {
+            // auto boost_histogram = py::module_::import("boost_histogram");
+            // py::object axis =
+            // boost_histogram.attr("axis").attr("Regular")(10, 0.0, 10.0);
+            // py::object histogram = boost_histogram.attr("Histogram")(axis);
+            // return histogram;
+            // return h;
+        },
+        R"(
+        Evaluate a 1D Gaussian function for all points in x using parameters par.
+
+        Parameters
+        ----------
+        x : array_like
+            The points at which to evaluate the Gaussian function.
+        par : array_like
+            The parameters of the Gaussian function. The first element is the amplitude, the second element is the mean, and the third element is the standard deviation.
+        )");
+}
--- a/python/src/module.cpp
+++ b/python/src/module.cpp
@ -9,6 +9,7 @@
 #include "cluster.hpp"
 #include "cluster_file.hpp"
 #include "fit.hpp"
+#include "interpolation.hpp"

 //Pybind stuff
 #include <pybind11/pybind11.h>
@ -31,5 +32,6 @@ PYBIND11_MODULE(_aare, m) {
    define_cluster_collector_bindings(m);
    define_cluster_file_sink_bindings(m);
    define_fit_bindings(m);
+    define_interpolation_bindings(m);

 }
--- a/python/src/np_helper.hpp
+++ b/python/src/np_helper.hpp
@ -40,25 +40,25 @@ template <typename T> py::array return_vector(std::vector<T> *vec) {
 }

 // todo rewrite generic
-template <class T, int Flags> auto get_shape_3d(py::array_t<T, Flags> arr) {
+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)};
 }

-template <class T, int Flags> auto make_view_3d(py::array_t<T, Flags> arr) {
+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));
 }

-template <class T, int Flags> auto get_shape_2d(py::array_t<T, Flags> 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(py::array_t<T, Flags> arr) {
+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) {
+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));
 }
-template <class T, int Flags> auto make_view_1d(py::array_t<T, Flags> 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));
 }
--- a/python/src/var_cluster.hpp
+++ b/python/src/var_cluster.hpp
@ -19,15 +19,24 @@ using namespace::aare;

 void define_var_cluster_finder_bindings(py::module &m) {
    PYBIND11_NUMPY_DTYPE(VarClusterFinder<double>::Hit, size, row, col,
-                         reserved, energy, max);
+                         reserved, energy, max, rows, cols, enes);

    py::class_<VarClusterFinder<double>>(m, "VarClusterFinder")
        .def(py::init<Shape<2>, double>())
        .def("labeled",
             [](VarClusterFinder<double> &self) {
-                 auto ptr = new NDArray<int, 2>(self.labeled());
+                 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_peripheralThresholdFactor",
+             &VarClusterFinder<double>::set_peripheralThresholdFactor)
        .def("find_clusters",
             [](VarClusterFinder<double> &self,
                py::array_t<double, py::array::c_style | py::array::forcecast>
@ -35,6 +44,30 @@ void define_var_cluster_finder_bindings(py::module &m) {
                 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("hits",
+             [](VarClusterFinder<double> &self) {
+                 auto ptr = new std::vector<VarClusterFinder<double>::Hit>(
+                     self.steal_hits());
+                 return return_vector(ptr);
+             })
+        .def("clear_hits",
+             [](VarClusterFinder<double> &self) {
+                 self.clear_hits();
+             })
        .def("steal_hits",
             [](VarClusterFinder<double> &self) {
                 auto ptr = new std::vector<VarClusterFinder<double>::Hit>(
--- a/src/ClusterFile.cpp
+++ b/src/ClusterFile.cpp
@ -108,6 +108,79 @@ ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters) {
    return clusters;
 }

+ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters, ROI roi) {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+    
+    ClusterVector<int32_t> clusters(3,3);
+    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();
+
+    Cluster3x3 tmp; //this would break if the cluster size changes
+
+    // if there are photons left from previous frame read them first
+    if (nph) {
+        if (nph > n_clusters) {
+            // if we have more photons left in the frame then photons to read we
+            // read directly the requested number
+            nn = n_clusters;
+        } 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));
+                nph_read++;
+            }
+        }
+
+        m_num_left = nph - nn; // write back the number of photons left
+    }
+
+    if (nph_read < n_clusters) {
+        // keep on reading frames and photons until reaching n_clusters
+        while (fread(&iframe, sizeof(iframe), 1, fp)) {
+            // read number of clusters in frame
+            if (fread(&nph, sizeof(nph), 1, fp)) {
+                if (nph > (n_clusters - nph_read))
+                    nn = n_clusters - nph_read;
+                else
+                    nn = nph;
+
+                // 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));
+                        nph_read++;
+                    }
+                }
+                m_num_left = nph - nn;
+            }
+            if (nph_read >= n_clusters)
+                break;
+        }
+    }
+
+    // Resize the vector to the number of clusters.
+    // No new allocation, only change bounds.
+    clusters.resize(nph_read);
+    return clusters;
+}
+
 ClusterVector<int32_t> ClusterFile::read_frame() {
    if (m_mode != "r") {
        throw std::runtime_error("File not opened for reading");
@ -268,11 +341,23 @@ ClusterVector<int32_t> ClusterFile::read_frame() {
 NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters) {
    //TOTO! make work with 2x2 clusters
    NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
-    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;
+    
+    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");
    }
+    
    return eta2;
 }

@ -290,7 +375,7 @@ Eta2 calculate_eta2(Cluster3x3 &cl) {
    tot2[3] = cl.data[4] + cl.data[5] + cl.data[7] + cl.data[8];

    auto c = std::max_element(tot2.begin(), tot2.end()) - tot2.begin();
-
+    eta.sum = tot2[c];
    switch (c) {
    case cBottomLeft:
        if ((cl.data[3] + cl.data[4]) != 0)
@ -333,6 +418,20 @@ Eta2 calculate_eta2(Cluster3x3 &cl) {
    return eta;
 }

+
+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_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad,
                    double *eta2x, double *eta2y, double *eta3x,
                    double *eta3y) {
--- a/src/Dtype.cpp
+++ b/src/Dtype.cpp
@ -70,7 +70,7 @@ uint8_t Dtype::bitdepth() const {
 /**
 * @brief Get the number of bytes of the data type
 */
-size_t Dtype::bytes() const { return bitdepth() / 8; }
+size_t Dtype::bytes() const { return bitdepth() / bits_per_byte; }

 /**
 * @brief Construct a DType object from a TypeIndex
--- a/src/File.cpp
+++ b/src/File.cpp
@ -73,7 +73,7 @@ size_t File::tell() const { return file_impl->tell(); }
 size_t File::rows() const { return file_impl->rows(); }
 size_t File::cols() const { return file_impl->cols(); }
 size_t File::bitdepth() const { return file_impl->bitdepth(); }
-size_t File::bytes_per_pixel() const { return file_impl->bitdepth() / 8; }
+size_t File::bytes_per_pixel() const { return file_impl->bitdepth() / bits_per_byte; }

 DetectorType File::detector_type() const { return file_impl->detector_type(); }

--- a/src/Interpolator.cpp
+++ b/src/Interpolator.cpp
@ -0,0 +1,144 @@
+#include "aare/Interpolator.hpp"
+#include "aare/algorithm.hpp"
+
+namespace aare {
+
+Interpolator::Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
+                           NDView<double, 1> ybins, NDView<double, 1> ebins)
+    : m_ietax(etacube), m_ietay(etacube), m_etabinsx(xbins), m_etabinsy(ybins), m_energy_bins(ebins) {
+    if (etacube.shape(0) != xbins.size() || etacube.shape(1) != ybins.size() ||
+        etacube.shape(2) != ebins.size()) {
+        throw std::invalid_argument(
+            "The shape of the etacube does not match the shape of the bins");
+    }
+
+    // Cumulative sum in the x direction
+    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 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 i = 0; i < m_ietax.shape(0); i++) {
+        for (ssize_t j = 0; j < m_ietax.shape(1); j++) {
+            for (ssize_t k = 0; k < m_ietax.shape(2); k++) {
+                auto val = m_ietax(m_ietax.shape(0) - 1, j, k);
+                double norm = val < 1 ? 1 : val;
+                m_ietax(i, j, k) /= norm;
+            }
+        }
+    }
+
+    // Cumulative sum in the y direction
+    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 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 i = 0; i < m_ietay.shape(0); i++) {
+        for (ssize_t j = 0; j < m_ietay.shape(1); j++) {
+            for (ssize_t k = 0; k < m_ietay.shape(2); k++) {
+                auto val = m_ietay(i, m_ietay.shape(1) - 1, k);
+                double norm = val < 1 ? 1 : val;
+                m_ietay(i, j, k) /= norm;
+            }
+        }
+    }
+}
+
+std::vector<Photon> Interpolator::interpolate(const ClusterVector<int32_t>& clusters) {
+    std::vector<Photon> photons;
+    photons.reserve(clusters.size());
+
+    if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
+        for (size_t i = 0; i<clusters.size(); i++){
+
+            auto cluster = clusters.at<Cluster3x3>(i);
+            Eta2 eta= calculate_eta2(cluster);
+    
+            Photon photon;
+            photon.x = cluster.x;
+            photon.y = cluster.y;
+            photon.energy = eta.sum;
+    
+            // auto ie = nearest_index(m_energy_bins, photon.energy)-1;
+            // auto ix = nearest_index(m_etabinsx, eta.x)-1;
+            // auto iy = nearest_index(m_etabinsy, eta.y)-1;   
+            //Finding the index of the last element that is smaller
+            //should work fine as long as we have many bins
+            auto ie = last_smaller(m_energy_bins, photon.energy);
+            auto ix = last_smaller(m_etabinsx, eta.x);
+            auto iy = last_smaller(m_etabinsy, eta.y); 
+
+            // fmt::print("ex: {}, ix: {}, iy: {}\n", ie, ix, iy);
+            
+            double dX, dY;
+            int ex, ey;
+            // cBottomLeft = 0,
+            // cBottomRight = 1,
+            // cTopLeft = 2,
+            // cTopRight = 3
+            switch (eta.c) {
+            case cTopLeft:
+                dX = -1.;
+                dY = 0.;
+                break;
+            case cTopRight:;
+                dX = 0.;
+                dY = 0.;
+                break;
+            case cBottomLeft:
+                dX = -1.;
+                dY = -1.;
+                break;
+            case cBottomRight:
+                dX = 0.; 
+                dY = -1.;
+                break;
+            }
+            photon.x += m_ietax(ix, iy, ie)*2 + dX;
+            photon.y += m_ietay(ix, iy, ie)*2 + dY;
+            photons.push_back(photon);
+        }
+    }else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 2){
+        for (size_t i = 0; i<clusters.size(); i++){
+            auto cluster = clusters.at<Cluster2x2>(i);
+            Eta2 eta= calculate_eta2(cluster);
+    
+            Photon photon;
+            photon.x = cluster.x;
+            photon.y = cluster.y;
+            photon.energy = eta.sum;
+    
+            //Now do some actual interpolation. 
+            //Find which energy bin the cluster is in
+            // auto ie = nearest_index(m_energy_bins, photon.energy)-1;
+            // auto ix = nearest_index(m_etabinsx, eta.x)-1;
+            // auto iy = nearest_index(m_etabinsy, eta.y)-1;  
+            //Finding the index of the last element that is smaller
+            //should work fine as long as we have many bins
+            auto ie = last_smaller(m_energy_bins, photon.energy);
+            auto ix = last_smaller(m_etabinsx, eta.x);
+            auto iy = last_smaller(m_etabinsy, eta.y); 
+
+            photon.x += m_ietax(ix, iy, ie)*2; //eta goes between 0 and 1 but we could move the hit anywhere in the 2x2
+            photon.y += m_ietay(ix, iy, ie)*2;
+            photons.push_back(photon);
+        }
+        
+    }else{
+        throw std::runtime_error("Only 3x3 and 2x2 clusters are supported for interpolation");
+    }
+  
+
+    return photons;
+}
+
+} // namespace aare
--- a/src/NDArray.test.cpp
+++ b/src/NDArray.test.cpp
@ -2,6 +2,7 @@
 #include <array>
 #include <catch2/benchmark/catch_benchmark.hpp>
 #include <catch2/catch_test_macros.hpp>
+#include <numeric>

 using aare::NDArray;
 using aare::NDView;
@ -34,6 +35,24 @@ TEST_CASE("Construct from an NDView") {
    }
 }

+TEST_CASE("3D NDArray from NDView"){
+    std::vector<int> data(27);
+    std::iota(data.begin(), data.end(), 0);
+    NDView<int, 3> view(data.data(), Shape<3>{3, 3, 3});
+    NDArray<int, 3> image(view);
+    REQUIRE(image.shape() == view.shape());
+    REQUIRE(image.size() == view.size());
+    REQUIRE(image.data() != view.data());
+
+    for(int64_t i=0; i<image.shape(0); i++){
+        for(int64_t j=0; j<image.shape(1); j++){
+            for(int64_t k=0; k<image.shape(2); k++){
+                REQUIRE(image(i, j, k) == view(i, j, k));
+            }
+        }
+    }
+}
+
 TEST_CASE("1D image") {
    std::array<int64_t, 1> shape{{20}};
    NDArray<short, 1> img(shape, 3);
--- a/src/RawFile.cpp
+++ b/src/RawFile.cpp
@ -76,8 +76,7 @@ size_t RawFile::n_mod() const { return n_subfile_parts; }


 size_t RawFile::bytes_per_frame() {
-    // return m_rows * m_cols * m_master.bitdepth() / 8;
-    return m_geometry.pixels_x * m_geometry.pixels_y * m_master.bitdepth() / 8;
+    return m_geometry.pixels_x * m_geometry.pixels_y * m_master.bitdepth() / bits_per_byte;
 }
 size_t RawFile::pixels_per_frame() { 
    // return m_rows * m_cols; 
--- a/src/algorithm.test.cpp
+++ b/src/algorithm.test.cpp
@ -0,0 +1,73 @@
+
+
+#include <catch2/catch_test_macros.hpp>
+#include <aare/algorithm.hpp>
+
+
+TEST_CASE("Find the closed index in a 1D array", "[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::nearest_index(arr, 2.3) == 2);
+    REQUIRE(aare::nearest_index(arr, 2.6) == 3);
+    REQUIRE(aare::nearest_index(arr, 45.0) == 4);
+    REQUIRE(aare::nearest_index(arr, 0.0) == 0);
+    REQUIRE(aare::nearest_index(arr, -1.0) == 0);
+}
+
+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;
+    }
+    // arr 0, 1, 2, 3, 4
+    REQUIRE(aare::nearest_index(arr, 2) == 2);
+    REQUIRE(aare::nearest_index(arr, 3) == 3);
+    REQUIRE(aare::nearest_index(arr, 45) == 4);
+    REQUIRE(aare::nearest_index(arr, 0) == 0);
+    REQUIRE(aare::nearest_index(arr, -1) == 0);
+}
+
+
+TEST_CASE("nearest_index works with std::vector", "[algorithm]"){
+    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);
+    REQUIRE(aare::nearest_index(vec, 4555555.0) == 4);
+    REQUIRE(aare::nearest_index(vec, 0.0) == 0);
+    REQUIRE(aare::nearest_index(vec, -10.0) == 0);
+}
+
+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);
+    REQUIRE(aare::nearest_index(arr, 4555555.0) == 4);
+    REQUIRE(aare::nearest_index(arr, 0.0) == 0);
+    REQUIRE(aare::nearest_index(arr, -10.0) == 0);
+}
+
+
+TEST_CASE("last 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, -10.0) == 0);
+    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);
+
+}