From 5614cb4673e115622267685268bbd409bb28aef2 Mon Sep 17 00:00:00 2001 From: froejdh_e Date: Wed, 5 Mar 2025 17:40:08 +0100 Subject: [PATCH] WIP --- CMakeLists.txt | 1 + include/aare/ClusterFile.hpp | 10 +++ include/aare/ClusterVector.hpp | 4 + include/aare/Interpolator.hpp | 29 ++++++ python/aare/__init__.py | 1 + python/examples/play.py | 56 +++++------- python/src/aare.code-workspace | 98 ++++++++++++++++++++ python/src/cluster.hpp | 10 ++- python/src/cluster_file.hpp | 5 ++ python/src/file.hpp | 2 + python/src/interpolation.hpp | 58 ++++++++++++ python/src/module.cpp | 2 + src/ClusterFile.cpp | 108 ++++++++++++++++++++-- src/Interpolator.cpp | 159 +++++++++++++++++++++++++++++++++ src/NDArray.test.cpp | 19 ++++ 15 files changed, 523 insertions(+), 39 deletions(-) create mode 100644 include/aare/Interpolator.hpp create mode 100644 python/src/aare.code-workspace create mode 100644 python/src/interpolation.hpp create mode 100644 src/Interpolator.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index cff4c75..bff2afe 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -346,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 diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index b796763..5bea342 100644 --- 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 read_clusters(size_t n_clusters); + ClusterVector 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 calculate_eta2(ClusterVector &clusters); Eta2 calculate_eta2(Cluster3x3 &cl); +Eta2 calculate_eta2(Cluster2x2 &cl); } // namespace aare diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp index febf06c..1c15a22 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -231,6 +231,10 @@ template class ClusterVector { return *reinterpret_cast(element_ptr(i)); } + template const V &at(size_t i) const { + return *reinterpret_cast(element_ptr(i)); + } + const std::string_view fmt_base() const { // TODO! how do we match on coord_t? return m_fmt_base; diff --git a/include/aare/Interpolator.hpp b/include/aare/Interpolator.hpp new file mode 100644 index 0000000..4905bce --- /dev/null +++ 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 m_ietax; + NDArray m_ietay; + + NDArray m_etabinsx; + NDArray m_etabinsy; + NDArray m_energy_bins; + public: + Interpolator(NDView etacube, NDView xbins, NDView ybins, NDView ebins); + NDArray get_ietax(){return m_ietax;} + NDArray get_ietay(){return m_ietay;} + + std::vector interpolate(const ClusterVector& clusters); +}; + +} // namespace aare \ No newline at end of file diff --git a/python/aare/__init__.py b/python/aare/__init__.py index f4c19cc..41deb6c 100644 --- a/python/aare/__init__.py +++ b/python/aare/__init__.py @@ -7,6 +7,7 @@ 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 diff --git a/python/examples/play.py b/python/examples/play.py index 37754df..05f3c82 100644 --- a/python/examples/play.py +++ b/python/examples/play.py @@ -1,50 +1,40 @@ 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 mpl_toolkits.axes_grid1 import make_axes_locatable import numpy as np import boost_histogram as bh import time +import tifffile -import aare - -data = np.random.normal(10, 1, 1000) - -hist = bh.Histogram(bh.axis.Regular(10, 0, 20)) -hist.fill(data) +#Directly import what we need from aare +from aare import File, ClusterFile, hitmap +from aare._aare import calculate_eta2, ClusterFinderMT, ClusterCollector -x = hist.axes[0].centers -y = hist.values() -y_err = np.sqrt(y)+1 -res = aare.fit_gaus(x, y, y_err, chi2 = True) +base = Path('/mnt/sls_det_storage/moench_data/tomcat_nanoscope_21042020/09_Moench_650um/') +# for f in base.glob('*'): +# print(f.name) - -t_elapsed = time.perf_counter()-t0 -print(f'Histogram filling took: {t_elapsed:.3f}s {total_clusters/t_elapsed/1e6:.3f}M clusters/s') +cluster_fname = base/'acq_interp_center_3.8Mfr_200V.clust' +flatfield_fname = base/'flatfield_center_200_d0_f000000000000_0.clust' -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) - -res = fit_gaus(y,x) -res2 = fit_gaus(y,x, y_err) -print(res) -print(res2) +cluster_fname.stat().st_size/1e6/4 +image = np.zeros((400,400)) +with ClusterFile(cluster_fname, chunk_size = 1000000) as f: + for clusters in f: + test = hitmap(image.shape, clusters) + break + # image += hitmap(image.shape, clusters) + # break +print('We are back in python') +# fig, ax = plt.subplots(figsize = (7,7)) +# im = ax.imshow(image) +# im.set_clim(0,1) \ No newline at end of file diff --git a/python/src/aare.code-workspace b/python/src/aare.code-workspace new file mode 100644 index 0000000..01045a6 --- /dev/null +++ b/python/src/aare.code-workspace @@ -0,0 +1,98 @@ +{ + "folders": [ + { + "path": "../../.." + }, + { + "path": "../../../../slsDetectorPackage" + } + ], + "settings": { + "files.associations": { + "compare": "cpp", + "cstdint": "cpp", + "cctype": "cpp", + "clocale": "cpp", + "cmath": "cpp", + "csignal": "cpp", + "cstdarg": "cpp", + "cstddef": "cpp", + "cstdio": "cpp", + "cstdlib": "cpp", + "cstring": "cpp", + "ctime": "cpp", + "cwchar": "cpp", + "cwctype": "cpp", + "any": "cpp", + "array": "cpp", + "atomic": "cpp", + "strstream": "cpp", + "bit": "cpp", + "*.tcc": "cpp", + "bitset": "cpp", + "cfenv": "cpp", + "charconv": "cpp", + "chrono": "cpp", + "codecvt": "cpp", + "complex": "cpp", + "concepts": "cpp", + "condition_variable": "cpp", + "deque": "cpp", + "forward_list": "cpp", + "list": "cpp", + "map": "cpp", + "set": "cpp", + "string": "cpp", + "unordered_map": "cpp", + "unordered_set": "cpp", + "vector": "cpp", + "exception": "cpp", + "algorithm": "cpp", + "functional": "cpp", + "iterator": "cpp", + "memory": "cpp", + "memory_resource": "cpp", + "numeric": "cpp", + "optional": "cpp", + "random": "cpp", + "ratio": "cpp", + "source_location": "cpp", + "string_view": "cpp", + "system_error": "cpp", + "tuple": "cpp", + "type_traits": "cpp", + "utility": "cpp", + "format": "cpp", + "fstream": "cpp", + "future": "cpp", + "initializer_list": "cpp", + "iomanip": "cpp", + "iosfwd": "cpp", + "iostream": "cpp", + "istream": "cpp", + "limits": "cpp", + "mutex": "cpp", + "new": "cpp", + "numbers": "cpp", + "ostream": "cpp", + "ranges": "cpp", + "semaphore": "cpp", + "shared_mutex": "cpp", + "span": "cpp", + "sstream": "cpp", + "stdexcept": "cpp", + "stdfloat": "cpp", + "stop_token": "cpp", + "streambuf": "cpp", + "text_encoding": "cpp", + "thread": "cpp", + "cinttypes": "cpp", + "typeindex": "cpp", + "typeinfo": "cpp", + "valarray": "cpp", + "variant": "cpp", + "regex": "cpp", + "*.ipp": "cpp" + } + } +} \ No newline at end of file diff --git a/python/src/cluster.hpp b/python/src/cluster.hpp index 792b7e6..3db816a 100644 --- a/python/src/cluster.hpp +++ b/python/src/cluster.hpp @@ -20,7 +20,13 @@ template void define_cluster_vector(py::module &m, const std::string &typestr) { auto class_name = fmt::format("ClusterVector_{}", typestr); py::class_>(m, class_name.c_str(), py::buffer_protocol()) - .def(py::init()) + .def(py::init(), + py::arg("cluster_size_x") = 3, py::arg("cluster_size_y") = 3) + .def("push_back", + [](ClusterVector &self, int x, int y, py::array_t data) { + // auto view = make_view_2d(data); + self.push_back(x, y, reinterpret_cast(data.data())); + }) .def_property_readonly("size", &ClusterVector::size) .def("item_size", &ClusterVector::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(self.sum_2x2()); return return_vector(vec); }) + .def_property_readonly("cluster_size_x", &ClusterVector::cluster_size_x) + .def_property_readonly("cluster_size_y", &ClusterVector::cluster_size_y) .def_property_readonly("capacity", &ClusterVector::capacity) .def_property("frame_number", &ClusterVector::frame_number, &ClusterVector::set_frame_number) diff --git a/python/src/cluster_file.hpp b/python/src/cluster_file.hpp index 8a431b5..f587443 100644 --- 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(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(self.read_clusters(n_clusters, roi)); + return v; + },py::return_value_policy::take_ownership) .def("read_frame", [](ClusterFile &self) { auto v = new ClusterVector(self.read_frame()); diff --git a/python/src/file.hpp b/python/src/file.hpp index c3c800c..0d64e16 100644 --- a/python/src/file.hpp +++ b/python/src/file.hpp @@ -195,6 +195,8 @@ void define_file_io_bindings(py::module &m) { py::class_(m, "ROI") .def(py::init<>()) + .def(py::init(), 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) diff --git a/python/src/interpolation.hpp b/python/src/interpolation.hpp new file mode 100644 index 0000000..02742e1 --- /dev/null +++ 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 +#include +#include +#include + +namespace py = pybind11; +void define_interpolation_bindings(py::module &m) { + + PYBIND11_NUMPY_DTYPE(aare::Photon, x,y,energy); + + py::class_(m, "Interpolator") + .def(py::init([](py::array_t etacube, py::array_t xbins, + py::array_t ybins, py::array_t 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{}; + *ptr = self.get_ietax(); + return return_image_data(ptr); + }) + .def("get_ietay", [](Interpolator& self){ + auto*ptr = new NDArray{}; + *ptr = self.get_ietay(); + return return_image_data(ptr); + }) + .def("interpolate", [](Interpolator& self, const ClusterVector& clusters){ + auto photons = self.interpolate(clusters); + auto* ptr = new std::vector{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. + )"); +} \ No newline at end of file diff --git a/python/src/module.cpp b/python/src/module.cpp index 70d143f..43f48ba 100644 --- 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 @@ -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); } \ No newline at end of file diff --git a/src/ClusterFile.cpp b/src/ClusterFile.cpp index 2928d26..37b5e89 100644 --- a/src/ClusterFile.cpp +++ b/src/ClusterFile.cpp @@ -108,6 +108,79 @@ ClusterVector ClusterFile::read_clusters(size_t n_clusters) { return clusters; } +ClusterVector ClusterFile::read_clusters(size_t n_clusters, ROI roi) { + if (m_mode != "r") { + throw std::runtime_error("File not opened for reading"); + } + + ClusterVector 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(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(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(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 ClusterFile::read_frame() { if (m_mode != "r") { throw std::runtime_error("File not opened for reading"); @@ -268,11 +341,23 @@ ClusterVector ClusterFile::read_frame() { NDArray calculate_eta2(ClusterVector &clusters) { //TOTO! make work with 2x2 clusters NDArray eta2({static_cast(clusters.size()), 2}); - for (size_t i = 0; i < clusters.size(); i++) { - auto e = calculate_eta2(clusters.at(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(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(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,19 @@ Eta2 calculate_eta2(Cluster3x3 &cl) { return eta; } + +Eta2 calculate_eta2(Cluster2x2 &cl) { + Eta2 eta{}; + + eta.x = static_cast(cl.data[0]) / (cl.data[0] + cl.data[1]); + eta.y = static_cast(cl.data[0]) / (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 + return eta; +} + + + int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad, double *eta2x, double *eta2y, double *eta3x, double *eta3y) { diff --git a/src/Interpolator.cpp b/src/Interpolator.cpp new file mode 100644 index 0000000..a4ecb94 --- /dev/null +++ b/src/Interpolator.cpp @@ -0,0 +1,159 @@ +#include "aare/Interpolator.hpp" + +namespace aare { + +Interpolator::Interpolator(NDView etacube, NDView xbins, + NDView ybins, NDView 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, can maybe be combined with a copy? + for (ssize_t k = 0; k < m_ietax.shape(2); k++) { + for (ssize_t j = 0; j < m_ietax.shape(1); j++) { + for (ssize_t i = 1; i < m_ietax.shape(0); i++) { + 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 j = 0; j < m_ietax.shape(1); j++) { + auto val = m_ietax(m_ietax.shape(0) - 1, j, k); + double norm = val < 1 ? 1 : val; + for (ssize_t i = 0; i < m_ietax.shape(0); i++) { + 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 j = 1; j < m_ietay.shape(1); j++) { + 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++) { + auto val = m_ietay(i, m_ietay.shape(1) - 1, k); + double norm = val < 1 ? 1 : val; + for (ssize_t j = 0; j < m_ietay.shape(1); j++) { + m_ietay(i, j, k) /= norm; + } + } + } + +} + +std::vector Interpolator::interpolate(const ClusterVector& clusters) { + std::vector photons; + photons.reserve(clusters.size()); + + if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) { + for (size_t i = 0; i(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 + //TODO! Could we use boost-histogram Axis.index here? + ssize_t idx = std::lower_bound(m_energy_bins.begin(), m_energy_bins.end(), photon.energy)-m_energy_bins.begin(); + auto ix = std::lower_bound(m_etabinsx.begin(), m_etabinsx.end(), eta.x)- m_etabinsx.begin(); + auto iy = std::lower_bound(m_etabinsy.begin(), m_etabinsy.end(), eta.y)- m_etabinsy.begin(); + + + // ibx=(np.abs(etabinsx - ex)).argmin() #Find out which bin the eta should land in + // iby=(np.abs(etabinsy - ey)).argmin() + 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, idx) + dX + 0.5; + photon.y += m_ietay(ix, iy, idx) + dY + 0.5; + + + // fmt::print("x: {}, y: {}, energy: {}\n", photon.x, photon.y, photon.energy); + photons.push_back(photon); + } + }else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 2){ + //TODO! Implement 2x2 interpolation + for (size_t i = 0; i(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 + //TODO! Could we use boost-histogram Axis.index here? + ssize_t idx = std::lower_bound(m_energy_bins.begin(), m_energy_bins.end(), photon.energy)-m_energy_bins.begin(); + // auto ix = std::lower_bound(m_etabinsx.begin(), m_etabinsx.end(), eta.x)- m_etabinsx.begin(); + // auto iy = std::lower_bound(m_etabinsy.begin(), m_etabinsy.end(), eta.y)- m_etabinsy.begin(); + // if(ix<0) ix=0; + // if(iy<0) iy=0; + + auto find_index = [](NDArray& etabins, double val){ + auto iter = std::min_element(etabins.begin(), etabins.end(), + [val,etabins](double a, double b) { + return std::abs(a - val) < std::abs(b - val); + }); + return std::distance(etabins.begin(), iter); + }; + auto ix = find_index(m_etabinsx, eta.x)-1; + auto iy = find_index(m_etabinsy, eta.y)-1; + + photon.x += (1-m_ietax(ix, iy, idx))*2; //eta goes between 0 and 1 but we could move the hit anywhere in the 2x2 + photon.y += (1-m_ietay(ix, iy, idx))*2; + + // photon.x = ix; + // photon.y = iy; + photons.push_back(photon); + } + + }else{ + throw std::runtime_error("Only 3x3 and 2x2 clusters are supported for interpolation"); + } + + + return photons; +} + +} // namespace aare \ No newline at end of file diff --git a/src/NDArray.test.cpp b/src/NDArray.test.cpp index 942481c..eff3e2c 100644 --- a/src/NDArray.test.cpp +++ b/src/NDArray.test.cpp @@ -2,6 +2,7 @@ #include #include #include +#include using aare::NDArray; using aare::NDView; @@ -34,6 +35,24 @@ TEST_CASE("Construct from an NDView") { } } +TEST_CASE("3D NDArray from NDView"){ + std::vector data(27); + std::iota(data.begin(), data.end(), 0); + NDView view(data.data(), Shape<3>{3, 3, 3}); + NDArray image(view); + REQUIRE(image.shape() == view.shape()); + REQUIRE(image.size() == view.size()); + REQUIRE(image.data() != view.data()); + + for(int64_t i=0; i shape{{20}}; NDArray img(shape, 3);