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

complete mess but need to install RedHat 9

Browse Source
This commit is contained in:
mazzol_a 2025-03-21 16:32:54 +01:00
parent e59a361b51
commit 6e7e81b36b
3 changed files with 162 additions and 118 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

48
include/aare/ClusterFile.hpp
View File

@ -8,16 +8,12 @@
namespace aare {
//TODO! Template this?
struct Cluster3x3 {
int16_t x;
int16_t y;
int32_t data[9];
};
struct Cluster2x2 {
int16_t x;
int16_t y;
int32_t data[4];
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = int16_t>
struct Cluster {
CoordType x;
CoordType y;
T data[ClusterSizeX * ClusterSizeY];
};
typedef enum {
@ -93,8 +89,7 @@ class ClusterFile {
*/
ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000,
const std::string &mode = "r");
~ClusterFile();
/**
@ -109,26 +104,26 @@ class ClusterFile {
/**
* @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
* @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);
// 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
* @brief Close the file. If not closed the file will be closed in the
* destructor
*/
void close();
};
@ -138,8 +133,17 @@ 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,
double *eta2x, double *eta2y, double *eta3x, double *eta3y);
NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters);
Eta2 calculate_eta2(Cluster3x3 &cl);
template <typename ClusterType>
NDArray<double, 2> calculate_eta2(ClusterVector<ClusterType> &clusters);
template <typename T> Eta2 calculate_eta2(Cluster<T, 3, 3> &cl);
Eta2 calculate_eta2(Cluster2x2 &cl);
template <typename ClusterType> Eta2 calculate_eta2(ClusterType &cl);
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
Eta2 calculate_eta2(Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl);
} // namespace aare

48
include/aare/ClusterVector.hpp
View File

