From 0876b6891ad641d534c7b39d6d5dba3256e70af5 Mon Sep 17 00:00:00 2001 From: AliceMazzoleni99 Date: Tue, 25 Mar 2025 21:42:50 +0100 Subject: [PATCH] cpp Cluster and ClusterVector and ClusterFile are templated now, they support generic cluster types --- include/aare/Cluster.hpp | 38 +++ include/aare/ClusterFile.hpp | 499 +++++++++++++++++++++++++++++++-- include/aare/ClusterVector.hpp | 13 +- include/aare/Interpolator.hpp | 26 +- src/ClusterFile.cpp | 54 ++-- src/Interpolator.cpp | 78 +++--- 6 files changed, 619 insertions(+), 89 deletions(-) create mode 100644 include/aare/Cluster.hpp diff --git a/include/aare/Cluster.hpp b/include/aare/Cluster.hpp new file mode 100644 index 0000000..74f5281 --- /dev/null +++ b/include/aare/Cluster.hpp @@ -0,0 +1,38 @@ + +/************************************************ + * @file Cluster.hpp + * @short definition of cluster, where CoordType (x,y) give + * the cluster center coordinates and data the actual cluster data + * cluster size is given as template parameters + ***********************************************/ + +#pragma once + +#include +#include + +namespace aare { + +// requires clause c++20 maybe update +template && + std::is_integral_v>> +struct Cluster { + CoordType x; + CoordType y; + T data[ClusterSizeX * ClusterSizeY]; +}; + +// Type Traits for is_cluster_type +template +struct is_cluster : std::false_type {}; // Default case: Not a Cluster + +// TODO: Do i need the require clause here as well? +template +struct is_cluster> : std::true_type {}; // Cluster + +// helper +template constexpr bool is_cluster_v = is_cluster::value; + +} // namespace aare diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index d35f362..d61ee2b 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -1,5 +1,6 @@ #pragma once +#include "aare/Cluster.hpp" #include "aare/ClusterVector.hpp" #include "aare/NDArray.hpp" #include "aare/defs.hpp" @@ -8,14 +9,6 @@ namespace aare { -template -struct Cluster { - CoordType x; - CoordType y; - T data[ClusterSizeX * ClusterSizeY]; -}; - typedef enum { cBottomLeft = 0, cBottomRight = 1, @@ -59,6 +52,8 @@ uint32_t number_of_clusters .... */ +// TODO: change to support any type of clusters, e.g. header line with +// clsuter_size_x, cluster_size_y, /** * @brief Class to read and write cluster files * Expects data to be laid out as: @@ -71,6 +66,8 @@ uint32_t number_of_clusters * uint32_t number_of_clusters * etc. */ +template , bool>> class ClusterFile { FILE *fp{}; uint32_t m_num_left{}; @@ -97,9 +94,9 @@ class ClusterFile { * If EOF is reached the returned vector will have less than n_clusters * clusters */ - ClusterVector read_clusters(size_t n_clusters); + ClusterVector read_clusters(size_t n_clusters); - ClusterVector read_clusters(size_t n_clusters, ROI roi); + ClusterVector read_clusters(size_t n_clusters, ROI roi); /** * @brief Read a single frame from the file and return the clusters. The @@ -107,9 +104,9 @@ class ClusterFile { * @throws std::runtime_error if the file is not opened for reading or the * file pointer not at the beginning of a frame */ - ClusterVector read_frame(); + ClusterVector read_frame(); - void write_frame(const ClusterVector &clusters); + void write_frame(const ClusterVector &clusters); // Need to be migrated to support NDArray and return a ClusterVector // std::vector @@ -130,20 +127,484 @@ class ClusterFile { 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); +int analyze_cluster(Cluster &cl, int32_t *t2, int32_t *t3, + char *quad, double *eta2x, double *eta2y, double *eta3x, + double *eta3y); -template +template >> NDArray calculate_eta2(ClusterVector &clusters); -template Eta2 calculate_eta2(Cluster &cl); +// TODO: do we need rquire clauses? +template Eta2 calculate_eta2(const Cluster &cl); -Eta2 calculate_eta2(Cluster2x2 &cl); +template Eta2 calculate_eta2(const Cluster &cl); -template Eta2 calculate_eta2(ClusterType &cl); +template >> +Eta2 calculate_eta2(const ClusterType &cl); template -Eta2 calculate_eta2(Cluster &cl); +Eta2 calculate_eta2( + const Cluster &cl); + +template +ClusterFile::ClusterFile( + const std::filesystem::path &fname, size_t chunk_size, + const std::string &mode) + : m_chunk_size(chunk_size), m_mode(mode) { + + if (mode == "r") { + fp = fopen(fname.c_str(), "rb"); + if (!fp) { + 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()); + } + } else if (mode == "a") { + fp = fopen(fname.c_str(), "ab"); + if (!fp) { + throw std::runtime_error("Could not open file for appending: " + + fname.string()); + } + } else { + throw std::runtime_error("Unsupported mode: " + mode); + } +} + +template +ClusterFile::~ClusterFile() { + close(); +} + +template +void ClusterFile::close() { + if (fp) { + fclose(fp); + fp = nullptr; + } +} + +// TODO generally supported for all clsuter types +template +void ClusterFile::write_frame( + const ClusterVector &clusters) { + if (m_mode != "w" && m_mode != "a") { + throw std::runtime_error("File not opened for writing"); + } + if (!(clusters.cluster_size_x() == 3) && + !(clusters.cluster_size_y() == 3)) { + throw std::runtime_error("Only 3x3 clusters are supported"); + } + int32_t frame_number = clusters.frame_number(); + fwrite(&frame_number, sizeof(frame_number), 1, fp); + uint32_t n_clusters = clusters.size(); + fwrite(&n_clusters, sizeof(n_clusters), 1, fp); + fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp); +} + +template +ClusterVector +ClusterFile::read_clusters(size_t n_clusters) { + if (m_mode != "r") { + throw std::runtime_error("File not opened for reading"); + } + + ClusterVector 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 = clusters.data(); + // 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; + } + nph_read += fread((buf + nph_read * clusters.item_size()), + clusters.item_size(), nn, fp); + 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); + 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; +} + +template +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; + 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(); + + ClusterType 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; +} + +template +ClusterVector ClusterFile::read_frame() { + if (m_mode != "r") { + throw std::runtime_error("File not opened for reading"); + } + if (m_num_left) { + throw std::runtime_error( + "There are still photons left in the last frame"); + } + int32_t frame_number; + if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) { + throw std::runtime_error("Could not read frame number"); + } + + int32_t n_clusters; // Saved as 32bit integer in the cluster file + 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); + ClusterVector clusters(n_clusters); + clusters.set_frame_number(frame_number); + + if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) != + static_cast(n_clusters)) { + throw std::runtime_error("Could not read clusters"); + } + clusters.resize(n_clusters); + return clusters; +} + +template >> +NDArray calculate_eta2(const 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; + } + + return eta2; +} + +/** + * @brief Calculate the eta2 values for a generic sized cluster and return them + * in a Eta2 struct containing etay, etax and the index of the respective 2x2 + * subcluster. + */ +template +Eta2 calculate_eta2( + const Cluster &cl) { + Eta2 eta{}; + + // TODO loads of overhead for a 2x2 clsuter maybe keep 2x2 calculation + constexpr size_t num_2x2_subclusters = + (ClusterSizeX - 1) * (ClusterSizeY - 1); + std::array 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] = + cl.data[i * ClusterSizeX + j] + + cl.data[i * ClusterSizeX + j + 1] + + cl.data[(i + 1) * ClusterSizeX + j] + + cl.data[(i + 1) * ClusterSizeX + j + 1]; + } + + auto c = + std::max_element(sum_2x2_subcluster.begin(), sum_2x2_subcluster.end()) - + sum_2x2_subcluster.begin(); + + eta.sum = sum_2x2_subcluster[c]; + + eta.x = static_cast(cl.data[(c + 1) * ClusterSizeX + 1]) / + (cl.data[0] + cl.data[1]); + + size_t index_top_left_2x2_subcluster = + (int(c / (ClusterSizeX - 1)) + 1) * ClusterSizeX + + c % (ClusterSizeX - 1) * 2 + 1; + if ((cl.data[index_top_left_2x2_subcluster] + + cl.data[index_top_left_2x2_subcluster - 1]) != 0) + eta.x = + static_cast(cl.data[index_top_left_2x2_subcluster] / + (cl.data[index_top_left_2x2_subcluster] + + cl.data[index_top_left_2x2_subcluster - 1])); + + if ((cl.data[index_top_left_2x2_subcluster] + + cl.data[index_top_left_2x2_subcluster - ClusterSizeX]) != 0) + eta.y = static_cast( + cl.data[index_top_left_2x2_subcluster] / + (cl.data[index_top_left_2x2_subcluster] + + cl.data[index_top_left_2x2_subcluster - ClusterSizeX])); + + eta.c = c; // TODO only supported for 2x2 and 3x3 clusters -> at least no + // underyling enum class + return eta; +} + +/** + * @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2 + * struct containing etay, etax and the corner of the cluster. + */ +template Eta2 calculate_eta2(const Cluster &cl) { + Eta2 eta{}; + + 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(); + eta.sum = tot2[c]; + switch (c) { + case cBottomLeft: + if ((cl.data[3] + cl.data[4]) != 0) + eta.x = static_cast(cl.data[4]) / (cl.data[3] + cl.data[4]); + if ((cl.data[1] + cl.data[4]) != 0) + eta.y = static_cast(cl.data[4]) / (cl.data[1] + cl.data[4]); + eta.c = cBottomLeft; + break; + case cBottomRight: + if ((cl.data[2] + cl.data[5]) != 0) + eta.x = static_cast(cl.data[5]) / (cl.data[4] + cl.data[5]); + if ((cl.data[1] + cl.data[4]) != 0) + eta.y = static_cast(cl.data[4]) / (cl.data[1] + cl.data[4]); + eta.c = cBottomRight; + break; + case cTopLeft: + if ((cl.data[7] + cl.data[4]) != 0) + eta.x = static_cast(cl.data[4]) / (cl.data[3] + cl.data[4]); + if ((cl.data[7] + cl.data[4]) != 0) + eta.y = static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); + eta.c = cTopLeft; + break; + case cTopRight: + if ((cl.data[5] + cl.data[4]) != 0) + eta.x = static_cast(cl.data[5]) / (cl.data[5] + cl.data[4]); + if ((cl.data[7] + cl.data[4]) != 0) + eta.y = static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); + eta.c = cTopRight; + break; + } + return eta; +} + +template Eta2 calculate_eta2(const Cluster &cl) { + Eta2 eta{}; + + eta.x = static_cast(cl.data[1]) / (cl.data[0] + cl.data[1]); + eta.y = static_cast(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; +} + +// TODO complicated API simplify? +int analyze_cluster(Cluster &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 analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, + double *eta2x, double *eta2y, double *eta3x, double *eta3y) { + + int ok = 1; + + int32_t tot2[4]; + int32_t t2max = 0; + char c = 0; + int32_t val, tot3; + + tot3 = 0; + for (int i = 0; i < 4; i++) + tot2[i] = 0; + + for (int ix = 0; ix < 3; ix++) { + for (int iy = 0; iy < 3; iy++) { + val = data[iy * 3 + ix]; + // printf ("%d ",data[iy * 3 + ix]); + tot3 += val; + if (ix <= 1 && iy <= 1) + tot2[cBottomLeft] += val; + if (ix >= 1 && iy <= 1) + tot2[cBottomRight] += val; + if (ix <= 1 && iy >= 1) + tot2[cTopLeft] += val; + if (ix >= 1 && iy >= 1) + tot2[cTopRight] += val; + } + // printf ("\n"); + } + // printf ("\n"); + + if (t2 || quad) { + + t2max = tot2[0]; + c = cBottomLeft; + for (int i = 1; i < 4; i++) { + if (tot2[i] > t2max) { + t2max = tot2[i]; + c = i; + } + } + // printf("*** %d %d %d %d -- + // %d\n",tot2[0],tot2[1],tot2[2],tot2[3],t2max); + if (quad) + *quad = c; + if (t2) + *t2 = t2max; + } + + if (t3) + *t3 = tot3; + + if (eta2x || eta2y) { + if (eta2x) + *eta2x = 0; + if (eta2y) + *eta2y = 0; + switch (c) { + case cBottomLeft: + if (eta2x && (data[3] + data[4]) != 0) + *eta2x = static_cast(data[4]) / (data[3] + data[4]); + if (eta2y && (data[1] + data[4]) != 0) + *eta2y = static_cast(data[4]) / (data[1] + data[4]); + break; + case cBottomRight: + if (eta2x && (data[2] + data[5]) != 0) + *eta2x = static_cast(data[5]) / (data[4] + data[5]); + if (eta2y && (data[1] + data[4]) != 0) + *eta2y = static_cast(data[4]) / (data[1] + data[4]); + break; + case cTopLeft: + if (eta2x && (data[7] + data[4]) != 0) + *eta2x = static_cast(data[4]) / (data[3] + data[4]); + if (eta2y && (data[7] + data[4]) != 0) + *eta2y = static_cast(data[7]) / (data[7] + data[4]); + break; + case cTopRight: + if (eta2x && t2max != 0) + *eta2x = static_cast(data[5]) / (data[5] + data[4]); + if (eta2y && t2max != 0) + *eta2y = static_cast(data[7]) / (data[7] + data[4]); + break; + default:; + } + } + + if (eta3x || eta3y) { + if (eta3x && (data[3] + data[4] + data[5]) != 0) + *eta3x = static_cast(-data[3] + data[3 + 2]) / + (data[3] + data[4] + data[5]); + if (eta3y && (data[1] + data[4] + data[7]) != 0) + *eta3y = static_cast(-data[1] + data[2 * 3 + 1]) / + (data[1] + data[4] + data[7]); + } + + return ok; +} } // namespace aare diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp index c5e66b7..ec0fa40 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -1,4 +1,5 @@ #pragma once +#include "aare/Cluster.hpp" //TODO maybe store in seperate file !!! #include #include #include @@ -10,7 +11,9 @@ namespace aare { -template class ClusterVector; // Forward declaration +template >> +class ClusterVector; // Forward declaration /** * @brief ClusterVector is a container for clusters of various sizes. It uses a @@ -44,12 +47,10 @@ class ClusterVector> { constexpr static char m_fmt_base[] = "=h:x:\nh:y:\n({},{}){}:data:"; public: - using ClusterType = Cluster; + using ClusterType = Cluster; /** * @brief Construct a new ClusterVector object - * @param cluster_size_x size of the cluster in x direction - * @param cluster_size_y size of the cluster in y direction * @param capacity initial capacity of the buffer in number of clusters * @param frame_number frame number of the clusters. Default is 0, which is * also used to indicate that the clusters come from many frames @@ -184,6 +185,10 @@ class ClusterVector> { */ size_t size() const { return m_size; } + uint8_t cluster_size_x() const { return ClusterSizeX; } + + uint8_t cluster_size_y() const { return ClusterSizeY; } + /** * @brief Return the capacity of the buffer in number of clusters. This is * the number of clusters that can be stored in the current buffer without diff --git a/include/aare/Interpolator.hpp b/include/aare/Interpolator.hpp index 4905bce..5843046 100644 --- a/include/aare/Interpolator.hpp +++ b/include/aare/Interpolator.hpp @@ -1,29 +1,35 @@ #pragma once + +#include "aare/Cluster.hpp" +#include "aare/ClusterFile.hpp" //Cluster_3x3 +#include "aare/ClusterVector.hpp" #include "aare/NDArray.hpp" #include "aare/NDView.hpp" -#include "aare/ClusterVector.hpp" -#include "aare/ClusterFile.hpp" //Cluster_3x3 -namespace aare{ +namespace aare { -struct Photon{ +struct Photon { double x; double y; double energy; }; -class Interpolator{ +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); + public: + Interpolator(NDView etacube, NDView xbins, + NDView ybins, NDView ebins); + NDArray get_ietax() { return m_ietax; } + NDArray get_ietay() { return m_ietay; } + + template >> + std::vector interpolate(const ClusterVector &clusters); }; } // namespace aare \ No newline at end of file diff --git a/src/ClusterFile.cpp b/src/ClusterFile.cpp index c6ae470..0fc5764 100644 --- a/src/ClusterFile.cpp +++ b/src/ClusterFile.cpp @@ -4,8 +4,11 @@ namespace aare { -ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size, - const std::string &mode) +template >> +ClusterFile::ClusterFile(const std::filesystem::path &fname, + size_t chunk_size, + const std::string &mode) : m_chunk_size(chunk_size), m_mode(mode) { if (mode == "r") { @@ -31,16 +34,21 @@ ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size, } } -ClusterFile::~ClusterFile() { close(); } +template ClusterFile::~ClusterFile() { + close(); +} -void ClusterFile::close() { +template void ClusterFile::close() { if (fp) { fclose(fp); fp = nullptr; } } -void ClusterFile::write_frame(const ClusterVector &clusters) { +// TODO generally supported for all clsuter types +template +void ClusterFile::write_frame( + const ClusterVector &clusters) { if (m_mode != "w" && m_mode != "a") { throw std::runtime_error("File not opened for writing"); } @@ -55,12 +63,14 @@ void ClusterFile::write_frame(const ClusterVector &clusters) { fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp); } -ClusterVector ClusterFile::read_clusters(size_t n_clusters) { +template +ClusterVector +ClusterFile::read_clusters(size_t n_clusters) { if (m_mode != "r") { throw std::runtime_error("File not opened for reading"); } - ClusterVector clusters(3, 3, n_clusters); + ClusterVector clusters(n_clusters); int32_t iframe = 0; // frame number needs to be 4 bytes! size_t nph_read = 0; @@ -108,12 +118,14 @@ ClusterVector ClusterFile::read_clusters(size_t n_clusters) { return clusters; } -ClusterVector ClusterFile::read_clusters(size_t n_clusters, ROI roi) { +template +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); + ClusterVector clusters; clusters.reserve(n_clusters); int32_t iframe = 0; // frame number needs to be 4 bytes! @@ -124,7 +136,7 @@ ClusterVector ClusterFile::read_clusters(size_t n_clusters, ROI roi) { // auto buf = reinterpret_cast(clusters.data()); // auto buf = clusters.data(); - Cluster3x3 tmp; // this would break if the cluster size changes + ClusterType tmp; // this would break if the cluster size changes // if there are photons left from previous frame read them first if (nph) { @@ -186,7 +198,8 @@ ClusterVector ClusterFile::read_clusters(size_t n_clusters, ROI roi) { return clusters; } -ClusterVector ClusterFile::read_frame() { +template +ClusterVector ClusterFile::read_frame() { if (m_mode != "r") { throw std::runtime_error("File not opened for reading"); } @@ -204,7 +217,7 @@ ClusterVector ClusterFile::read_frame() { throw std::runtime_error("Could not read number of clusters"); } // std::vector clusters(n_clusters); - ClusterVector clusters(3, 3, n_clusters); + ClusterVector clusters(n_clusters); clusters.set_frame_number(frame_number); if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) != @@ -372,8 +385,9 @@ Eta2 calculate_eta2(Cluster &cl) { Eta2 eta{}; // TODO loads of overhead for a 2x2 clsuter maybe keep 2x2 calculation - size_t num_2x2_subclusters = (ClusterSizeX - 1) * (ClusterSizeY - 1); - std::array sum_2x2_subcluster; + constexpr size_t num_2x2_subclusters = + (ClusterSizeX - 1) * (ClusterSizeY - 1); + std::array 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] = @@ -383,9 +397,9 @@ Eta2 calculate_eta2(Cluster &cl) { cl.data[(i + 1) * ClusterSizeX + j + 1]; } - auto c = std::max_element(sum_2x2_subclusters.begin(), - sum_2x2_subcluster.end()) - - sum_2x2_subcluster.begin(); + auto c = + std::max_element(sum_2x2_subcluster.begin(), sum_2x2_subcluster.end()) - + sum_2x2_subcluster.begin(); eta.sum = sum_2x2_subcluster[c]; @@ -458,7 +472,6 @@ template Eta2 calculate_eta2(Cluster &cl) { eta.y = static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); eta.c = cTopRight; break; - // no default to allow compiler to warn about missing cases } return eta; } @@ -473,8 +486,9 @@ template Eta2 calculate_eta2(Cluster &cl) { return eta; } -int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad, - double *eta2x, double *eta2y, double *eta3x, +// TODO complicated API simplify? +int analyze_cluster(Cluster &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); diff --git a/src/Interpolator.cpp b/src/Interpolator.cpp index 85e0b5d..d95405a 100644 --- a/src/Interpolator.cpp +++ b/src/Interpolator.cpp @@ -5,7 +5,8 @@ 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) { + : 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( @@ -51,35 +52,37 @@ Interpolator::Interpolator(NDView etacube, NDView xbins, } } } - } -std::vector Interpolator::interpolate(const ClusterVector& clusters) { +template >> +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; - + // 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 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); + 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, @@ -100,44 +103,47 @@ std::vector Interpolator::interpolate(const ClusterVector& clus dY = -1.; break; case cBottomRight: - dX = 0.; + dX = 0.; dY = -1.; break; } - photon.x += m_ietax(ix, iy, 0)*2 + dX; - photon.y += m_ietay(ix, iy, 0)*2 + dY; + photon.x += m_ietax(ix, iy, 0) * 2 + dX; + photon.y += m_ietay(ix, iy, 0) * 2 + dY; photons.push_back(photon); } - }else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 2){ - for (size_t i = 0; i(i); - Eta2 eta= calculate_eta2(cluster); - + } 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); + 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 + + // 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); + auto iy = last_smaller(m_etabinsy, eta.y); - photon.x += m_ietax(ix, iy, 0)*2; //eta goes between 0 and 1 but we could move the hit anywhere in the 2x2 - photon.y += m_ietay(ix, iy, 0)*2; + photon.x += + m_ietax(ix, iy, 0) * 2; // eta goes between 0 and 1 but we could + // move the hit anywhere in the 2x2 + photon.y += m_ietay(ix, iy, 0) * 2; photons.push_back(photon); } - - }else{ - throw std::runtime_error("Only 3x3 and 2x2 clusters are supported for interpolation"); + + } else { + throw std::runtime_error( + "Only 3x3 and 2x2 clusters are supported for interpolation"); } - return photons; }