From e6098c02efdad6c006979c3604cc7333f903175f Mon Sep 17 00:00:00 2001 From: froejdh_e Date: Mon, 16 Dec 2024 14:24:46 +0100 Subject: [PATCH] bumped version --- conda-recipe/meta.yaml | 2 +- include/aare/ClusterFile.hpp | 29 ++++---- include/aare/ClusterVector.hpp | 8 ++- pyproject.toml | 2 +- python/src/cluster_file.hpp | 15 +++-- src/ClusterFile.cpp | 120 ++++++++++++++++++++++++++------- 6 files changed, 129 insertions(+), 47 deletions(-) diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index 86fc9a8..c3c823b 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -1,6 +1,6 @@ package: name: aare - version: 2024.11.28.dev0 #TODO! how to not duplicate this? + version: 2024.12.16.dev0 #TODO! how to not duplicate this? source: diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index edcb91e..a484560 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -1,13 +1,14 @@ #pragma once -#include "aare/defs.hpp" #include "aare/ClusterVector.hpp" +#include "aare/NDArray.hpp" +#include "aare/defs.hpp" #include #include namespace aare { -struct Cluster { +struct Cluster3x3 { int16_t x; int16_t y; int32_t data[9]; @@ -58,21 +59,23 @@ class ClusterFile { ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000, const std::string &mode = "r"); ~ClusterFile(); - std::vector read_clusters(size_t n_clusters); - std::vector read_frame(int32_t &out_fnum); - void write_frame(int32_t frame_number, const ClusterVector& clusters); - std::vector + std::vector read_clusters(size_t n_clusters); + std::vector read_frame(int32_t &out_fnum); + void write_frame(int32_t frame_number, + const ClusterVector &clusters); + std::vector read_cluster_with_cut(size_t n_clusters, double *noise_map, int nx, int ny); - int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, - double *eta2x, double *eta2y, double *eta3x, - double *eta3y); - int analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3, char *quad, - double *eta2x, double *eta2y, double *eta3x, - double *eta3y); - size_t chunk_size() const { return m_chunk_size; } void close(); }; +int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, + double *eta2x, double *eta2y, double *eta3x, double *eta3y); +int analyze_cluster(Cluster3x3& cl, int32_t *t2, int32_t *t3, char *quad, + double *eta2x, double *eta2y, double *eta3x, double *eta3y); + +NDArray calculate_eta2( ClusterVector& clusters); +std::array calculate_eta2( Cluster3x3& cl); + } // namespace aare diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp index ce8d935..98d4b37 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -148,12 +148,18 @@ template class ClusterVector { * @brief Return a pointer to the i-th cluster */ std::byte *element_ptr(size_t i) { return m_data + element_offset(i); } + const std::byte * element_ptr(size_t i) const { return m_data + element_offset(i); } size_t cluster_size_x() const { return m_cluster_size_x; } size_t cluster_size_y() const { return m_cluster_size_y; } std::byte *data() { return m_data; } - const std::byte *data() const { return m_data; } + std::byte const *data() const { return m_data; } + + template + V& at(size_t i) { + return *reinterpret_cast(element_ptr(i)); + } const std::string_view fmt_base() const { //TODO! how do we match on coord_t? diff --git a/pyproject.toml b/pyproject.toml index f194c68..b839003 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "scikit_build_core.build" [project] name = "aare" -version = "2024.11.28.dev0" +version = "2024.12.16.dev0" [tool.scikit-build] cmake.verbose = true diff --git a/python/src/cluster_file.hpp b/python/src/cluster_file.hpp index aa7fd23..82870c4 100644 --- a/python/src/cluster_file.hpp +++ b/python/src/cluster_file.hpp @@ -20,7 +20,7 @@ namespace py = pybind11; using namespace ::aare; void define_cluster_file_io_bindings(py::module &m) { - PYBIND11_NUMPY_DTYPE(Cluster, x, y, data); + PYBIND11_NUMPY_DTYPE(Cluster3x3, x, y, data); py::class_(m, "ClusterFile") .def(py::init(self.read_clusters(n_clusters)); + new std::vector(self.read_clusters(n_clusters)); return return_vector(vec); }) .def("read_frame", [](ClusterFile &self) { int32_t frame_number; auto *vec = - new std::vector(self.read_frame(frame_number)); + new std::vector(self.read_frame(frame_number)); return py::make_tuple(frame_number, return_vector(vec)); }) .def("write_frame", &ClusterFile::write_frame) @@ -45,7 +45,7 @@ void define_cluster_file_io_bindings(py::module &m) { py::array_t noise_map, int nx, int ny) { auto view = make_view_2d(noise_map); auto *vec = - new std::vector(self.read_cluster_with_cut( + new std::vector(self.read_cluster_with_cut( n_clusters, view.data(), nx, ny)); return return_vector(vec); }) @@ -60,12 +60,17 @@ void define_cluster_file_io_bindings(py::module &m) { .def("__iter__", [](ClusterFile &self) { return &self; }) .def("__next__", [](ClusterFile &self) { auto vec = - new std::vector(self.read_clusters(self.chunk_size())); + new std::vector(self.read_clusters(self.chunk_size())); if (vec->size() == 0) { throw py::stop_iteration(); } return return_vector(vec); }); + + m.def("calculate_eta2", []( aare::ClusterVector &clusters) { + auto eta2 = new NDArray(calculate_eta2(clusters)); + return return_image_data(eta2); + }); } #pragma GCC diagnostic pop \ No newline at end of file diff --git a/src/ClusterFile.cpp b/src/ClusterFile.cpp index 182726b..855e0e7 100644 --- a/src/ClusterFile.cpp +++ b/src/ClusterFile.cpp @@ -1,5 +1,7 @@ #include "aare/ClusterFile.hpp" +#include + namespace aare { ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size, @@ -9,17 +11,18 @@ ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size, if (mode == "r") { fp = fopen(fname.c_str(), "rb"); if (!fp) { - throw std::runtime_error("Could not open file for reading: " + fname.string()); + throw std::runtime_error("Could not open file for reading: " + + fname.string()); } } else if (mode == "w") { fp = fopen(fname.c_str(), "wb"); if (!fp) { - throw std::runtime_error("Could not open file for writing: " + fname.string()); + throw std::runtime_error("Could not open file for writing: " + + fname.string()); } } else { throw std::runtime_error("Unsupported mode: " + mode); } - } ClusterFile::~ClusterFile() { close(); } @@ -31,11 +34,13 @@ void ClusterFile::close() { } } -void ClusterFile::write_frame(int32_t frame_number, const ClusterVector& clusters){ +void ClusterFile::write_frame(int32_t frame_number, + const ClusterVector &clusters) { if (m_mode != "w") { throw std::runtime_error("File not opened for writing"); } - if(!(clusters.cluster_size_x()==3) && !(clusters.cluster_size_y()==3)){ + if (!(clusters.cluster_size_x() == 3) && + !(clusters.cluster_size_y() == 3)) { throw std::runtime_error("Only 3x3 clusters are supported"); } fwrite(&frame_number, sizeof(frame_number), 1, fp); @@ -46,18 +51,18 @@ void ClusterFile::write_frame(int32_t frame_number, const ClusterVector // fwrite(clusters.data(), sizeof(Cluster), clusters.size(), fp); } -std::vector ClusterFile::read_clusters(size_t n_clusters) { +std::vector ClusterFile::read_clusters(size_t n_clusters) { if (m_mode != "r") { throw std::runtime_error("File not opened for reading"); } - std::vector clusters(n_clusters); + std::vector clusters(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 = reinterpret_cast(clusters.data()); // if there are photons left from previous frame read them first if (nph) { if (nph > n_clusters) { @@ -68,7 +73,7 @@ std::vector ClusterFile::read_clusters(size_t n_clusters) { nn = nph; } nph_read += fread(reinterpret_cast(buf + nph_read), - sizeof(Cluster), nn, fp); + sizeof(Cluster3x3), nn, fp); m_num_left = nph - nn; // write back the number of photons left } @@ -83,7 +88,7 @@ std::vector ClusterFile::read_clusters(size_t n_clusters) { nn = nph; nph_read += fread(reinterpret_cast(buf + nph_read), - sizeof(Cluster), nn, fp); + sizeof(Cluster3x3), nn, fp); m_num_left = nph - nn; } if (nph_read >= n_clusters) @@ -97,7 +102,7 @@ std::vector ClusterFile::read_clusters(size_t n_clusters) { return clusters; } -std::vector ClusterFile::read_frame(int32_t &out_fnum) { +std::vector ClusterFile::read_frame(int32_t &out_fnum) { if (m_mode != "r") { throw std::runtime_error("File not opened for reading"); } @@ -114,22 +119,22 @@ std::vector ClusterFile::read_frame(int32_t &out_fnum) { if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) { throw std::runtime_error("Could not read number of clusters"); } - std::vector clusters(n_clusters); + std::vector clusters(n_clusters); - if (fread(clusters.data(), sizeof(Cluster), n_clusters, fp) != + if (fread(clusters.data(), sizeof(Cluster3x3), n_clusters, fp) != static_cast(n_clusters)) { throw std::runtime_error("Could not read clusters"); } return clusters; } -std::vector ClusterFile::read_cluster_with_cut(size_t n_clusters, - double *noise_map, - int nx, int ny) { +std::vector ClusterFile::read_cluster_with_cut(size_t n_clusters, + double *noise_map, + int nx, int ny) { if (m_mode != "r") { throw std::runtime_error("File not opened for reading"); } - std::vector clusters(n_clusters); + std::vector clusters(n_clusters); // size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf, // uint32_t *n_left, double *noise_map, int // nx, int ny) { @@ -143,7 +148,7 @@ std::vector ClusterFile::read_cluster_with_cut(size_t n_clusters, int32_t t2max, tot1; int32_t tot3; // Cluster *ptr = buf; - Cluster *ptr = clusters.data(); + Cluster3x3 *ptr = clusters.data(); int good = 1; double noise; // read photons left from previous frame @@ -161,7 +166,7 @@ std::vector ClusterFile::read_cluster_with_cut(size_t n_clusters, for (size_t iph = 0; iph < nn; iph++) { // read photons 1 by 1 size_t n_read = - fread(reinterpret_cast(ptr), sizeof(Cluster), 1, fp); + fread(reinterpret_cast(ptr), sizeof(Cluster3x3), 1, fp); if (n_read != 1) { clusters.resize(nph_read); return clusters; @@ -207,7 +212,7 @@ std::vector ClusterFile::read_cluster_with_cut(size_t n_clusters, for (size_t iph = 0; iph < nph; iph++) { // // read photons 1 by 1 size_t n_read = fread(reinterpret_cast(ptr), - sizeof(Cluster), 1, fp); + sizeof(Cluster3x3), 1, fp); if (n_read != 1) { clusters.resize(nph_read); return clusters; @@ -250,16 +255,78 @@ std::vector ClusterFile::read_cluster_with_cut(size_t n_clusters, return clusters; } -int ClusterFile::analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3, - char *quad, double *eta2x, double *eta2y, - double *eta3x, double *eta3y) { +NDArray calculate_eta2(ClusterVector &clusters) { + NDArray eta2({clusters.size(), 2}); + for (size_t i = 0; i < clusters.size(); i++) { + // int32_t t2; + // auto* ptr = reinterpret_cast (clusters.element_ptr(i) + 2 * + // sizeof(int16_t)); analyze_cluster(clusters.at(i), &t2, + // nullptr, nullptr, &eta2(i,0), &eta2(i,1) , nullptr, nullptr); + auto [x, y] = calculate_eta2(clusters.at(i)); + eta2(i, 0) = x; + eta2(i, 1) = y; + } + return eta2; +} + +std::array calculate_eta2(Cluster3x3 &cl) { + std::array eta2{}; + + std::array tot2; + tot2[0] = cl.data[0] + cl.data[1] + cl.data[3] + cl.data[4]; + tot2[1] = cl.data[1] + cl.data[2] + cl.data[4] + cl.data[5]; + tot2[2] = cl.data[3] + cl.data[4] + cl.data[6] + cl.data[7]; + 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(); + + switch (c) { + case cBottomLeft: + if ((cl.data[3] + cl.data[4]) != 0) + eta2[0] = + static_cast(cl.data[4]) / (cl.data[3] + cl.data[4]); + if ((cl.data[1] + cl.data[4]) != 0) + eta2[1] = + static_cast(cl.data[4]) / (cl.data[1] + cl.data[4]); + break; + case cBottomRight: + if ((cl.data[2] + cl.data[5]) != 0) + eta2[0] = + static_cast(cl.data[5]) / (cl.data[4] + cl.data[5]); + if ((cl.data[1] + cl.data[4]) != 0) + eta2[1] = + static_cast(cl.data[4]) / (cl.data[1] + cl.data[4]); + break; + case cTopLeft: + if ((cl.data[7] + cl.data[4]) != 0) + eta2[0] = + static_cast(cl.data[4]) / (cl.data[3] + cl.data[4]); + if ((cl.data[7] + cl.data[4]) != 0) + eta2[1] = + static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); + break; + case cTopRight: + if ((cl.data[5] + cl.data[4]) != 0) + eta2[0] = + static_cast(cl.data[5]) / (cl.data[5] + cl.data[4]); + if ((cl.data[7] + cl.data[4]) != 0) + eta2[1] = + static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); + break; + // default:; + } + return eta2; +} + +int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad, + double *eta2x, double *eta2y, double *eta3x, + double *eta3y) { return analyze_data(cl.data, t2, t3, quad, eta2x, eta2y, eta3x, eta3y); } -int ClusterFile::analyze_data(int32_t *data, int32_t *t2, int32_t *t3, - char *quad, double *eta2x, double *eta2y, - double *eta3x, double *eta3y) { +int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, + double *eta2x, double *eta2y, double *eta3x, double *eta3y) { int ok = 1; @@ -307,6 +374,7 @@ int ClusterFile::analyze_data(int32_t *data, int32_t *t2, int32_t *t3, if (t2) *t2 = t2max; } + if (t3) *t3 = tot3;