@ -10,6 +10,8 @@
namespace aare {
template <typename ClusterType> class ClusterVector; // Forward declaration
/**
* @brief ClusterVector is a container for clusters of various sizes. It uses a
* contiguous memory buffer to store the clusters. It is templated on the data
@ -21,10 +23,12 @@ namespace aare {
* @tparam CoordType data type of the x and y coordinates of the cluster
* (normally int16_t)
*/
template <typename T, typename CoordType = int16_t> class ClusterVector {
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
using value_type = T;
size_t m_cluster_size_x;
size_t m_cluster_size_y;
// size_t m_cluster_size_x;
// size_t m_cluster_size_y;
std::byte *m_data{};
size_t m_size{0};
size_t m_capacity;
@ -40,6 +44,8 @@ template <typename T, typename CoordType = int16_t> class ClusterVector {
constexpr static char m_fmt_base[] = "=h:x:\nh:y:\n({},{}){}:data:";
public:
using ClusterType = Cluster<T, SizeX, SizeY>;
/**
* @brief Construct a new ClusterVector object
* @param cluster_size_x size of the cluster in x direction
@ -48,10 +54,8 @@ template <typename T, typename CoordType = int16_t> class ClusterVector {
* @param frame_number frame number of the clusters. Default is 0, which is
* also used to indicate that the clusters come from many frames
*/
ClusterVector(size_t cluster_size_x = 3, size_t cluster_size_y = 3,
size_t capacity = 1024, uint64_t frame_number = 0)
: m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y),
m_capacity(capacity), m_frame_number(frame_number) {
ClusterVector(size_t capacity = 1024, uint64_t frame_number = 0)
: m_capacity(capacity), m_frame_number(frame_number) {
allocate_buffer(capacity);
}
@ -59,10 +63,8 @@ template <typename T, typename CoordType = int16_t> class ClusterVector {
// Move constructor
ClusterVector(ClusterVector &&other) noexcept
: m_cluster_size_x(other.m_cluster_size_x),
m_cluster_size_y(other.m_cluster_size_y), m_data(other.m_data),
m_size(other.m_size), m_capacity(other.m_capacity),
m_frame_number(other.m_frame_number) {
: m_data(other.m_data), m_size(other.m_size),
m_capacity(other.m_capacity), m_frame_number(other.m_frame_number) {
other.m_data = nullptr;
other.m_size = 0;
other.m_capacity = 0;
@ -72,8 +74,6 @@ template <typename T, typename CoordType = int16_t> class ClusterVector {
ClusterVector &operator=(ClusterVector &&other) noexcept {
if (this != &other) {
delete[] m_data;
m_cluster_size_x = other.m_cluster_size_x;
m_cluster_size_y = other.m_cluster_size_y;
m_data = other.m_data;
m_size = other.m_size;
m_capacity = other.m_capacity;
@ -116,8 +116,7 @@ template <typename T, typename CoordType = int16_t> class ClusterVector {
*reinterpret_cast<CoordType *>(ptr) = y;
ptr += sizeof(CoordType);
std::copy(data, data + m_cluster_size_x * m_cluster_size_y * sizeof(T),
ptr);
std::copy(data, data + ClusterSizeX * ClusterSizeY * sizeof(T), ptr);
m_size++;
}
ClusterVector &operator+=(const ClusterVector &other) {
@ -137,7 +136,7 @@ template <typename T, typename CoordType = int16_t> class ClusterVector {
std::vector<T> sum() {
std::vector<T> sums(m_size);
const size_t stride = item_size();
const size_t n_pixels = m_cluster_size_x * m_cluster_size_y;
const size_t n_pixels = ClusterSizeX * ClusterSizeY;
std::byte *ptr = m_data + 2 * sizeof(CoordType); // skip x and y
for (size_t i = 0; i < m_size; i++) {
@ -159,7 +158,7 @@ template <typename T, typename CoordType = int16_t> class ClusterVector {
std::vector<T> sums(m_size);
const size_t stride = item_size();
if (m_cluster_size_x != 3 || m_cluster_size_y != 3) {
if (ClusterSizeX != 3 || ClusterSizeY != 3) {
throw std::runtime_error(
"Only 3x3 clusters are supported for the 2x2 sum.");
}
@ -196,8 +195,7 @@ template <typename T, typename CoordType = int16_t> class ClusterVector {
* @brief Return the size in bytes of a single cluster
*/
size_t item_size() const {
return 2 * sizeof(CoordType) +
m_cluster_size_x * m_cluster_size_y * sizeof(T);
return 2 * sizeof(CoordType) + ClusterSizeX * ClusterSizeY * sizeof(T);
}
/**
@ -217,8 +215,8 @@ template <typename T, typename CoordType = int16_t> 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; }
// 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; }
@ -227,12 +225,12 @@ template <typename T, typename CoordType = int16_t> class ClusterVector {
* @brief Return a reference to the i-th cluster casted to type V
* @tparam V type of the cluster
*/
template <typename V> V &at(size_t i) {
return *reinterpret_cast<V *>(element_ptr(i));
ClusterType &at(size_t i) {
return *reinterpret_cast<ClusterType *>(element_ptr(i));
}
template <typename V> const V &at(size_t i) const {
return *reinterpret_cast<const V *>(element_ptr(i));
const ClusterType &at(size_t i) const {
return *reinterpret_cast<const ClusterType *>(element_ptr(i));
}
const std::string_view fmt_base() const {

184
src/ClusterFile.cpp
View File

@ -59,8 +59,8 @@ ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters) {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
ClusterVector<int32_t> clusters(3,3, n_clusters);
ClusterVector<int32_t> clusters(3, 3, n_clusters);
int32_t iframe = 0; // frame number needs to be 4 bytes!
size_t nph_read = 0;
@ -78,7 +78,7 @@ ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters) {
} else {
nn = nph;
}
nph_read += fread((buf + nph_read*clusters.item_size()),
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
}
@ -93,7 +93,7 @@ ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters) {
else
nn = nph;
nph_read += fread((buf + nph_read*clusters.item_size()),
nph_read += fread((buf + nph_read * clusters.item_size()),
clusters.item_size(), nn, fp);
m_num_left = nph - nn;
}
@ -112,8 +112,8 @@ ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters, ROI roi) {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
ClusterVector<int32_t> clusters(3,3);
ClusterVector<int32_t> clusters(3, 3);
clusters.reserve(n_clusters);
int32_t iframe = 0; // frame number needs to be 4 bytes!
@ -124,7 +124,7 @@ ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters, ROI roi) {
// auto buf = reinterpret_cast<Cluster3x3 *>(clusters.data());
// auto buf = clusters.data();
Cluster3x3 tmp; //this would break if the cluster size changes
Cluster3x3 tmp; // this would break if the cluster size changes
// if there are photons left from previous frame read them first
if (nph) {
@ -135,13 +135,15 @@ ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters, ROI roi) {
} 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++){
// 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<std::byte*>(tmp.data));
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<std::byte *>(tmp.data));
nph_read++;
}
}
@ -161,10 +163,13 @@ ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters, ROI roi) {
// nph_read += fread((buf + nph_read*clusters.item_size()),
// clusters.item_size(), nn, fp);
for(size_t i = 0; i < nn; i++){
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<std::byte*>(tmp.data));
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<std::byte *>(tmp.data));
nph_read++;
}
}
@ -210,7 +215,6 @@ ClusterVector<int32_t> ClusterFile::read_frame() {
return clusters;
}
// std::vector<Cluster3x3> ClusterFile::read_cluster_with_cut(size_t n_clusters,
// double *noise_map,
// int nx, int ny) {
@ -218,7 +222,8 @@ ClusterVector<int32_t> ClusterFile::read_frame() {
// throw std::runtime_error("File not opened for reading");
// }
// 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
// // nx, int ny) {
// int iframe = 0;
@ -249,7 +254,8 @@ ClusterVector<int32_t> ClusterFile::read_frame() {
// 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);
// fread(reinterpret_cast<void *>(ptr), sizeof(Cluster3x3), 1,
// fp);
// if (n_read != 1) {
// clusters.resize(nph_read);
// return clusters;
@ -257,12 +263,15 @@ ClusterVector<int32_t> ClusterFile::read_frame() {
// // TODO! error handling on read
// good = 1;
// if (noise_map) {
// if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && ptr->y < ny) {
// 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,
// 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) {
// if (tot1 > noise || t2max > 2 * noise || tot3 > 3 *
// noise) {
// ;
// } else {
// good = 0;
@ -316,8 +325,8 @@ ClusterVector<int32_t> ClusterFile::read_frame() {
// } else
// good = 0;
// } else {
// printf("Bad pixel number %d %d\n", ptr->x, ptr->y);
// good = 0;
// printf("Bad pixel number %d %d\n", ptr->x,
// ptr->y); good = 0;
// }
// }
// if (good) {
@ -338,37 +347,81 @@ ClusterVector<int32_t> ClusterFile::read_frame() {
// return clusters;
// }
NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters) {
//TOTO! make work with 2x2 clusters
template <typename ClusterType>
NDArray<double, 2> calculate_eta2(ClusterVector<ClusterType> &clusters) {
// TOTO! make work with 2x2 clusters
NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
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<Cluster3x3>(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<Cluster2x2>(i));
eta2(i, 0) = e.x;
eta2(i, 1) = e.y;
}
}else{
throw std::runtime_error("Only 3x3 and 2x2 clusters are supported");
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_eta2<ClusterType>(clusters.at(i));
eta2(i, 0) = e.x;
eta2(i, 1) = e.y;
}
return eta2;
}
/**
* @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.
*/
Eta2 calculate_eta2(Cluster3x3 &cl) {
/**
* @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 <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
Eta2 calculate_eta2(Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
Eta2 eta{};
std::array<int32_t, 4> tot2;
// TODO loads of overhead for a 2x2 clsuter maybe keep 2x2 calculation
size_t num_2x2_subclusters = (ClusterSizeX - 1) * (ClusterSizeY - 1);
std::array<int32_t, num_2x2_subclusters> 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_subclusters.begin(),
sum_2x2_subcluster.end()) -
sum_2x2_subcluster.begin();
eta.sum = sum_2x2_subcluster[c];
eta.x = static_cast<double>(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<double>(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<double>(
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 <typename T> Eta2 calculate_eta2(Cluster<T, 3, 3> &cl) {
Eta2 eta{};
std::array<T, 4> 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];
@ -379,58 +432,47 @@ Eta2 calculate_eta2(Cluster3x3 &cl) {
switch (c) {
case cBottomLeft:
if ((cl.data[3] + cl.data[4]) != 0)
eta.x =
static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
eta.x = static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[1] + cl.data[4]) != 0)
eta.y =
static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
eta.y = static_cast<double>(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<double>(cl.data[5]) / (cl.data[4] + cl.data[5]);
eta.x = static_cast<double>(cl.data[5]) / (cl.data[4] + cl.data[5]);
if ((cl.data[1] + cl.data[4]) != 0)
eta.y =
static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
eta.y = static_cast<double>(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<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
eta.x = static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0)
eta.y =
static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
eta.y = static_cast<double>(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<double>(cl.data[5]) / (cl.data[5] + cl.data[4]);
eta.x = static_cast<double>(cl.data[5]) / (cl.data[5] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0)
eta.y =
static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
eta.y = static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
eta.c = cTopRight;
break;
// no default to allow compiler to warn about missing cases
// no default to allow compiler to warn about missing cases
}
return eta;
}
Eta2 calculate_eta2(Cluster2x2 &cl) {
template <typename T> Eta2 calculate_eta2(Cluster<T, 2, 2> &cl) {
Eta2 eta{};
eta.x = static_cast<double>(cl.data[1]) / (cl.data[0] + cl.data[1]);
eta.y = static_cast<double>(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
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;
}
int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x,
double *eta3y) {