Fix/adapt and test interpolation (#231)

Adapted eta interpolation: 

### Issues with previous interpolation: 

## Eta Calculation: 
- previously assumed photon hit to be in bottom left pixel of cluster
(photon hit assumed in bottom right pixel of cluster)
- clusters are filled from top left to bottom right (previously assumed:
bottom left to top right)

## Actual Interpolation: 
- photon hits are given in pixel coordinates (previous interpolation
assumed euclidean coordinates, e.g. positive distance in y coordinate
becomes negative distance in row pixels)
- removed *2 of calculated distance 

## General Adaption: 
- max_sum_2x2 return subcluster index relative to cluster center e.g.
bottomleft, bottomright

## Added proper test case 
- simulated photon hit with normal energy distribution 
- Note: Test case for 2x2 cluster fails - Think uniform photon hit
distribution cant be modeled by normalized eta distribution for 2x2
clusters

include/aare/CalculateEta.hpp

@@ -7,10 +7,10 @@
 namespace aare {
 
 enum class corner : int {
-    cBottomLeft = 0,
-    cBottomRight = 1,
-    cTopLeft = 2,
-    cTopRight = 3
+    cTopLeft = 0,
+    cTopRight = 1,
+    cBottomLeft = 2,
+    cBottomRight = 3
 };
 
 enum class pixel : int {
@@ -58,90 +58,126 @@ template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType>
 Eta2<T>
 calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
-    Eta2<T> eta{};
 
-    auto max_sum = cl.max_sum_2x2();
-    eta.sum = max_sum.first;
-    auto c = max_sum.second;
+    static_assert(ClusterSizeX > 1 && ClusterSizeY > 1);
+    Eta2<T> eta{};
 
     size_t cluster_center_index =
         (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
 
-    size_t index_bottom_left_max_2x2_subcluster =
-        (int(c / (ClusterSizeX - 1))) * ClusterSizeX + c % (ClusterSizeX - 1);
+    auto max_sum = cl.max_sum_2x2();
+    eta.sum = max_sum.first;
+    int c = max_sum.second;
 
-    // calculate direction of gradient
-
-    // check that cluster center is in max subcluster
-    if (cluster_center_index != index_bottom_left_max_2x2_subcluster &&
-        cluster_center_index != index_bottom_left_max_2x2_subcluster + 1 &&
-        cluster_center_index !=
-            index_bottom_left_max_2x2_subcluster + ClusterSizeX &&
-        cluster_center_index !=
-            index_bottom_left_max_2x2_subcluster + ClusterSizeX + 1)
-        throw std::runtime_error("Photon center is not in max 2x2_subcluster");
-
-    if ((cluster_center_index - index_bottom_left_max_2x2_subcluster) %
-            ClusterSizeX ==
-        0) {
-        if ((cl.data[cluster_center_index + 1] +
+    // subcluster top right from center
+    switch (static_cast<corner>(c)) {
+    case (corner::cTopLeft):
+        if ((cl.data[cluster_center_index - 1] +
              cl.data[cluster_center_index]) != 0)
+            eta.x = static_cast<double>(cl.data[cluster_center_index - 1]) /
+                    static_cast<double>(cl.data[cluster_center_index - 1] +
+                                        cl.data[cluster_center_index]);
+        if ((cl.data[cluster_center_index - ClusterSizeX] +
+             cl.data[cluster_center_index]) != 0)
+            eta.y = static_cast<double>(
+                        cl.data[cluster_center_index - ClusterSizeX]) /
+                    static_cast<double>(
+                        cl.data[cluster_center_index - ClusterSizeX] +
+                        cl.data[cluster_center_index]);
 
-            eta.x = static_cast<double>(cl.data[cluster_center_index + 1]) /
-                    static_cast<double>((cl.data[cluster_center_index + 1] +
-                                         cl.data[cluster_center_index]));
-    } else {
-        if ((cl.data[cluster_center_index] +
-             cl.data[cluster_center_index - 1]) != 0)
-
+        // dx = 0
+        // dy = 0
+        break;
+    case (corner::cTopRight):
+        if (cl.data[cluster_center_index] + cl.data[cluster_center_index + 1] !=
+            0)
             eta.x = static_cast<double>(cl.data[cluster_center_index]) /
-                    static_cast<double>((cl.data[cluster_center_index - 1] +
-                                         cl.data[cluster_center_index]));
-    }
-    if ((cluster_center_index - index_bottom_left_max_2x2_subcluster) /
-            ClusterSizeX <
-        1) {
-        assert(cluster_center_index + ClusterSizeX <
-               ClusterSizeX * ClusterSizeY); // suppress warning
+                    static_cast<double>(cl.data[cluster_center_index] +
+                                        cl.data[cluster_center_index + 1]);
+        if ((cl.data[cluster_center_index - ClusterSizeX] +
+             cl.data[cluster_center_index]) != 0)
+            eta.y = static_cast<double>(
+                        cl.data[cluster_center_index - ClusterSizeX]) /
+                    static_cast<double>(
+                        cl.data[cluster_center_index - ClusterSizeX] +
+                        cl.data[cluster_center_index]);
+        // dx = 1
+        // dy = 0
+        break;
+    case (corner::cBottomLeft):
+        if ((cl.data[cluster_center_index - 1] +
+             cl.data[cluster_center_index]) != 0)
+            eta.x = static_cast<double>(cl.data[cluster_center_index - 1]) /
+                    static_cast<double>(cl.data[cluster_center_index - 1] +
+                                        cl.data[cluster_center_index]);
         if ((cl.data[cluster_center_index] +
              cl.data[cluster_center_index + ClusterSizeX]) != 0)
-            eta.y = static_cast<double>(
-                        cl.data[cluster_center_index + ClusterSizeX]) /
-                    static_cast<double>(
-                        (cl.data[cluster_center_index] +
-                         cl.data[cluster_center_index + ClusterSizeX]));
-    } else {
-        if ((cl.data[cluster_center_index] +
-             cl.data[cluster_center_index - ClusterSizeX]) != 0)
             eta.y = static_cast<double>(cl.data[cluster_center_index]) /
                     static_cast<double>(
-                        (cl.data[cluster_center_index] +
-                         cl.data[cluster_center_index - ClusterSizeX]));
+                        cl.data[cluster_center_index] +
+                        cl.data[cluster_center_index + ClusterSizeX]);
+        // dx = 0
+        // dy = 1
+        break;
+    case (corner::cBottomRight):
+        if (cl.data[cluster_center_index] + cl.data[cluster_center_index + 1] !=
+            0)
+            eta.x = static_cast<double>(cl.data[cluster_center_index]) /
+                    static_cast<double>(cl.data[cluster_center_index] +
+                                        cl.data[cluster_center_index + 1]);
+        if ((cl.data[cluster_center_index] +
+             cl.data[cluster_center_index + ClusterSizeX]) != 0)
+            eta.y = static_cast<double>(cl.data[cluster_center_index]) /
+                    static_cast<double>(
+                        cl.data[cluster_center_index] +
+                        cl.data[cluster_center_index + ClusterSizeX]);
+        // dx = 1
+        // dy = 1
+        break;
     }
 
-    eta.c = c; // TODO only supported for 2x2 and 3x3 clusters -> at least no
-               // underyling enum class
+    eta.c = c;
+
     return eta;
 }
 
-// TODO! Look up eta2 calculation - photon center should be top right corner
+// TODO! Look up eta2 calculation - photon center should be bottom right corner
 template <typename T>
 Eta2<T> calculate_eta2(const Cluster<T, 2, 2, int16_t> &cl) {
     Eta2<T> eta{};
 
     if ((cl.data[0] + cl.data[1]) != 0)
-        eta.x = static_cast<double>(cl.data[1]) /
-                (cl.data[0] + cl.data[1]); // between (0,1) the closer to zero
+        eta.x = static_cast<double>(cl.data[2]) /
+                (cl.data[2] + cl.data[3]); // between (0,1) the closer to zero
                                            // left value probably larger
     if ((cl.data[0] + cl.data[2]) != 0)
-        eta.y = static_cast<double>(cl.data[2]) /
-                (cl.data[0] + cl.data[2]); // between (0,1) the closer to zero
+        eta.y = static_cast<double>(cl.data[1]) /
+                (cl.data[1] + cl.data[3]); // between (0,1) the closer to zero
                                            // bottom value probably larger
     eta.sum = cl.sum();
 
     return eta;
 }
 
+// TODO generalize
+template <typename T>
+Eta2<T> calculate_eta2(const Cluster<T, 1, 2, int16_t> &cl) {
+    Eta2<T> eta{};
+
+    eta.x = 0;
+    eta.y = static_cast<double>(cl.data[0]) / cl.data[1];
+    eta.sum = cl.sum();
+}
+
+template <typename T>
+Eta2<T> calculate_eta2(const Cluster<T, 2, 1, int16_t> &cl) {
+    Eta2<T> eta{};
+
+    eta.x = static_cast<double>(cl.data[0]) / cl.data[1];
+    eta.y = 0;
+    eta.sum = cl.sum();
+}
+
 // calculates Eta3 for 3x3 cluster based on code from analyze_cluster
 // TODO only supported for 3x3 Clusters
 template <typename T> Eta2<T> calculate_eta3(const Cluster<T, 3, 3> &cl) {

include/aare/Cluster.hpp

@@ -8,7 +8,6 @@
 
 #pragma once
 
-#include "logger.hpp"
 #include <algorithm>
 #include <array>
 #include <cstdint>
@@ -19,7 +18,7 @@ namespace aare {
 
 // requires clause c++20 maybe update
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType = int16_t>
+          typename CoordType = uint16_t>
 struct Cluster {
 
     static_assert(std::is_arithmetic_v<T>, "T needs to be an arithmetic type");
@@ -39,6 +38,13 @@ struct Cluster {
 
     T sum() const { return std::accumulate(data.begin(), data.end(), T{}); }
 
+    // TODO: handle 1 dimensional clusters
+    // TODO: change int to corner
+    /**
+     * @brief sum of 2x2 subcluster with highest energy
+     * @return photon energy of subcluster, 2x2 subcluster index relative to
+     * cluster center
+     */
     std::pair<T, int> max_sum_2x2() const {
 
         if constexpr (cluster_size_x == 3 && cluster_size_y == 3) {
@@ -54,17 +60,38 @@ struct Cluster {
         } else if constexpr (cluster_size_x == 2 && cluster_size_y == 2) {
             return std::make_pair(data[0] + data[1] + data[2] + data[3], 0);
         } else {
-            constexpr size_t num_2x2_subclusters =
-                (ClusterSizeX - 1) * (ClusterSizeY - 1);
+            constexpr size_t cluster_center_index =
+                (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
 
-            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];
+            std::array<T, 4> sum_2x2_subcluster{0};
+            // subcluster top left from center
+            sum_2x2_subcluster[0] =
+                data[cluster_center_index] + data[cluster_center_index - 1] +
+                data[cluster_center_index - ClusterSizeX] +
+                data[cluster_center_index - 1 - ClusterSizeX];
+            // subcluster top right from center
+            if (ClusterSizeX > 2) {
+                sum_2x2_subcluster[1] =
+                    data[cluster_center_index] +
+                    data[cluster_center_index + 1] +
+                    data[cluster_center_index - ClusterSizeX] +
+                    data[cluster_center_index - ClusterSizeX + 1];
+            }
+            // subcluster bottom left from center
+            if (ClusterSizeY > 2) {
+                sum_2x2_subcluster[2] =
+                    data[cluster_center_index] +
+                    data[cluster_center_index - 1] +
+                    data[cluster_center_index + ClusterSizeX] +
+                    data[cluster_center_index + ClusterSizeX - 1];
+            }
+            // subcluster bottom right from center
+            if (ClusterSizeX > 2 && ClusterSizeY > 2) {
+                sum_2x2_subcluster[3] =
+                    data[cluster_center_index] +
+                    data[cluster_center_index + 1] +
+                    data[cluster_center_index + ClusterSizeX] +
+                    data[cluster_center_index + ClusterSizeX + 1];
             }
 
             int index = std::max_element(sum_2x2_subcluster.begin(),

include/aare/ClusterFinder.hpp

@@ -136,7 +136,7 @@ class ClusterFinder {
                     // don't have a photon
                     int i = 0;
                     for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
-                        for (int ic = -dx; ic < dx + has_center_pixel_y; ic++) {
+                        for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) {
                             if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
                                 iy + ir >= 0 && iy + ir < frame.shape(0)) {
                                 CT tmp =
@@ -144,7 +144,7 @@ class ClusterFinder {
                                     static_cast<CT>(
                                         m_pedestal.mean(iy + ir, ix + ic));
                                 cluster.data[i] =
-                                    tmp; // Watch for out of bounds access                                
+                                    tmp; // Watch for out of bounds access
                             }
                             i++;
                         }

include/aare/Interpolator.hpp

@@ -69,26 +69,27 @@ Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
             // cBottomRight = 1,
             // cTopLeft = 2,
             // cTopRight = 3
+            // TODO: could also chaneg the sign of the eta calculation
             switch (static_cast<corner>(eta.c)) {
             case corner::cTopLeft:
-                dX = -1.;
-                dY = 0;
+                dX = 0.0;
+                dY = 0.0;
                 break;
             case corner::cTopRight:;
-                dX = 0;
-                dY = 0;
+                dX = 1.0;
+                dY = 0.0;
                 break;
             case corner::cBottomLeft:
-                dX = -1.;
-                dY = -1.;
+                dX = 0.0;
+                dY = 1.0;
                 break;
             case corner::cBottomRight:
-                dX = 0.;
-                dY = -1.;
+                dX = 1.0;
+                dY = 1.0;
                 break;
             }
-            photon.x += m_ietax(ix, iy, ie) * 2 + dX;
-            photon.y += m_ietay(ix, iy, ie) * 2 + dY;
+            photon.x -= m_ietax(ix, iy, ie) - dX;
+            photon.y -= m_ietay(ix, iy, ie) - dY;
             photons.push_back(photon);
         }
     } else if (clusters.cluster_size_x() == 2 ||
@@ -112,10 +113,11 @@ Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
             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;
+            // TODO: why 2?
+            photon.x -=
+                m_ietax(ix, iy, ie); // eta goes between 0 and 1 but we could
+                                     // move the hit anywhere in the 2x2
+            photon.y -= m_ietay(ix, iy, ie);
             photons.push_back(photon);
         }
 

python/src/bind_Cluster.hpp

@@ -24,7 +24,7 @@ void define_Cluster(py::module &m, const std::string &typestr) {
     py::class_<Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType>>(
         m, class_name.c_str(), py::buffer_protocol())
 
-        .def(py::init([](uint8_t x, uint8_t y,
+        .def(py::init([](CoordType x, CoordType y,
                          py::array_t<Type, py::array::forcecast> data) {
             py::buffer_info buf_info = data.request();
             Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType> cluster;

python/src/bind_ClusterFile.hpp

@@ -44,14 +44,14 @@ void define_ClusterFile(py::module &m, const std::string &typestr) {
                  auto v = new ClusterVector<ClusterType>(self.read_frame());
                  return v;
              })
-        .def("set_roi", &ClusterFile<ClusterType>::set_roi,
-             py::arg("roi"))
+        .def("set_roi", &ClusterFile<ClusterType>::set_roi, py::arg("roi"))
         .def(
             "set_noise_map",
             [](ClusterFile<ClusterType> &self, py::array_t<int32_t> noise_map) {
                 auto view = make_view_2d(noise_map);
                 self.set_noise_map(view);
-            }, py::arg("noise_map"))
+            },
+            py::arg("noise_map"))
 
         .def("set_gain_map",
              [](ClusterFile<ClusterType> &self, py::array_t<double> gain_map) {
@@ -84,11 +84,19 @@ template <typename Type, uint8_t CoordSizeX, uint8_t CoordSizeY,
           typename CoordType = uint16_t>
 void register_calculate_eta(py::module &m) {
     using ClusterType = Cluster<Type, CoordSizeX, CoordSizeY, CoordType>;
+
     m.def("calculate_eta2",
           [](const aare::ClusterVector<ClusterType> &clusters) {
               auto eta2 = new NDArray<double, 2>(calculate_eta2(clusters));
               return return_image_data(eta2);
           });
+
+    m.def("calculate_eta2", [](const aare::Cluster<Type, CoordSizeX, CoordSizeY,
+                                                   CoordType> &cluster) {
+        auto eta2 = calculate_eta2(cluster);
+        // TODO return proper eta class
+        return py::make_tuple(eta2.x, eta2.y, eta2.sum);
+    });
 }
 
 #pragma GCC diagnostic pop

python/src/module.cpp

@@ -87,6 +87,8 @@ PYBIND11_MODULE(_aare, m) {
     DEFINE_CLUSTER_BINDINGS(double, 9, 9, uint16_t, d);
     DEFINE_CLUSTER_BINDINGS(float, 9, 9, uint16_t, f);
 
+    // DEFINE_CLUSTER_BINDINGS(double, 2, 1, uint16_t, d);
+
     define_3x3_reduction<int, 3, 3, uint16_t>(m);
     define_3x3_reduction<double, 3, 3, uint16_t>(m);
     define_3x3_reduction<float, 3, 3, uint16_t>(m);

python/tests/test_Cluster.py

@@ -53,8 +53,8 @@ def test_Interpolator():
 
     assert interpolated_photons.size == 1
 
-    assert interpolated_photons[0]["x"] == -1 
-    assert interpolated_photons[0]["y"] == -1
+    assert interpolated_photons[0]["x"] == 0 
+    assert interpolated_photons[0]["y"] == 0
     assert interpolated_photons[0]["energy"] == 4 #eta_sum = 4, dx, dy = -1,-1 m_ietax = 0, m_ietay = 0
 
     clustervector = _aare.ClusterVector_Cluster2x2i()
@@ -84,7 +84,7 @@ def test_calculate_eta():
     assert eta2[0,0] == 0.5
     assert eta2[0,1] == 0.5
     assert eta2[1,0] == 0.5
-    assert eta2[1,1] == 0.6 #1/5
+    assert eta2[1,1] == 0.4 #2/5
 
 def test_cluster_finder(): 
     """Test ClusterFinder""" 

python/tests/test_Interpolation.py

@@ -0,0 +1,148 @@
+import pytest
+import numpy as np 
+import boost_histogram as bh 
+import pickle
+from scipy.stats import multivariate_normal
+
+from aare import Interpolator, calculate_eta2
+from aare._aare import ClusterVector_Cluster2x2d, Cluster2x2d, Cluster3x3d, ClusterVector_Cluster3x3d
+
+from conftest import test_data_path
+
+pixel_width = 1e-4
+values = np.arange(0.5*pixel_width, 0.1, pixel_width)
+num_pixels = values.size
+X, Y = np.meshgrid(values, values)
+data_points = np.stack([X.ravel(), Y.ravel()], axis=1)
+variance = 10*pixel_width
+covariance_matrix = np.array([[variance, 0],[0, variance]])
+
+
+
+def create_photon_hit_with_gaussian_distribution(mean, covariance_matrix, data_points):
+    gaussian = multivariate_normal(mean=mean, cov=covariance_matrix)
+    probability_values = gaussian.pdf(data_points)
+    return (probability_values.reshape(X.shape)).round() #python bindings only support frame types of uint16_t
+
+def photon_hit_in_euclidean_space(cluster_center, pixels_per_superpixel, photon_hit): 
+    scaled_photon_hit_x = cluster_center - (1 - photon_hit[0][0])*pixels_per_superpixel*pixel_width 
+    scaled_photon_hit_y = cluster_center - (1 - photon_hit[0][1])*pixels_per_superpixel*pixel_width
+    return (scaled_photon_hit_x, scaled_photon_hit_y)
+
+def create_2x2cluster_from_frame(frame, pixels_per_superpixel):
+    return Cluster2x2d(1, 1, np.array([frame[0:pixels_per_superpixel, 0:pixels_per_superpixel].sum(), 
+                                            frame[0:pixels_per_superpixel, pixels_per_superpixel:2*pixels_per_superpixel].sum(), 
+                                            frame[pixels_per_superpixel:2*pixels_per_superpixel, 0:pixels_per_superpixel].sum(), 
+                                            frame[pixels_per_superpixel:2*pixels_per_superpixel, pixels_per_superpixel:2*pixels_per_superpixel].sum()], dtype=np.float64))
+
+
+def create_3x3cluster_from_frame(frame, pixels_per_superpixel):
+    return Cluster3x3d(1, 1, np.array([frame[0:pixels_per_superpixel, 0:pixels_per_superpixel].sum(), 
+                                    frame[0:pixels_per_superpixel, pixels_per_superpixel:2*pixels_per_superpixel].sum(), 
+                                    frame[0:pixels_per_superpixel, 2*pixels_per_superpixel:3*pixels_per_superpixel].sum(), 
+                                    frame[pixels_per_superpixel:2*pixels_per_superpixel, 0:pixels_per_superpixel].sum(), 
+                                    frame[pixels_per_superpixel:2*pixels_per_superpixel, pixels_per_superpixel:2*pixels_per_superpixel].sum(), 
+                                    frame[pixels_per_superpixel:2*pixels_per_superpixel, 2*pixels_per_superpixel:3*pixels_per_superpixel].sum(), 
+                                    frame[2*pixels_per_superpixel:3*pixels_per_superpixel, 0:pixels_per_superpixel].sum(), 
+                                    frame[2*pixels_per_superpixel:3*pixels_per_superpixel, pixels_per_superpixel:2*pixels_per_superpixel].sum(), 
+                                    frame[2*pixels_per_superpixel:3*pixels_per_superpixel, 2*pixels_per_superpixel:3*pixels_per_superpixel].sum()], dtype=np.float64))
+
+
+def calculate_eta_distribution(num_frames, pixels_per_superpixel, random_number_generator, cluster_2x2 = True): 
+    hist = bh.Histogram(
+    bh.axis.Regular(100, -0.2, 1.2), 
+    bh.axis.Regular(100, -0.2, 1.2), bh.axis.Regular(1, 0, num_pixels*num_pixels*1/(variance*2*np.pi))) 
+
+    for _ in range(0, num_frames):
+        mean_x = random_number_generator.uniform(pixels_per_superpixel*pixel_width, 2*pixels_per_superpixel*pixel_width)
+        mean_y = random_number_generator.uniform(pixels_per_superpixel*pixel_width, 2*pixels_per_superpixel*pixel_width)
+        frame = create_photon_hit_with_gaussian_distribution(np.array([mean_x, mean_y]), variance, data_points)
+
+        cluster = None
+
+        if cluster_2x2:
+            cluster = create_2x2cluster_from_frame(frame, pixels_per_superpixel)
+        else: 
+            cluster = create_3x3cluster_from_frame(frame, pixels_per_superpixel)
+
+        eta2 = calculate_eta2(cluster)
+        hist.fill(eta2[0], eta2[1], eta2[2])
+
+    return hist
+
+@pytest.mark.withdata
+def test_interpolation_of_2x2_cluster(test_data_path): 
+    """Test Interpolation of 2x2 cluster from Photon hit with Gaussian Distribution""" 
+
+    #TODO maybe better to compute in test instead of loading - depends on eta 
+    """ 
+    filename = test_data_path/"eta_distributions"/"eta_distribution_2x2cluster_gaussian.pkl"
+    with open(filename, "rb") as f:
+        eta_distribution = pickle.load(f)
+    """
+
+    num_frames = 1000
+    pixels_per_superpixel = int(num_pixels*0.5)
+    random_number_generator = np.random.default_rng(42)
+
+    eta_distribution = calculate_eta_distribution(num_frames, pixels_per_superpixel, random_number_generator)
+
+    interpolation = Interpolator(eta_distribution, eta_distribution.axes[0].edges[:-1], eta_distribution.axes[1].edges[:-1], eta_distribution.axes[2].edges[:-1])
+
+    #actual photon hit
+    mean = 1.2*pixels_per_superpixel*pixel_width
+    mean = np.array([mean, mean])
+    frame = create_photon_hit_with_gaussian_distribution(mean, covariance_matrix, data_points)
+    cluster = create_2x2cluster_from_frame(frame, pixels_per_superpixel)
+
+    clustervec = ClusterVector_Cluster2x2d()
+    clustervec.push_back(cluster)
+
+    interpolated_photon = interpolation.interpolate(clustervec)
+
+    assert interpolated_photon.size == 1
+
+    cluster_center = 1.5*pixels_per_superpixel*pixel_width
+
+    scaled_photon_hit = photon_hit_in_euclidean_space(cluster_center, pixels_per_superpixel, interpolated_photon)
+
+    assert (np.linalg.norm(scaled_photon_hit - mean) < np.linalg.norm(np.array([cluster_center, cluster_center] - mean)))
+
+
+@pytest.mark.withdata
+def test_interpolation_of_3x3_cluster(test_data_path): 
+    """Test Interpolation of 3x3 Cluster from Photon hit with Gaussian Distribution""" 
+
+    #TODO maybe better to compute in test instead of loading - depends on eta 
+    """ 
+    filename = test_data_path/"eta_distributions"/"eta_distribution_3x3cluster_gaussian.pkl"
+    with open(filename, "rb") as f:
+        eta_distribution = pickle.load(f)
+    """
+
+    num_frames = 1000
+    pixels_per_superpixel = int(num_pixels/3)
+    random_number_generator = np.random.default_rng(42)
+    eta_distribution = calculate_eta_distribution(num_frames, pixels_per_superpixel, random_number_generator, False)
+
+    interpolation = Interpolator(eta_distribution, eta_distribution.axes[0].edges[:-1], eta_distribution.axes[1].edges[:-1], eta_distribution.axes[2].edges[:-1])
+
+    #actual photon hit
+    mean = 1.2*pixels_per_superpixel*pixel_width
+    mean = np.array([mean, mean])
+    frame = create_photon_hit_with_gaussian_distribution(mean, covariance_matrix, data_points)
+    cluster = create_3x3cluster_from_frame(frame, pixels_per_superpixel)
+
+    clustervec = ClusterVector_Cluster3x3d()
+    clustervec.push_back(cluster)
+
+    interpolated_photon = interpolation.interpolate(clustervec)
+
+    assert interpolated_photon.size == 1
+
+    cluster_center = 1.5*pixels_per_superpixel*pixel_width
+
+    scaled_photon_hit = photon_hit_in_euclidean_space(cluster_center, pixels_per_superpixel, interpolated_photon)
+
+    assert (np.linalg.norm(scaled_photon_hit - mean) < np.linalg.norm(np.array([cluster_center, cluster_center] - mean)))
+

src/ClusterFile.test.cpp

@@ -279,7 +279,8 @@ TEST_CASE("Read cluster from multiple frame file", "[.with-data]") {
     }
 }
 
-TEST_CASE("Write cluster with potential padding", "[.with-data][.ClusterFile]") {
+TEST_CASE("Write cluster with potential padding",
+          "[.with-data][.ClusterFile]") {
 
     using ClusterType = Cluster<double, 3, 3>;
 
@@ -290,7 +291,7 @@ TEST_CASE("Write cluster with potential padding", "[.with-data][.ClusterFile]")
     ClusterFile<ClusterType> file(fpath, 1000, "w");
 
     ClusterVector<ClusterType> clustervec(2);
-    int16_t coordinate = 5;
+    uint16_t coordinate = 5;
     clustervec.push_back(ClusterType{
         coordinate, coordinate, {0., 0., 0., 0., 0., 0., 0., 0., 0.}});
     clustervec.push_back(ClusterType{