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

bumped version

Browse Source
This commit is contained in:
froejdh_e 2024-12-16 14:24:46 +01:00
parent 29b1dc8df3
commit e6098c02ef
6 changed files with 129 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
conda-recipe/meta.yaml
View File

@ -1,6 +1,6 @@
package:
name: aare
version: 2024.11.28.dev0 #TODO! how to not duplicate this?
version: 2024.12.16.dev0 #TODO! how to not duplicate this?
source:

29
include/aare/ClusterFile.hpp
View File

@ -1,13 +1,14 @@
#pragma once
#include "aare/defs.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/NDArray.hpp"
#include "aare/defs.hpp"
#include <filesystem>
#include <fstream>
namespace aare {
struct Cluster {
struct Cluster3x3 {
int16_t x;
int16_t y;
int32_t data[9];
@ -58,21 +59,23 @@ class ClusterFile {
ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000,
const std::string &mode = "r");
~ClusterFile();
std::vector<Cluster> read_clusters(size_t n_clusters);
std::vector<Cluster> read_frame(int32_t &out_fnum);
void write_frame(int32_t frame_number, const ClusterVector<int32_t>& clusters);
std::vector<Cluster>
std::vector<Cluster3x3> read_clusters(size_t n_clusters);
std::vector<Cluster3x3> read_frame(int32_t &out_fnum);
void write_frame(int32_t frame_number,
const ClusterVector<int32_t> &clusters);
std::vector<Cluster3x3>
read_cluster_with_cut(size_t n_clusters, double *noise_map, int nx, int ny);
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(Cluster cl, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x,
double *eta3y);
size_t chunk_size() const { return m_chunk_size; }
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,
double *eta2x, double *eta2y, double *eta3x, double *eta3y);
NDArray<double, 2> calculate_eta2( ClusterVector<int>& clusters);
std::array<double,2> calculate_eta2( Cluster3x3& cl);
} // namespace aare

8
include/aare/ClusterVector.hpp
View File

@ -148,12 +148,18 @@ template <typename T, typename CoordType=int16_t> class ClusterVector {
* @brief Return a pointer to the i-th cluster
*/
std::byte *element_ptr(size_t i) { return m_data + element_offset(i); }
const std::byte * element_ptr(size_t i) const { 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; }
const std::byte *data() const { return m_data; }
std::byte const *data() const { return m_data; }
template<typename V>
V& at(size_t i) {
return *reinterpret_cast<V*>(element_ptr(i));
}
const std::string_view fmt_base() const {
//TODO! how do we match on coord_t?

2
pyproject.toml
View File

@ -4,7 +4,7 @@ build-backend = "scikit_build_core.build"
[project]
name = "aare"
version = "2024.11.28.dev0"
version = "2024.12.16.dev0"
[tool.scikit-build]
cmake.verbose = true

15
python/src/cluster_file.hpp
View File

@ -20,7 +20,7 @@ namespace py = pybind11;
using namespace ::aare;
void define_cluster_file_io_bindings(py::module &m) {
PYBIND11_NUMPY_DTYPE(Cluster, x, y, data);
PYBIND11_NUMPY_DTYPE(Cluster3x3, x, y, data);
py::class_<ClusterFile>(m, "ClusterFile")
.def(py::init<const std::filesystem::path &, size_t,
@ -29,14 +29,14 @@ void define_cluster_file_io_bindings(py::module &m) {
.def("read_clusters",
[](ClusterFile &self, size_t n_clusters) {
auto *vec =
new std::vector<Cluster>(self.read_clusters(n_clusters));
new std::vector<Cluster3x3>(self.read_clusters(n_clusters));
return return_vector(vec);
})
.def("read_frame",
[](ClusterFile &self) {
int32_t frame_number;
auto *vec =
new std::vector<Cluster>(self.read_frame(frame_number));
new std::vector<Cluster3x3>(self.read_frame(frame_number));
return py::make_tuple(frame_number, return_vector(vec));
})
.def("write_frame", &ClusterFile::write_frame)
@ -45,7 +45,7 @@ void define_cluster_file_io_bindings(py::module &m) {
py::array_t<double> noise_map, int nx, int ny) {
auto view = make_view_2d(noise_map);
auto *vec =
new std::vector<Cluster>(self.read_cluster_with_cut(
new std::vector<Cluster3x3>(self.read_cluster_with_cut(
n_clusters, view.data(), nx, ny));
return return_vector(vec);
})
@ -60,12 +60,17 @@ void define_cluster_file_io_bindings(py::module &m) {
.def("__iter__", [](ClusterFile &self) { return &self; })
.def("__next__", [](ClusterFile &self) {
auto vec =
new std::vector<Cluster>(self.read_clusters(self.chunk_size()));
new std::vector<Cluster3x3>(self.read_clusters(self.chunk_size()));
if (vec->size() == 0) {
throw py::stop_iteration();
}
return return_vector(vec);
});
m.def("calculate_eta2", []( aare::ClusterVector<int32_t> &clusters) {
auto eta2 = new NDArray<double, 2>(calculate_eta2(clusters));
return return_image_data(eta2);
});
}
#pragma GCC diagnostic pop

120
src/ClusterFile.cpp
View File

@ -1,5 +1,7 @@
#include "aare/ClusterFile.hpp"
#include <algorithm>
namespace aare {
ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size,
@ -9,17 +11,18 @@ ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size,
if (mode == "r") {
fp = fopen(fname.c_str(), "rb");
if (!fp) {
throw std::runtime_error("Could not open file for reading: " + fname.string());
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());
throw std::runtime_error("Could not open file for writing: " +
fname.string());
}
} else {
throw std::runtime_error("Unsupported mode: " + mode);
}
}
ClusterFile::~ClusterFile() { close(); }
@ -31,11 +34,13 @@ void ClusterFile::close() {
}
}
void ClusterFile::write_frame(int32_t frame_number, const ClusterVector<int32_t>& clusters){
void ClusterFile::write_frame(int32_t frame_number,
const ClusterVector<int32_t> &clusters) {
if (m_mode != "w") {
throw std::runtime_error("File not opened for writing");
}
if(!(clusters.cluster_size_x()==3) && !(clusters.cluster_size_y()==3)){
if (!(clusters.cluster_size_x() == 3) &&
!(clusters.cluster_size_y() == 3)) {
throw std::runtime_error("Only 3x3 clusters are supported");
}
fwrite(&frame_number, sizeof(frame_number), 1, fp);
@ -46,18 +51,18 @@ void ClusterFile::write_frame(int32_t frame_number, const ClusterVector<int32_t>
// fwrite(clusters.data(), sizeof(Cluster), clusters.size(), fp);
}
std::vector<Cluster> ClusterFile::read_clusters(size_t n_clusters) {
std::vector<Cluster3x3> ClusterFile::read_clusters(size_t n_clusters) {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
std::vector<Cluster> clusters(n_clusters);
std::vector<Cluster3x3> 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<Cluster *>(clusters.data());
auto buf = reinterpret_cast<Cluster3x3 *>(clusters.data());
// if there are photons left from previous frame read them first
if (nph) {
if (nph > n_clusters) {
@ -68,7 +73,7 @@ std::vector<Cluster> ClusterFile::read_clusters(size_t n_clusters) {
nn = nph;
}
nph_read += fread(reinterpret_cast<void *>(buf + nph_read),
sizeof(Cluster), nn, fp);
sizeof(Cluster3x3), nn, fp);
m_num_left = nph - nn; // write back the number of photons left
}
@ -83,7 +88,7 @@ std::vector<Cluster> ClusterFile::read_clusters(size_t n_clusters) {
nn = nph;
nph_read += fread(reinterpret_cast<void *>(buf + nph_read),
sizeof(Cluster), nn, fp);
sizeof(Cluster3x3), nn, fp);
m_num_left = nph - nn;
}
if (nph_read >= n_clusters)
@ -97,7 +102,7 @@ std::vector<Cluster> ClusterFile::read_clusters(size_t n_clusters) {
return clusters;
}
std::vector<Cluster> ClusterFile::read_frame(int32_t &out_fnum) {
std::vector<Cluster3x3> ClusterFile::read_frame(int32_t &out_fnum) {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
@ -114,22 +119,22 @@ std::vector<Cluster> ClusterFile::read_frame(int32_t &out_fnum) {
if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
throw std::runtime_error("Could not read number of clusters");
}
std::vector<Cluster> clusters(n_clusters);
std::vector<Cluster3x3> clusters(n_clusters);
if (fread(clusters.data(), sizeof(Cluster), n_clusters, fp) !=
if (fread(clusters.data(), sizeof(Cluster3x3), n_clusters, fp) !=
static_cast<size_t>(n_clusters)) {
throw std::runtime_error("Could not read clusters");
}
return clusters;
}
std::vector<Cluster> ClusterFile::read_cluster_with_cut(size_t n_clusters,
double *noise_map,
int nx, int ny) {
std::vector<Cluster3x3> 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<Cluster> clusters(n_clusters);
std::vector<Cluster3x3> 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) {
@ -143,7 +148,7 @@ std::vector<Cluster> ClusterFile::read_cluster_with_cut(size_t n_clusters,
int32_t t2max, tot1;
int32_t tot3;
// Cluster *ptr = buf;
Cluster *ptr = clusters.data();
Cluster3x3 *ptr = clusters.data();
int good = 1;
double noise;
// read photons left from previous frame
@ -161,7 +166,7 @@ std::vector<Cluster> ClusterFile::read_cluster_with_cut(size_t n_clusters,
for (size_t iph = 0; iph < nn; iph++) {
// read photons 1 by 1
size_t n_read =
fread(reinterpret_cast<void *>(ptr), sizeof(Cluster), 1, fp);
fread(reinterpret_cast<void *>(ptr), sizeof(Cluster3x3), 1, fp);
if (n_read != 1) {
clusters.resize(nph_read);
return clusters;
@ -207,7 +212,7 @@ std::vector<Cluster> ClusterFile::read_cluster_with_cut(size_t n_clusters,
for (size_t iph = 0; iph < nph; iph++) {
// // read photons 1 by 1
size_t n_read = fread(reinterpret_cast<void *>(ptr),
sizeof(Cluster), 1, fp);
sizeof(Cluster3x3), 1, fp);
if (n_read != 1) {
clusters.resize(nph_read);
return clusters;
@ -250,16 +255,78 @@ std::vector<Cluster> ClusterFile::read_cluster_with_cut(size_t n_clusters,
return clusters;
}
int ClusterFile::analyze_cluster(Cluster 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) {
NDArray<double, 2> eta2({clusters.size(), 2});
for (size_t i = 0; i < clusters.size(); i++) {
// int32_t t2;
// auto* ptr = reinterpret_cast<int32_t*> (clusters.element_ptr(i) + 2 *
// sizeof(int16_t)); analyze_cluster(clusters.at<Cluster3x3>(i), &t2,
// nullptr, nullptr, &eta2(i,0), &eta2(i,1) , nullptr, nullptr);
auto [x, y] = calculate_eta2(clusters.at<Cluster3x3>(i));
eta2(i, 0) = x;
eta2(i, 1) = y;
}
return eta2;
}
std::array<double, 2> calculate_eta2(Cluster3x3 &cl) {
std::array<double, 2> eta2{};
std::array<int32_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];
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();
switch (c) {
case cBottomLeft:
if ((cl.data[3] + cl.data[4]) != 0)
eta2[0] =
static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[1] + cl.data[4]) != 0)
eta2[1] =
static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
break;
case cBottomRight:
if ((cl.data[2] + cl.data[5]) != 0)
eta2[0] =
static_cast<double>(cl.data[5]) / (cl.data[4] + cl.data[5]);
if ((cl.data[1] + cl.data[4]) != 0)
eta2[1] =
static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
break;
case cTopLeft:
if ((cl.data[7] + cl.data[4]) != 0)
eta2[0] =
static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0)
eta2[1] =
static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
break;
case cTopRight:
if ((cl.data[5] + cl.data[4]) != 0)
eta2[0] =
static_cast<double>(cl.data[5]) / (cl.data[5] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0)
eta2[1] =
static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
break;
// default:;
}
return eta2;
}
int analyze_cluster(Cluster3x3 &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 ClusterFile::analyze_data(int32_t *data, int32_t *t2, int32_t *t3,
char *quad, double *eta2x, double *eta2y,
double *eta3x, double *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;
@ -307,6 +374,7 @@ int ClusterFile::analyze_data(int32_t *data, int32_t *t2, int32_t *t3,
if (t2)
*t2 = t2max;
}
if (t3)
*t3 = tot3;