From f6d736facdacdece5c6cad3e9c3ce0eddd42dfd6 Mon Sep 17 00:00:00 2001 From: froejdh_e Date: Wed, 15 Jan 2025 09:15:41 +0100 Subject: [PATCH] docs for ClusterFile --- include/aare/ClusterFile.hpp | 66 +++++++-- python/src/cluster_file.hpp | 18 +-- src/ClusterFile.cpp | 255 ++++++++++++++++++----------------- 3 files changed, 197 insertions(+), 142 deletions(-) diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index 8a0a907..b796763 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -55,6 +55,19 @@ int32_t frame_number uint32_t number_of_clusters .... */ + +/** + * @brief Class to read and write cluster files + * Expects data to be laid out as: + * + * + * int32_t frame_number + * uint32_t number_of_clusters + * int16_t x, int16_t y, int32_t data[9] x number_of_clusters + * int32_t frame_number + * uint32_t number_of_clusters + * etc. + */ class ClusterFile { FILE *fp{}; uint32_t m_num_left{}; @@ -62,26 +75,61 @@ class ClusterFile { const std::string m_mode; public: + /** + * @brief Construct a new Cluster File object + * @param fname path to the file + * @param chunk_size number of clusters to read at a time when iterating + * over the file + * @param mode mode to open the file in. "r" for reading, "w" for writing, + * "a" for appending + * @throws std::runtime_error if the file could not be opened + */ ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000, const std::string &mode = "r"); + + ~ClusterFile(); - ClusterVector read_clusters(size_t n_clusters); - ClusterVector read_frame(); - void write_frame(const ClusterVector &clusters); - std::vector - read_cluster_with_cut(size_t n_clusters, double *noise_map, int nx, int ny); + /** + * @brief Read n_clusters clusters from the file discarding frame numbers. + * If EOF is reached the returned vector will have less than n_clusters + * clusters + */ + ClusterVector read_clusters(size_t n_clusters); + + /** + * @brief Read a single frame from the file and return the clusters. The + * cluster vector will have the frame number set. + * @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(); + + + void write_frame(const ClusterVector &clusters); + + // Need to be migrated to support NDArray and return a ClusterVector + // std::vector + // read_cluster_with_cut(size_t n_clusters, double *noise_map, int nx, int ny); + + /** + * @brief Return the chunk size + */ size_t chunk_size() const { return m_chunk_size; } + + + /** + * @brief Close the file. If not closed the file will be closed in the destructor + */ 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, +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); -Eta2 calculate_eta2( Cluster3x3& cl); +NDArray calculate_eta2(ClusterVector &clusters); +Eta2 calculate_eta2(Cluster3x3 &cl); } // namespace aare diff --git a/python/src/cluster_file.hpp b/python/src/cluster_file.hpp index baae7a1..8a431b5 100644 --- a/python/src/cluster_file.hpp +++ b/python/src/cluster_file.hpp @@ -37,15 +37,15 @@ void define_cluster_file_io_bindings(py::module &m) { return v; }) .def("write_frame", &ClusterFile::write_frame) - .def("read_cluster_with_cut", - [](ClusterFile &self, size_t n_clusters, - 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( - n_clusters, view.data(), nx, ny)); - return return_vector(vec); - }) + // .def("read_cluster_with_cut", + // [](ClusterFile &self, size_t n_clusters, + // 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( + // n_clusters, view.data(), nx, ny)); + // return return_vector(vec); + // }) .def("__enter__", [](ClusterFile &self) { return &self; }) .def("__exit__", [](ClusterFile &self, diff --git a/src/ClusterFile.cpp b/src/ClusterFile.cpp index 48f5c9a..2928d26 100644 --- a/src/ClusterFile.cpp +++ b/src/ClusterFile.cpp @@ -20,6 +20,12 @@ ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size, 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); } @@ -35,7 +41,7 @@ void ClusterFile::close() { } void ClusterFile::write_frame(const ClusterVector &clusters) { - if (m_mode != "w") { + if (m_mode != "w" && m_mode != "a") { throw std::runtime_error("File not opened for writing"); } if (!(clusters.cluster_size_x() == 3) && @@ -132,134 +138,135 @@ ClusterVector ClusterFile::read_frame() { } -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; +// 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; - int32_t t2max, tot1; - int32_t tot3; - // Cluster *ptr = buf; - Cluster3x3 *ptr = clusters.data(); - int good = 1; - double noise; - // read photons left from previous frame - if (noise_map) - printf("Using noise map\n"); +// int32_t t2max, tot1; +// int32_t tot3; +// // Cluster *ptr = buf; +// Cluster3x3 *ptr = clusters.data(); +// int good = 1; +// double noise; +// // read photons left from previous frame +// if (noise_map) +// printf("Using noise map\n"); - 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 (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; -} +// 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; +// } 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));