detectors/aare
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/aare.git synced 2025-04-21 22:30:02 +02:00
Code Issues Actions Packages Releases Activity

docs for ClusterFile

Browse Source
This commit is contained in:
froejdh_e 2025-01-15 09:15:41 +01:00
parent e1cc774d6c
commit f6d736facd
3 changed files with 197 additions and 142 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

62
include/aare/ClusterFile.hpp
View File

@ -55,6 +55,19 @@ int32_t frame_number
uint32_t number_of_clusters 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 { class ClusterFile {
FILE *fp{}; FILE *fp{};
uint32_t m_num_left{}; uint32_t m_num_left{};
@ -62,16 +75,52 @@ class ClusterFile {
const std::string m_mode; const std::string m_mode;
public: 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, ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000,
const std::string &mode = "r"); const std::string &mode = "r");
~ClusterFile();
ClusterVector<int32_t> read_clusters(size_t n_clusters);
ClusterVector<int32_t> read_frame();
void write_frame(const ClusterVector<int32_t> &clusters);
std::vector<Cluster3x3>
read_cluster_with_cut(size_t n_clusters, double *noise_map, int nx, int ny);
~ClusterFile();
/**
* @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<int32_t> 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<int32_t> read_frame();
void write_frame(const ClusterVector<int32_t> &clusters);
// Need to be migrated to support NDArray and return a ClusterVector
// std::vector<Cluster3x3>
// 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; } 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(); void close();
}; };
@ -80,7 +129,6 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
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); double *eta2x, double *eta2y, double *eta3x, double *eta3y);
NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters); NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters);
Eta2 calculate_eta2(Cluster3x3 &cl); Eta2 calculate_eta2(Cluster3x3 &cl);

18
python/src/cluster_file.hpp
View File

@ -37,15 +37,15 @@ void define_cluster_file_io_bindings(py::module &m) {
return v; return v;
}) })
.def("write_frame", &ClusterFile::write_frame) .def("write_frame", &ClusterFile::write_frame)
.def("read_cluster_with_cut", // .def("read_cluster_with_cut",
[](ClusterFile &self, size_t n_clusters, // [](ClusterFile &self, size_t n_clusters,
py::array_t<double> noise_map, int nx, int ny) { // py::array_t<double> noise_map, int nx, int ny) {
auto view = make_view_2d(noise_map); // auto view = make_view_2d(noise_map);
auto *vec = // auto *vec =
new std::vector<Cluster3x3>(self.read_cluster_with_cut( // new std::vector<Cluster3x3>(self.read_cluster_with_cut(
n_clusters, view.data(), nx, ny)); // n_clusters, view.data(), nx, ny));
return return_vector(vec); // return return_vector(vec);
}) // })
.def("__enter__", [](ClusterFile &self) { return &self; }) .def("__enter__", [](ClusterFile &self) { return &self; })
.def("__exit__", .def("__exit__",
[](ClusterFile &self, [](ClusterFile &self,

253
src/ClusterFile.cpp
View File

@ -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: " + throw std::runtime_error("Could not open file for writing: " +
fname.string()); 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 { } else {
throw std::runtime_error("Unsupported mode: " + mode); throw std::runtime_error("Unsupported mode: " + mode);
} }
@ -35,7 +41,7 @@ void ClusterFile::close() {
} }
void ClusterFile::write_frame(const ClusterVector<int32_t> &clusters) { void ClusterFile::write_frame(const ClusterVector<int32_t> &clusters) {
if (m_mode != "w") { if (m_mode != "w" && m_mode != "a") {
throw std::runtime_error("File not opened for writing"); throw std::runtime_error("File not opened for writing");
} }
if (!(clusters.cluster_size_x() == 3) && if (!(clusters.cluster_size_x() == 3) &&
@ -132,134 +138,135 @@ ClusterVector<int32_t> ClusterFile::read_frame() {
} }
std::vector<Cluster3x3> ClusterFile::read_cluster_with_cut(size_t n_clusters, // std::vector<Cluster3x3> ClusterFile::read_cluster_with_cut(size_t n_clusters,
double *noise_map, // double *noise_map,
int nx, int ny) { // int nx, int ny) {
if (m_mode != "r") { // if (m_mode != "r") {
throw std::runtime_error("File not opened for reading"); // throw std::runtime_error("File not opened for reading");
} // }
std::vector<Cluster3x3> clusters(n_clusters); // std::vector<Cluster3x3> clusters(n_clusters);
// size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf, // // size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
// uint32_t *n_left, double *noise_map, int // // uint32_t *n_left, double *noise_map, int
// nx, int ny) { // // nx, int ny) {
int iframe = 0; // int iframe = 0;
// uint32_t nph = *n_left; // // uint32_t nph = *n_left;
uint32_t nph = m_num_left; // uint32_t nph = m_num_left;
// uint32_t nn = *n_left; // // uint32_t nn = *n_left;
uint32_t nn = m_num_left; // uint32_t nn = m_num_left;
size_t nph_read = 0; // size_t nph_read = 0;
int32_t t2max, tot1; // int32_t t2max, tot1;
int32_t tot3; // int32_t tot3;
// Cluster *ptr = buf; // // Cluster *ptr = buf;
Cluster3x3 *ptr = clusters.data(); // Cluster3x3 *ptr = clusters.data();
int good = 1; // int good = 1;
double noise; // double noise;
// read photons left from previous frame // // read photons left from previous frame
if (noise_map) // if (noise_map)
printf("Using noise map\n"); // printf("Using noise map\n");
if (nph) { // if (nph) {
if (nph > n_clusters) { // if (nph > n_clusters) {
// if we have more photons left in the frame then photons to // // if we have more photons left in the frame then photons to
// read we read directly the requested number // // read we read directly the requested number
nn = n_clusters; // nn = n_clusters;
} else { // } else {
nn = nph; // nn = nph;
} // }
for (size_t iph = 0; iph < nn; iph++) { // for (size_t iph = 0; iph < nn; iph++) {
// read photons 1 by 1
size_t n_read =
fread(reinterpret_cast<void *>(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 // // read photons 1 by 1
size_t n_read = fread(reinterpret_cast<void *>(ptr), // size_t n_read =
sizeof(Cluster3x3), 1, fp); // fread(reinterpret_cast<void *>(ptr), sizeof(Cluster3x3), 1, fp);
if (n_read != 1) { // if (n_read != 1) {
clusters.resize(nph_read); // clusters.resize(nph_read);
return clusters; // return clusters;
// return nph_read; // }
} // // TODO! error handling on read
good = 1; // good = 1;
if (noise_map) { // if (noise_map) {
if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && // if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && ptr->y < ny) {
ptr->y < ny) { // tot1 = ptr->data[4];
tot1 = ptr->data[4]; // analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, NULL, NULL,
analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, // NULL);
NULL, NULL, NULL);
// noise = noise_map[ptr->y * nx + ptr->x]; // 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) {
if (tot1 > noise || t2max > 2 * noise || // ;
tot3 > 3 * noise) { // } else {
; // good = 0;
} else // printf("%d %d %f %d %d %d\n", ptr->x, ptr->y, noise,
good = 0; // tot1, t2max, tot3);
} else { // }
printf("Bad pixel number %d %d\n", ptr->x, ptr->y); // } else {
good = 0; // printf("Bad pixel number %d %d\n", ptr->x, ptr->y);
} // good = 0;
} // }
if (good) { // }
ptr++; // if (good) {
nph_read++; // ptr++;
} // nph_read++;
(m_num_left)--; // }
if (nph_read >= n_clusters) // (m_num_left)--;
break; // if (nph_read >= n_clusters)
} // break;
} // }
if (nph_read >= n_clusters) // }
break; // if (nph_read < n_clusters) {
} // // // keep on reading frames and photons until reaching
} // // n_clusters
// printf("%d\n",nph_read); // while (fread(&iframe, sizeof(iframe), 1, fp)) {
clusters.resize(nph_read); // // // printf("%d\n",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<void *>(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<double, 2> calculate_eta2(ClusterVector<int> &clusters) { NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters) {
//TOTO! make work with 2x2 clusters
NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2}); NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
for (size_t i = 0; i < clusters.size(); i++) { for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_eta2(clusters.at<Cluster3x3>(i)); auto e = calculate_eta2(clusters.at<Cluster3x3>(i));