From 04728929cba0a61831748f3179e08f9967011380 Mon Sep 17 00:00:00 2001 From: Mazzoleni Alice Francesca Date: Tue, 1 Apr 2025 18:29:08 +0200 Subject: [PATCH] implemented sum_2x2() for general clusters, only one calculate_eta2 function for all clusters --- include/aare/CalculateEta.hpp | 80 +--- include/aare/Cluster.hpp | 22 +- include/aare/ClusterFile.hpp | 8 +- include/aare/ClusterVector.hpp | 37 +- src/Cluster.test.cpp | 3 - src/ClusterFile.cpp | 704 --------------------------------- 6 files changed, 22 insertions(+), 832 deletions(-) delete mode 100644 src/ClusterFile.cpp diff --git a/include/aare/CalculateEta.hpp b/include/aare/CalculateEta.hpp index 57c1460..1a0797a 100644 --- a/include/aare/CalculateEta.hpp +++ b/include/aare/CalculateEta.hpp @@ -60,23 +60,9 @@ Eta2 calculate_eta2( const Cluster &cl) { Eta2 eta{}; - 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]; + auto max_sum = cl.max_sum_2x2(); + eta.sum = max_sum.first; + auto c = max_sum.second; size_t index_bottom_left_max_2x2_subcluster = (int(c / (ClusterSizeX - 1))) * ClusterSizeX + c % (ClusterSizeX - 1); @@ -101,66 +87,6 @@ Eta2 calculate_eta2( 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{}; - - if ((cl.data[0] + cl.data[1]) != 0) - eta.x = static_cast(cl.data[1]) / (cl.data[0] + cl.data[1]); - if ((cl.data[0] + cl.data[2]) != 0) - 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; -} - // calculates Eta3 for 3x3 cluster based on code from analyze_cluster // TODO only supported for 3x3 Clusters template Eta2 calculate_eta3(const Cluster &cl) { diff --git a/include/aare/Cluster.hpp b/include/aare/Cluster.hpp index da756dc..cc102c4 100644 --- a/include/aare/Cluster.hpp +++ b/include/aare/Cluster.hpp @@ -35,7 +35,7 @@ struct Cluster { return std::accumulate(data, data + ClusterSizeX * ClusterSizeY, 0); } - T max_sum_2x2() const { + std::pair max_sum_2x2() const { constexpr size_t num_2x2_subclusters = (ClusterSizeX - 1) * (ClusterSizeY - 1); @@ -49,8 +49,10 @@ struct Cluster { data[(i + 1) * ClusterSizeX + j + 1]; } - return *std::max_element(sum_2x2_subcluster.begin(), - sum_2x2_subcluster.end()); + int index = std::max_element(sum_2x2_subcluster.begin(), + sum_2x2_subcluster.end()) - + sum_2x2_subcluster.begin(); + return std::make_pair(sum_2x2_subcluster[index], index); } }; @@ -62,9 +64,9 @@ template struct Cluster { T sum() const { return std::accumulate(data, data + 4, 0); } - T max_sum_2x2() const { - return data[0] + data[1] + data[2] + - data[3]; // Only one possible 2x2 sum + std::pair max_sum_2x2() const { + return std::make_pair(data[0] + data[1] + data[2] + data[3], + 0); // Only one possible 2x2 sum } }; @@ -76,14 +78,16 @@ template struct Cluster { T sum() const { return std::accumulate(data, data + 9, 0); } - T max_sum_2x2() const { + std::pair max_sum_2x2() const { std::array sum_2x2_subclusters; sum_2x2_subclusters[0] = data[0] + data[1] + data[3] + data[4]; sum_2x2_subclusters[1] = data[1] + data[2] + data[4] + data[5]; sum_2x2_subclusters[2] = data[3] + data[4] + data[6] + data[7]; sum_2x2_subclusters[3] = data[4] + data[5] + data[7] + data[8]; - return *std::max_element(sum_2x2_subclusters.begin(), - sum_2x2_subclusters.end()); + int index = std::max_element(sum_2x2_subclusters.begin(), + sum_2x2_subclusters.end()) - + sum_2x2_subclusters.begin(); + return std::make_pair(sum_2x2_subclusters[index], index); } }; diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index 3fc855a..9c43326 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -415,10 +415,12 @@ bool ClusterFile::is_selected(ClusterType &cl) { return false; } } + // TODO types are wrong generalize if (m_noise_map) { - int32_t sum_1x1 = cl.data[4]; // central pixel - int32_t sum_2x2 = cl.max_sum_2x2(); // highest sum of 2x2 subclusters - int32_t sum_3x3 = cl.sum(); // sum of all pixels + int32_t sum_1x1 = cl.data[4]; // central pixel + int32_t sum_2x2 = + cl.max_sum_2x2().first; // highest sum of 2x2 subclusters + int32_t sum_3x3 = cl.sum(); // sum of all pixels auto noise = (*m_noise_map)(cl.y, cl.x); // TODO! check if this is correct diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp index 188b018..30be5eb 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -148,38 +148,6 @@ class ClusterVector> { return sums; } - /** - * @brief Return the maximum sum of the 2x2 subclusters in each cluster - * @return std::vector vector of sums for each cluster - * @throws std::runtime_error if the cluster size is not 3x3 - * @warning Only 3x3 clusters are supported for the 2x2 sum. - */ - /* only needed to calculate eta TODO: in previous PR already added calculate - sum in PR std::vector sum_2x2() { std::vector sums(m_size); const - size_t stride = item_size(); - - if (ClusterSizeX != 3 || ClusterSizeY != 3) { - throw std::runtime_error( - "Only 3x3 clusters are supported for the 2x2 sum."); - } - std::byte *ptr = m_data + 2 * sizeof(CoordType); // skip x and y - - for (size_t i = 0; i < m_size; i++) { - std::array total; - auto T_ptr = reinterpret_cast(ptr); - total[0] = T_ptr[0] + T_ptr[1] + T_ptr[3] + T_ptr[4]; - total[1] = T_ptr[1] + T_ptr[2] + T_ptr[4] + T_ptr[5]; - total[2] = T_ptr[3] + T_ptr[4] + T_ptr[6] + T_ptr[7]; - total[3] = T_ptr[4] + T_ptr[5] + T_ptr[7] + T_ptr[8]; - - sums[i] = *std::max_element(total.begin(), total.end()); - ptr += stride; - } - - return sums; - } - */ - /** * @brief Return the number of clusters in the vector */ @@ -220,9 +188,6 @@ class ClusterVector> { 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; } std::byte const *data() const { return m_data; } @@ -272,7 +237,7 @@ class ClusterVector> { m_size = new_size; } - // TODO: Generalize !!!! + // TODO: Generalize !!!! Maybe move somewhere else void apply_gain_map(const NDView gain_map) { // in principle we need to know the size of the image for this lookup // TODO! check orientations diff --git a/src/Cluster.test.cpp b/src/Cluster.test.cpp index 20c3948..7918d72 100644 --- a/src/Cluster.test.cpp +++ b/src/Cluster.test.cpp @@ -33,9 +33,6 @@ using ClusterTypes = TEST_CASE("calculate_eta2", "[.cluster][.eta_calculation]") { - // weird expect cluster_start to be in bottom_left corner -> row major -> - // check how its used should be an image!! - auto [cluster, expected_eta] = GENERATE( std::make_tuple(ClusterTypes{Cluster{0, 0, {1, 2, 3, 1}}}, Eta2{2. / 3, 3. / 4, corner::cBottomLeft, 7}), diff --git a/src/ClusterFile.cpp b/src/ClusterFile.cpp deleted file mode 100644 index 221c15d..0000000 --- a/src/ClusterFile.cpp +++ /dev/null @@ -1,704 +0,0 @@ -#include "aare/ClusterFile.hpp" - -#include - -namespace aare { - -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); - } -} - -<<<<<<< HEAD -template ClusterFile::~ClusterFile() { - close(); -} -======= -void ClusterFile::set_roi(ROI roi){ - m_roi = roi; -} - -void ClusterFile::set_noise_map(const NDView noise_map){ - m_noise_map = NDArray(noise_map); -} - -void ClusterFile::set_gain_map(const NDView gain_map){ - m_gain_map = NDArray(gain_map); -} - -ClusterFile::~ClusterFile() { close(); } ->>>>>>> developer - -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"); - } - //First write the frame number - 4 bytes - int32_t frame_number = clusters.frame_number(); - if(fwrite(&frame_number, sizeof(frame_number), 1, fp)!=1){ - throw std::runtime_error(LOCATION + "Could not write frame number"); - } - - //Then write the number of clusters - 4 bytes - uint32_t n_clusters = clusters.size(); - if(fwrite(&n_clusters, sizeof(n_clusters), 1, fp)!=1){ - throw std::runtime_error(LOCATION + "Could not write number of clusters"); - } - - //Now write the clusters in the frame - if(fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp)!=clusters.size()){ - throw std::runtime_error(LOCATION + "Could not write clusters"); - } -} - -<<<<<<< HEAD -template -ClusterVector -ClusterFile::read_clusters(size_t n_clusters) { -======= - -ClusterVector ClusterFile::read_clusters(size_t n_clusters){ - if (m_mode != "r") { - throw std::runtime_error("File not opened for reading"); - } - if (m_noise_map || m_roi){ - return read_clusters_with_cut(n_clusters); - }else{ - return read_clusters_without_cut(n_clusters); - } -} - -ClusterVector ClusterFile::read_clusters_without_cut(size_t n_clusters) { ->>>>>>> developer - 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)) { - clusters.set_frame_number(iframe); - // 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); - if(m_gain_map) - clusters.apply_gain_map(m_gain_map->view()); - return clusters; -} - -<<<<<<< HEAD -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() { -======= - - -ClusterVector ClusterFile::read_clusters_with_cut(size_t n_clusters) { - ClusterVector clusters(3,3); - clusters.reserve(n_clusters); - - // if there are photons left from previous frame read them first - if (m_num_left) { - while(m_num_left && clusters.size() < n_clusters){ - Cluster3x3 c = read_one_cluster(); - if(is_selected(c)){ - clusters.push_back(c.x, c.y, reinterpret_cast(c.data)); - } - } - } - - // we did not have enough clusters left in the previous frame - // keep on reading frames until reaching n_clusters - if (clusters.size() < n_clusters) { - // sanity check - if (m_num_left) { - throw std::runtime_error(LOCATION + "Entered second loop with clusters left\n"); - } - - int32_t frame_number = 0; // frame number needs to be 4 bytes! - while (fread(&frame_number, sizeof(frame_number), 1, fp)) { - if (fread(&m_num_left, sizeof(m_num_left), 1, fp)) { - clusters.set_frame_number(frame_number); //cluster vector will hold the last frame number - while(m_num_left && clusters.size() < n_clusters){ - Cluster3x3 c = read_one_cluster(); - if(is_selected(c)){ - clusters.push_back(c.x, c.y, reinterpret_cast(c.data)); - } - } - } - - // we have enough clusters, break out of the outer while loop - if (clusters.size() >= n_clusters) - break; - } - - } - if(m_gain_map) - clusters.apply_gain_map(m_gain_map->view()); - - return clusters; -} - -Cluster3x3 ClusterFile::read_one_cluster(){ - Cluster3x3 c; - auto rc = fread(&c, sizeof(c), 1, fp); - if (rc != 1) { - throw std::runtime_error(LOCATION + "Could not read cluster"); - } - --m_num_left; - return c; -} - -ClusterVector ClusterFile::read_frame(){ - if (m_mode != "r") { - throw std::runtime_error(LOCATION + "File not opened for reading"); - } - if (m_noise_map || m_roi){ - return read_frame_with_cut(); - }else{ - return read_frame_without_cut(); - } -} - -ClusterVector ClusterFile::read_frame_without_cut() { - 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(LOCATION + "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(LOCATION + "Could not read number of clusters"); - } - - ClusterVector clusters(3, 3, 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(LOCATION + "Could not read clusters"); - } - clusters.resize(n_clusters); - if (m_gain_map) - clusters.apply_gain_map(m_gain_map->view()); - return clusters; -} - -ClusterVector ClusterFile::read_frame_with_cut() { ->>>>>>> developer - 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"); - } - - - if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) { - throw std::runtime_error("Could not read number of clusters"); - } -<<<<<<< HEAD - // std::vector clusters(n_clusters); - ClusterVector clusters(n_clusters); -======= - - ClusterVector clusters(3, 3); - clusters.reserve(m_num_left); ->>>>>>> developer - clusters.set_frame_number(frame_number); - while(m_num_left){ - Cluster3x3 c = read_one_cluster(); - if(is_selected(c)){ - clusters.push_back(c.x, c.y, reinterpret_cast(c.data)); - } - } - if (m_gain_map) - clusters.apply_gain_map(m_gain_map->view()); - return clusters; -} - -<<<<<<< HEAD -// 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); -// // 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) { -// int iframe = 0; -// // uint32_t nph = *n_left; -// uint32_t nph = m_num_left; -// // uint32_t nn = *n_left; -// uint32_t nn = m_num_left; -// size_t nph_read = 0; -======= - ->>>>>>> developer - -bool ClusterFile::is_selected(Cluster3x3 &cl) { - //Should fail fast - if (m_roi) { - if (!(m_roi->contains(cl.x, cl.y))) { - return false; - } - } - if (m_noise_map){ - int32_t sum_1x1 = cl.data[4]; // central pixel - int32_t sum_2x2 = cl.sum_2x2(); // highest sum of 2x2 subclusters - int32_t sum_3x3 = cl.sum(); // sum of all pixels - -<<<<<<< HEAD -// 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; -// } -// for (size_t iph = 0; iph < nn; iph++) { -// // read photons 1 by 1 -// size_t n_read = -// fread(reinterpret_cast(ptr), sizeof(Cluster3x3), 1, -// fp); -// if (n_read != 1) { -// clusters.resize(nph_read); -// return clusters; -// } -// // TODO! error handling on read -// good = 1; -// if (noise_map) { -// if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && ptr->y < ny) -// { -// tot1 = ptr->data[4]; -// analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, NULL, -// NULL, -// NULL); -// noise = noise_map[ptr->y * nx + ptr->x]; -// if (tot1 > noise || t2max > 2 * noise || tot3 > 3 * -// noise) { -// ; -// } else { -// good = 0; -// printf("%d %d %f %d %d %d\n", ptr->x, ptr->y, noise, -// tot1, t2max, tot3); -// } -// } else { -// printf("Bad pixel number %d %d\n", ptr->x, ptr->y); -// good = 0; -// } -// } -// if (good) { -// ptr++; -// nph_read++; -// } -// (m_num_left)--; -// if (nph_read >= n_clusters) -// break; -// } -// } -// if (nph_read < n_clusters) { -// // // keep on reading frames and photons until reaching -// // n_clusters -// while (fread(&iframe, sizeof(iframe), 1, fp)) { -// // // printf("%d\n",nph_read); - -// if (fread(&nph, sizeof(nph), 1, fp)) { -// // // printf("** %d\n",nph); -// m_num_left = nph; -// for (size_t iph = 0; iph < nph; iph++) { -// // // read photons 1 by 1 -// size_t n_read = fread(reinterpret_cast(ptr), -// sizeof(Cluster3x3), 1, fp); -// if (n_read != 1) { -// clusters.resize(nph_read); -// return clusters; -// // return nph_read; -// } -// good = 1; -// if (noise_map) { -// if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && -// ptr->y < ny) { -// tot1 = ptr->data[4]; -// analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, -// NULL, NULL, NULL); -// // noise = noise_map[ptr->y * nx + ptr->x]; -// noise = noise_map[ptr->y + ny * ptr->x]; -// if (tot1 > noise || t2max > 2 * noise || -// tot3 > 3 * noise) { -// ; -// } else -// good = 0; -// } else { -// printf("Bad pixel number %d %d\n", ptr->x, -// ptr->y); good = 0; -// } -// } -// if (good) { -// ptr++; -// nph_read++; -// } -// (m_num_left)--; -// if (nph_read >= n_clusters) -// break; -// } -// } -// if (nph_read >= n_clusters) -// break; -// } -// } -// // printf("%d\n",nph_read); -// clusters.resize(nph_read); -// return clusters; -// } -======= - auto noise = (*m_noise_map)(cl.y, cl.x); //TODO! check if this is correct - if (sum_1x1 <= noise || sum_2x2 <= 2 * noise || sum_3x3 <= 3 * noise) { - return false; - } - } - //we passed all checks - return true; -} ->>>>>>> developer - -template -NDArray calculate_eta2(ClusterVector &clusters) { - // TOTO! make work with 2x2 clusters - NDArray eta2({static_cast(clusters.size()), 2}); -<<<<<<< HEAD - - 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"); - } - ->>>>>>> developer - 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(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(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; -} - -<<<<<<< HEAD -template Eta2 calculate_eta2(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) { -======= ->>>>>>> developer - -Eta2 calculate_eta2(Cluster2x2 &cl) { - Eta2 eta{}; - if ((cl.data[0] + cl.data[1]) != 0) - eta.x = static_cast(cl.data[1]) / (cl.data[0] + cl.data[1]); - if ((cl.data[0] + cl.data[2]) != 0) - 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; -} - - -} // namespace aare \ No newline at end of file