added some python tests
Some checks failed
Build the package using cmake then documentation / build (ubuntu-latest, 3.12) (push) Failing after 41s

This commit is contained in:
Mazzoleni Alice Francesca 2025-04-04 17:19:15 +02:00
parent 885309d97c
commit 9de84a7f87
12 changed files with 264 additions and 151 deletions

View File

@ -40,6 +40,7 @@ struct Cluster {
constexpr size_t num_2x2_subclusters = constexpr size_t num_2x2_subclusters =
(ClusterSizeX - 1) * (ClusterSizeY - 1); (ClusterSizeX - 1) * (ClusterSizeY - 1);
std::array<T, num_2x2_subclusters> sum_2x2_subcluster; std::array<T, num_2x2_subclusters> sum_2x2_subcluster;
for (size_t i = 0; i < ClusterSizeY - 1; ++i) { for (size_t i = 0; i < ClusterSizeY - 1; ++i) {
for (size_t j = 0; j < ClusterSizeX - 1; ++j) for (size_t j = 0; j < ClusterSizeX - 1; ++j)

View File

@ -36,7 +36,7 @@ uint32_t number_of_clusters
* etc. * etc.
*/ */
template <typename ClusterType, template <typename ClusterType,
typename Enable = std::enable_if_t<is_cluster_v<ClusterType>, bool>> typename Enable = std::enable_if_t<is_cluster_v<ClusterType>>>
class ClusterFile { class ClusterFile {
FILE *fp{}; FILE *fp{};
uint32_t m_num_left{}; /*Number of photons left in frame*/ uint32_t m_num_left{}; /*Number of photons left in frame*/
@ -70,8 +70,6 @@ class ClusterFile {
*/ */
ClusterVector<ClusterType> read_clusters(size_t n_clusters); ClusterVector<ClusterType> read_clusters(size_t n_clusters);
ClusterVector<ClusterType> read_clusters(size_t n_clusters, ROI roi);
/** /**
* @brief Read a single frame from the file and return the clusters. The * @brief Read a single frame from the file and return the clusters. The
* cluster vector will have the frame number set. * cluster vector will have the frame number set.

View File

@ -133,6 +133,7 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
* @brief Sum the pixels in each cluster * @brief Sum the pixels in each cluster
* @return std::vector<T> vector of sums for each cluster * @return std::vector<T> vector of sums for each cluster
*/ */
/*
std::vector<T> sum() { std::vector<T> sum() {
std::vector<T> sums(m_size); std::vector<T> sums(m_size);
const size_t stride = item_size(); const size_t stride = item_size();
@ -147,12 +148,14 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
} }
return sums; return sums;
} }
*/
/** /**
* @brief Sum the pixels in the 2x2 subcluster with the biggest pixel sum in * @brief Sum the pixels in the 2x2 subcluster with the biggest pixel sum in
* each cluster * each cluster
* @return std::vector<T> vector of sums for each cluster * @return std::vector<T> vector of sums for each cluster
*/ //TODO if underlying container is a vector use std::for_each */ //TODO if underlying container is a vector use std::for_each
/*
std::vector<T> sum_2x2() { std::vector<T> sum_2x2() {
std::vector<T> sums_2x2(m_size); std::vector<T> sums_2x2(m_size);
@ -161,6 +164,7 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
} }
return sums_2x2; return sums_2x2;
} }
*/
/** /**
* @brief Return the number of clusters in the vector * @brief Return the number of clusters in the vector

View File

@ -1,10 +1,12 @@
#pragma once #pragma once
#include "aare/CalculateEta.hpp"
#include "aare/Cluster.hpp" #include "aare/Cluster.hpp"
#include "aare/ClusterFile.hpp" //Cluster_3x3 #include "aare/ClusterFile.hpp" //Cluster_3x3
#include "aare/ClusterVector.hpp" #include "aare/ClusterVector.hpp"
#include "aare/NDArray.hpp" #include "aare/NDArray.hpp"
#include "aare/NDView.hpp" #include "aare/NDView.hpp"
#include "aare/algorithm.hpp"
namespace aare { namespace aare {
@ -28,8 +30,103 @@ class Interpolator {
NDArray<double, 3> get_ietax() { return m_ietax; } NDArray<double, 3> get_ietax() { return m_ietax; }
NDArray<double, 3> get_ietay() { return m_ietay; } NDArray<double, 3> get_ietay() { return m_ietay; }
template <typename ClusterType> template <typename ClusterType,
typename Eanble = std::enable_if_t<is_cluster_v<ClusterType>>>
std::vector<Photon> interpolate(const ClusterVector<ClusterType> &clusters); std::vector<Photon> interpolate(const ClusterVector<ClusterType> &clusters);
}; };
// TODO: generalize to support any clustertype!!! otherwise add std::enable_if_t
// to only take Cluster2x2 and Cluster3x3
template <typename ClusterType, typename Enable>
std::vector<Photon>
Interpolator::interpolate(const ClusterVector<ClusterType> &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(i);
auto 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;
// 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(i);
auto 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 } // namespace aare

View File

@ -18,20 +18,81 @@ using pd_type = double;
using namespace aare; using namespace aare;
template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
void define_cluster(py::module &m, const std::string &typestr) {
auto class_name = fmt::format("Cluster{}", typestr);
using ClusterType =
Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType, void>;
py::class_<Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType, void>>(
m, class_name.c_str())
.def(py::init([](uint8_t x, uint8_t y, py::array_t<Type> data) {
py::buffer_info buf_info = data.request();
Type *ptr = static_cast<Type *>(buf_info.ptr);
Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType, void> cluster;
cluster.x = x;
cluster.y = y;
std::copy(ptr, ptr + ClusterSizeX * ClusterSizeY,
cluster.data); // Copy array contents
return cluster;
}))
//.def(py::init<>())
.def_readwrite("x", &ClusterType::x)
.def_readwrite("y", &ClusterType::y)
.def_property(
"data",
[](ClusterType &c) -> py::array {
return py::array(py::buffer_info(
c.data, sizeof(Type),
py::format_descriptor<Type>::format(), // Type
// format
1, // Number of dimensions
{static_cast<ssize_t>(ClusterSizeX *
ClusterSizeY)}, // Shape (flattened)
{sizeof(Type)} // Stride (step size between elements)
));
},
[](ClusterType &c, py::array_t<Type> arr) {
py::buffer_info buf_info = arr.request();
Type *ptr = static_cast<Type *>(buf_info.ptr);
std::copy(ptr, ptr + ClusterSizeX * ClusterSizeY, c.data);
});
}
template <typename ClusterType> template <typename ClusterType>
void define_cluster_vector(py::module &m, const std::string &typestr) { void define_cluster_vector(py::module &m, const std::string &typestr) {
auto class_name = fmt::format("ClusterVector_{}", typestr); auto class_name = fmt::format("ClusterVector_{}", typestr);
py::class_<ClusterVector<ClusterType>>(m, class_name.c_str(), py::class_<ClusterVector<ClusterType>>(m, class_name.c_str(),
py::buffer_protocol()) py::buffer_protocol())
.def(py::init<int, int>(), py::arg("cluster_size_x") = 3, .def(py::init()) // TODO change!!!
py::arg("cluster_size_y") = 3) // TODO change!!! /*
.def("push_back",
[](ClusterVector<ClusterType> &self, ClusterType &cl) {
// auto view = make_view_2d(data);
self.push_back(cl);
})
*/
/*
.def(
"push_back",
[](ClusterVector<ClusterType> &self, py::object obj) {
ClusterType &cl = py::cast<ClusterType &>(obj);
self.push_back(cl);
},
py::arg("cluster"))
*/
/*
.def("push_back", .def("push_back",
[](ClusterVector<ClusterType> &self, ClusterType &cl) { [](ClusterVector<ClusterType> &self, const ClusterType &cluster) {
// auto view = make_view_2d(data); self.push_back(cluster);
self.push_back(cl);
}) })
*/
//.def("push_back", &ClusterVector<ClusterType>::push_back) //TODO
//implement push_back
.def_property_readonly("size", &ClusterVector<ClusterType>::size) .def_property_readonly("size", &ClusterVector<ClusterType>::size)
.def("item_size", &ClusterVector<ClusterType>::item_size) .def("item_size", &ClusterVector<ClusterType>::item_size)
.def_property_readonly("fmt", .def_property_readonly("fmt",
@ -78,7 +139,6 @@ void define_cluster_vector(py::module &m, const std::string &typestr) {
template <typename ClusterType> template <typename ClusterType>
void define_cluster_finder_mt_bindings(py::module &m, void define_cluster_finder_mt_bindings(py::module &m,
const std::string &typestr) { const std::string &typestr) {
auto class_name = fmt::format("ClusterFinderMT_{}", typestr); auto class_name = fmt::format("ClusterFinderMT_{}", typestr);
py::class_<ClusterFinderMT<ClusterType, uint16_t, pd_type>>( py::class_<ClusterFinderMT<ClusterType, uint16_t, pd_type>>(
@ -129,7 +189,6 @@ void define_cluster_finder_mt_bindings(py::module &m,
template <typename ClusterType> template <typename ClusterType>
void define_cluster_collector_bindings(py::module &m, void define_cluster_collector_bindings(py::module &m,
const std::string &typestr) { const std::string &typestr) {
auto class_name = fmt::format("ClusterCollector_{}", typestr); auto class_name = fmt::format("ClusterCollector_{}", typestr);
py::class_<ClusterCollector<ClusterType>>(m, class_name.c_str()) py::class_<ClusterCollector<ClusterType>>(m, class_name.c_str())
@ -148,7 +207,6 @@ void define_cluster_collector_bindings(py::module &m,
template <typename ClusterType> template <typename ClusterType>
void define_cluster_file_sink_bindings(py::module &m, void define_cluster_file_sink_bindings(py::module &m,
const std::string &typestr) { const std::string &typestr) {
auto class_name = fmt::format("ClusterFileSink_{}", typestr); auto class_name = fmt::format("ClusterFileSink_{}", typestr);
py::class_<ClusterFileSink<ClusterType>>(m, class_name.c_str()) py::class_<ClusterFileSink<ClusterType>>(m, class_name.c_str())
@ -159,7 +217,6 @@ void define_cluster_file_sink_bindings(py::module &m,
template <typename ClusterType> template <typename ClusterType>
void define_cluster_finder_bindings(py::module &m, const std::string &typestr) { void define_cluster_finder_bindings(py::module &m, const std::string &typestr) {
auto class_name = fmt::format("ClusterFinder_{}", typestr); auto class_name = fmt::format("ClusterFinder_{}", typestr);
py::class_<ClusterFinder<ClusterType, uint16_t, pd_type>>( py::class_<ClusterFinder<ClusterType, uint16_t, pd_type>>(
@ -227,21 +284,4 @@ void define_cluster_finder_bindings(py::module &m, const std::string &typestr) {
} }
return hitmap; return hitmap;
}); });
py::class_<DynamicCluster>(m, "DynamicCluster", py::buffer_protocol())
.def(py::init<int, int, Dtype>())
.def("size", &DynamicCluster::size)
.def("begin", &DynamicCluster::begin)
.def("end", &DynamicCluster::end)
.def_readwrite("x", &DynamicCluster::x)
.def_readwrite("y", &DynamicCluster::y)
.def_buffer([](DynamicCluster &c) -> py::buffer_info {
return py::buffer_info(c.data(), c.dt.bytes(), c.dt.format_descr(),
1, {c.size()}, {c.dt.bytes()});
})
.def("__repr__", [](const DynamicCluster &a) {
return "<DynamicCluster: x: " + std::to_string(a.x) +
", y: " + std::to_string(a.y) + ">";
});
} }

View File

@ -39,14 +39,6 @@ void define_cluster_file_io_bindings(py::module &m,
return v; return v;
}, },
py::return_value_policy::take_ownership) py::return_value_policy::take_ownership)
.def(
"read_clusters",
[](ClusterFile<ClusterType> &self, size_t n_clusters, ROI roi) {
auto v = new ClusterVector<ClusterType>(
self.read_clusters(n_clusters, roi));
return v;
},
py::return_value_policy::take_ownership)
.def("read_frame", .def("read_frame",
[](ClusterFile<ClusterType> &self) { [](ClusterFile<ClusterType> &self) {
auto v = new ClusterVector<ClusterType>(self.read_frame()); auto v = new ClusterVector<ClusterType>(self.read_frame());

View File

@ -10,9 +10,9 @@
namespace py = pybind11; namespace py = pybind11;
template <typename ClusterType> template <typename ClusterType>
void register_interpolate(py::class_<aare::Interpolator> &interpolator) { void register_interpolate(py::class_<aare::Interpolator> &interpolator,
std::string name = const std::string &typestr) {
fmt::format("interpolate_{}", typeid(ClusterType).name()); auto name = fmt::format("interpolate_{}", typestr);
interpolator.def(name.c_str(), interpolator.def(name.c_str(),
[](aare::Interpolator &self, [](aare::Interpolator &self,
@ -50,12 +50,12 @@ void define_interpolation_bindings(py::module &m) {
return return_image_data(ptr); return return_image_data(ptr);
}); });
register_interpolate<Cluster<int, 3, 3>>(interpolator); register_interpolate<Cluster<int, 3, 3>>(interpolator, "Cluster3x3i");
register_interpolate<Cluster<float, 3, 3>>(interpolator); register_interpolate<Cluster<float, 3, 3>>(interpolator, "Cluster3x3f");
register_interpolate<Cluster<double, 3, 3>>(interpolator); register_interpolate<Cluster<double, 3, 3>>(interpolator, "Cluster3x3d");
register_interpolate<Cluster<int, 2, 2>>(interpolator); register_interpolate<Cluster<int, 2, 2>>(interpolator, "Cluster2x2i");
register_interpolate<Cluster<float, 2, 2>>(interpolator); register_interpolate<Cluster<float, 2, 2>>(interpolator, "Cluster2x2f");
register_interpolate<Cluster<double, 2, 2>>(interpolator); register_interpolate<Cluster<double, 2, 2>>(interpolator, "Cluster2x2d");
// TODO! Evaluate without converting to double // TODO! Evaluate without converting to double
m.def( m.def(

View File

@ -70,4 +70,11 @@ PYBIND11_MODULE(_aare, m) {
define_cluster_collector_bindings<Cluster<int, 2, 2>>(m, "Cluster2x2i"); define_cluster_collector_bindings<Cluster<int, 2, 2>>(m, "Cluster2x2i");
define_cluster_collector_bindings<Cluster<double, 2, 2>>(m, "Cluster2x2f"); define_cluster_collector_bindings<Cluster<double, 2, 2>>(m, "Cluster2x2f");
define_cluster_collector_bindings<Cluster<float, 2, 2>>(m, "Cluster2x2d"); define_cluster_collector_bindings<Cluster<float, 2, 2>>(m, "Cluster2x2d");
define_cluster<int, 3, 3, uint16_t>(m, "3x3i");
define_cluster<float, 3, 3, uint16_t>(m, "3x3f");
define_cluster<double, 3, 3, uint16_t>(m, "3x3d");
define_cluster<int, 2, 2, uint16_t>(m, "2x2i");
define_cluster<float, 2, 2, uint16_t>(m, "2x2f");
define_cluster<double, 2, 2, uint16_t>(m, "2x2d");
} }

View File

@ -40,25 +40,28 @@ template <typename T> py::array return_vector(std::vector<T> *vec) {
} }
// todo rewrite generic // todo rewrite generic
template <class T, int Flags> auto get_shape_3d(const 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)}; 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)); return aare::NDView<T, 3>(arr.mutable_data(), get_shape_3d<T, Flags>(arr));
} }
template <class T, int Flags> auto get_shape_2d(const 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)}; return aare::Shape<2>{arr.shape(0), arr.shape(1)};
} }
template <class T, int Flags> auto get_shape_1d(const 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)}; 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)); 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)); return aare::NDView<T, 1>(arr.mutable_data(), get_shape_1d<T, Flags>(arr));
} }

View File

@ -0,0 +1,64 @@
import pytest
import numpy as np
from _aare import ClusterVector_Cluster3x3i, Interpolator, Cluster3x3i, ClusterFinder_Cluster3x3i
def test_ClusterVector():
"""Test ClusterVector"""
clustervector = ClusterVector_Cluster3x3i()
assert clustervector.cluster_size_x == 3
assert clustervector.cluster_size_y == 3
assert clustervector.item_size() == 4+9*4
assert clustervector.frame_number == 0
assert clustervector.capacity == 1024
assert clustervector.size == 0
cluster = Cluster3x3i(0,0,np.ones(9, dtype=np.int32))
#clustervector.push_back(cluster)
#assert clustervector.size == 1
#push_back - check size
def test_Interpolator():
"""Test Interpolator"""
ebins = np.linspace(0,10, 20, dtype=np.float64)
xbins = np.linspace(0, 5, 30, dtype=np.float64)
ybins = np.linspace(0, 5, 30, dtype=np.float64)
etacube = np.zeros(shape=[30, 30, 20], dtype=np.float64)
interpolator = Interpolator(etacube, xbins, ybins, ebins)
assert interpolator.get_ietax().shape == (30,30,20)
assert interpolator.get_ietay().shape == (30,30,20)
clustervector = ClusterVector_Cluster3x3i()
#TODO clustervector is empty
cluster = Cluster3x3i(0,0, np.ones(9, dtype=np.int32))
#clustervector.push_back(cluster)
num_clusters = 1;
assert interpolator.interpolate_Cluster3x3i(clustervector).shape == (num_clusters, 3)
#def test_cluster_file():
#def test_cluster_finder():
#"""Test ClusterFinder"""
#clusterfinder = ClusterFinder_Cluster3x3i([100,100])
#clusterfinder.find_clusters()
#clusters = clusterfinder.steal_clusters()
#print("cluster size: ", clusters.size())

View File

@ -61,11 +61,13 @@ TEST_CASE("Summing 3x1 clusters of int64", "[.ClusterVector]") {
REQUIRE(cv.capacity() == 4); REQUIRE(cv.capacity() == 4);
REQUIRE(cv.size() == 3); REQUIRE(cv.size() == 3);
/*
auto sums = cv.sum(); auto sums = cv.sum();
REQUIRE(sums.size() == 3); REQUIRE(sums.size() == 3);
REQUIRE(sums[0] == 12); REQUIRE(sums[0] == 12);
REQUIRE(sums[1] == 27); REQUIRE(sums[1] == 27);
REQUIRE(sums[2] == 42); REQUIRE(sums[2] == 42);
*/
} }
TEST_CASE("Storing floats", "[.ClusterVector]") { TEST_CASE("Storing floats", "[.ClusterVector]") {
@ -87,10 +89,12 @@ TEST_CASE("Storing floats", "[.ClusterVector]") {
REQUIRE(cv.capacity() == 10); REQUIRE(cv.capacity() == 10);
REQUIRE(cv.size() == 2); REQUIRE(cv.size() == 2);
/*
auto sums = cv.sum(); auto sums = cv.sum();
REQUIRE(sums.size() == 2); REQUIRE(sums.size() == 2);
REQUIRE_THAT(sums[0], Catch::Matchers::WithinAbs(36.0, 1e-6)); REQUIRE_THAT(sums[0], Catch::Matchers::WithinAbs(36.0, 1e-6));
REQUIRE_THAT(sums[1], Catch::Matchers::WithinAbs(76.0, 1e-6)); REQUIRE_THAT(sums[1], Catch::Matchers::WithinAbs(76.0, 1e-6));
*/
} }
TEST_CASE("Push back more than initial capacity", "[.ClusterVector]") { TEST_CASE("Push back more than initial capacity", "[.ClusterVector]") {

View File

@ -1,6 +1,4 @@
#include "aare/Interpolator.hpp" #include "aare/Interpolator.hpp"
#include "aare/CalculateEta.hpp"
#include "aare/algorithm.hpp"
namespace aare { namespace aare {
@ -55,99 +53,4 @@ Interpolator::Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
} }
} }
// TODO: generalize to support any clustertype!!! otherwise add std::enable_if_t
// to only take Cluster2x2 and Cluster3x3
template <typename ClusterType>
std::vector<Photon>
Interpolator::interpolate(const ClusterVector<ClusterType> &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(i);
auto 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(i);
auto 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 } // namespace aare