detectors/aare
25
0
Fork 0
mirror of https://github.com/slsdetectorgroup/aare.git synced 2026-02-19 15:18:41 +01:00
Code Issues Actions Packages Releases Activity

added rounding

Browse Source
This commit is contained in:
froejdh_e
2026-01-20 15:40:26 +01:00
parent 055feb2290 f79ee55430
commit bf24122db8
183 changed files with 7154 additions and 1035 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
include/aare/ArrayExpr.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/defs.hpp"
#include <array>

455
include/aare/CalculateEta.hpp
View File

@@ -1,18 +1,13 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Cluster.hpp"
#include "aare/ClusterVector.hpp"
#include "aare/NDArray.hpp"
#include "aare/defs.hpp"
namespace aare {
enum class corner : int {
cBottomLeft = 0,
cBottomRight = 1,
cTopLeft = 2,
cTopRight = 3
};
enum class pixel : int {
pBottomLeft = 0,
pBottom = 1,
@@ -25,146 +20,428 @@ enum class pixel : int {
pTopRight = 8
};
// TODO: better to have sum after x,y
/**
* eta struct
*/
template <typename T> struct Eta2 {
double x;
double y;
int c;
T sum;
/// @brief eta in x direction
double x{};
/// @brief eta in y direction
double y{};
/// @brief index of subcluster with highest energy value (given as corner
/// relative to cluster center)
corner c{0};
/// @brief photon energy (cluster sum)
T sum{};
};
/**
* @brief Calculate the eta2 values for all clusters in a Clustervector
* @brief Calculate the eta2 values for all clusters in a ClusterVector
*/
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
NDArray<double, 2> calculate_eta2(const ClusterVector<ClusterType> &clusters) {
NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
std::vector<Eta2<typename ClusterType::value_type>>
calculate_eta2(const ClusterVector<ClusterType> &clusters) {
std::vector<Eta2<typename ClusterType::value_type>> eta2{};
eta2.reserve(clusters.size());
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_eta2(clusters[i]);
eta2(i, 0) = e.x;
eta2(i, 1) = e.y;
eta2.push_back(e);
}
return eta2;
}
/**
* @brief Calculate the full eta2 values for all clusters in a ClusterVector
*/
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
std::vector<Eta2<typename ClusterType::value_type>>
calculate_full_eta2(const ClusterVector<ClusterType> &clusters) {
std::vector<Eta2<typename ClusterType::value_type>> eta2{};
eta2.reserve(clusters.size());
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_full_eta2(clusters[i]);
eta2.push_back(e);
}
return eta2;
}
/**
* @brief Calculate eta3 for all 3x3 clusters in a ClusterVector
*/
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
std::vector<Eta2<typename ClusterType::value_type>>
calculate_eta3(const ClusterVector<ClusterType> &clusters) {
std::vector<Eta2<typename ClusterType::value_type>> eta2{};
eta2.reserve(clusters.size());
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_eta3(clusters[i]);
eta2.push_back(e);
}
return eta2;
}
/**
* @brief Calculate cross eta3 for all 3x3 clusters in a ClusterVector
*/
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
std::vector<Eta2<typename ClusterType::value_type>>
calculate_cross_eta3(const ClusterVector<ClusterType> &clusters) {
std::vector<Eta2<typename ClusterType::value_type>> eta2{};
eta2.reserve(clusters.size());
for (size_t i = 0; i < clusters.size(); i++) {
auto e = calculate_cross_eta3(clusters[i]);
eta2.push_back(e);
}
return eta2;
}
/**
* @brief helper function to calculate eta2 x and y values
* @param eta reference to the Eta2 object to update
* @param left_x value of the left pixel
* @param right_x value of the right pixel
* @param bottom_y value of the bottom pixel
* @param top_y value of the top pixel
*/
template <typename T>
inline void calculate_eta2(Eta2<T> &eta, const T left_x, const T right_x,
const T bottom_y, const T top_y) {
if ((right_x + left_x) != 0)
eta.x = static_cast<double>(right_x) /
static_cast<double>(right_x + left_x); // between (0,1) the
// closer to zero left
// value probably larger
if ((top_y + bottom_y) != 0)
eta.y = static_cast<double>(top_y) /
static_cast<double>(top_y + bottom_y); // between (0,1) the
// closer to zero bottom
// value probably larger
}
/**
* @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.
* in a Eta2 struct containing etay, etax and the index (as corner) of the
* respective 2x2 subcluster relative to the cluster center.
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
typename CoordType = uint16_t>
Eta2<T>
calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
Eta2<T> eta{};
auto max_sum = cl.max_sum_2x2();
eta.sum = max_sum.first;
auto c = max_sum.second;
static_assert(ClusterSizeX > 1 && ClusterSizeY > 1);
Eta2<T> eta{};
size_t cluster_center_index =
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
size_t index_bottom_left_max_2x2_subcluster =
(int(c / (ClusterSizeX - 1))) * ClusterSizeX + c % (ClusterSizeX - 1);
auto max_sum = cl.max_sum_2x2();
eta.sum = max_sum.sum;
corner c = max_sum.index;
// check that cluster center is in max subcluster
if (cluster_center_index != index_bottom_left_max_2x2_subcluster &&
cluster_center_index != index_bottom_left_max_2x2_subcluster + 1 &&
cluster_center_index !=
index_bottom_left_max_2x2_subcluster + ClusterSizeX &&
cluster_center_index !=
index_bottom_left_max_2x2_subcluster + ClusterSizeX + 1)
throw std::runtime_error("Photon center is not in max 2x2_subcluster");
if ((cluster_center_index - index_bottom_left_max_2x2_subcluster) %
ClusterSizeX ==
0) {
if ((cl.data[cluster_center_index + 1] +
cl.data[cluster_center_index]) != 0)
eta.x = static_cast<double>(cl.data[cluster_center_index + 1]) /
static_cast<double>((cl.data[cluster_center_index + 1] +
cl.data[cluster_center_index]));
} else {
if ((cl.data[cluster_center_index] +
cl.data[cluster_center_index - 1]) != 0)
eta.x = static_cast<double>(cl.data[cluster_center_index]) /
static_cast<double>((cl.data[cluster_center_index - 1] +
cl.data[cluster_center_index]));
// subcluster top right from center
switch (c) {
case (corner::cTopLeft):
calculate_eta2(eta, cl.data[cluster_center_index - 1],
cl.data[cluster_center_index],
cl.data[cluster_center_index - ClusterSizeX],
cl.data[cluster_center_index]);
// dx = -1
// dy = -1
break;
case (corner::cTopRight):
calculate_eta2(eta, cl.data[cluster_center_index],
cl.data[cluster_center_index + 1],
cl.data[cluster_center_index - ClusterSizeX],
cl.data[cluster_center_index]);
// dx = 0
// dy = -1
break;
case (corner::cBottomLeft):
calculate_eta2(eta, cl.data[cluster_center_index - 1],
cl.data[cluster_center_index],
cl.data[cluster_center_index],
cl.data[cluster_center_index + ClusterSizeX]);
// dx = -1
// dy = 0
break;
case (corner::cBottomRight):
calculate_eta2(eta, cl.data[cluster_center_index],
cl.data[cluster_center_index + 1],
cl.data[cluster_center_index],
cl.data[cluster_center_index + ClusterSizeX]);
// dx = 0
// dy = 0
break;
}
if ((cluster_center_index - index_bottom_left_max_2x2_subcluster) /
ClusterSizeX <
1) {
assert(cluster_center_index + ClusterSizeX <
ClusterSizeX * ClusterSizeY); // suppress warning
if ((cl.data[cluster_center_index] +
cl.data[cluster_center_index + ClusterSizeX]) != 0)
eta.y = static_cast<double>(
eta.c = c;
return eta;
}
/**
* @brief Calculate the eta2 values for a generic sized cluster and return them
* in a Eta2 struct containing etay, etax and the index (as corner) of the
* respective 2x2 subcluster relative to the cluster center.
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
Eta2<T> calculate_full_eta2(
const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
static_assert(ClusterSizeX > 1 && ClusterSizeY > 1);
Eta2<T> eta{};
constexpr size_t cluster_center_index =
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
auto max_sum = cl.max_sum_2x2();
eta.sum = max_sum.sum;
corner c = max_sum.index;
// subcluster top right from center
switch (c) {
case (corner::cTopLeft):
if (eta.sum != 0) {
eta.x = static_cast<double>(
cl.data[cluster_center_index] +
cl.data[cluster_center_index - ClusterSizeX]) /
static_cast<double>(eta.sum);
eta.y = static_cast<double>(cl.data[cluster_center_index - 1] +
cl.data[cluster_center_index]) /
static_cast<double>(eta.sum);
}
// dx = -1
// dy = -1
break;
case (corner::cTopRight):
if (eta.sum != 0) {
eta.x = static_cast<double>(
cl.data[cluster_center_index + 1] +
cl.data[cluster_center_index - ClusterSizeX + 1]) /
static_cast<double>(eta.sum);
eta.y = static_cast<double>(cl.data[cluster_center_index] +
cl.data[cluster_center_index + 1]) /
static_cast<double>(eta.sum);
}
// dx = 0
// dy = -1
break;
case (corner::cBottomLeft):
if (eta.sum != 0) {
eta.x = static_cast<double>(
cl.data[cluster_center_index] +
cl.data[cluster_center_index + ClusterSizeX]) /
static_cast<double>(
(cl.data[cluster_center_index] +
cl.data[cluster_center_index + ClusterSizeX]));
} else {
if ((cl.data[cluster_center_index] +
cl.data[cluster_center_index - ClusterSizeX]) != 0)
eta.y = static_cast<double>(cl.data[cluster_center_index]) /
static_cast<double>(
(cl.data[cluster_center_index] +
cl.data[cluster_center_index - ClusterSizeX]));
static_cast<double>(eta.sum);
eta.y = static_cast<double>(
cl.data[cluster_center_index + ClusterSizeX] +
cl.data[cluster_center_index + ClusterSizeX - 1]) /
static_cast<double>(eta.sum);
}
// dx = -1
// dy = 0
break;
case (corner::cBottomRight):
if (eta.sum != 0) {
eta.x = static_cast<double>(
cl.data[cluster_center_index + 1] +
cl.data[cluster_center_index + ClusterSizeX + 1]) /
static_cast<double>(eta.sum);
eta.y = static_cast<double>(
cl.data[cluster_center_index + ClusterSizeX] +
cl.data[cluster_center_index + ClusterSizeX + 1]) /
static_cast<double>(eta.sum);
}
// dx = 0
// dy = 0
break;
}
eta.c = c; // TODO only supported for 2x2 and 3x3 clusters -> at least no
// underyling enum class
eta.c = c;
return eta;
}
// TODO! Look up eta2 calculation - photon center should be top right corner
template <typename T>
Eta2<T> calculate_eta2(const Cluster<T, 2, 2, int16_t> &cl) {
Eta2<T> calculate_eta2(const Cluster<T, 2, 2, uint16_t> &cl) {
Eta2<T> eta{};
if ((cl.data[0] + cl.data[1]) != 0)
eta.x = static_cast<double>(cl.data[1]) / (cl.data[0] + cl.data[1]);
if ((cl.data[0] + cl.data[2]) != 0)
eta.y = static_cast<double>(cl.data[2]) / (cl.data[0] + cl.data[2]);
// TODO: maybe have as member function of cluster
const uint8_t photon_hit_index =
std::max_element(cl.data.begin(), cl.data.end()) - cl.data.begin();
eta.c = static_cast<corner>(3 - photon_hit_index);
switch (eta.c) {
case corner::cTopLeft:
calculate_eta2(eta, cl.data[2], cl.data[3], cl.data[1], cl.data[3]);
break;
case corner::cTopRight:
calculate_eta2(eta, cl.data[2], cl.data[3], cl.data[0], cl.data[2]);
break;
case corner::cBottomLeft:
calculate_eta2(eta, cl.data[0], cl.data[1], cl.data[1], cl.data[3]);
break;
case corner::cBottomRight:
calculate_eta2(eta, cl.data[0], cl.data[1], cl.data[0], cl.data[2]);
break;
}
eta.sum = cl.sum();
eta.c = static_cast<int>(corner::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 <typename T> Eta2<T> calculate_eta3(const Cluster<T, 3, 3> &cl) {
template <typename T>
Eta2<T> calculate_full_eta2(const Cluster<T, 2, 2, uint16_t> &cl) {
Eta2<T> eta{};
T sum = 0;
eta.sum = cl.sum();
std::for_each(std::begin(cl.data), std::end(cl.data),
[&sum](T x) { sum += x; });
const uint8_t photon_hit_index =
std::max_element(cl.data.begin(), cl.data.end()) - cl.data.begin();
eta.sum = sum;
eta.c = static_cast<corner>(3 - photon_hit_index);
eta.c = corner::cBottomLeft;
if (eta.sum != 0) {
eta.x = static_cast<double>(cl.data[1] + cl.data[3]) /
static_cast<double>(eta.sum);
eta.y = static_cast<double>(cl.data[2] + cl.data[3]) /
static_cast<double>(eta.sum);
}
return eta;
}
// TODO generalize
template <typename T>
Eta2<T> calculate_eta2(const Cluster<T, 1, 2, uint16_t> &cl) {
Eta2<T> eta{};
eta.x = 0;
eta.y = static_cast<double>(cl.data[1]) / cl.data[0];
eta.sum = cl.sum();
}
template <typename T>
Eta2<T> calculate_eta2(const Cluster<T, 2, 1, uint16_t> &cl) {
Eta2<T> eta{};
eta.x = static_cast<double>(cl.data[1]) / cl.data[0];
eta.y = 0;
eta.sum = cl.sum();
}
/**
* @brief calculates cross Eta3 for 3x3 cluster
* cross Eta3 calculates the eta by taking into account only the cross pixels
* {top, bottom, left, right, center}
*/
template <typename T, typename CoordType = uint16_t>
Eta2<T> calculate_cross_eta3(const Cluster<T, 3, 3, CoordType> &cl) {
Eta2<T> eta{};
T photon_energy = cl.sum();
eta.sum = photon_energy;
if ((cl.data[3] + cl.data[4] + cl.data[5]) != 0)
eta.x = static_cast<double>(-cl.data[3] + cl.data[3 + 2]) /
eta.x =
static_cast<double>(-cl.data[3] + cl.data[3 + 2]) /
(cl.data[3] + cl.data[4] + cl.data[5]);
static_cast<double>(cl.data[3] + cl.data[4] + cl.data[5]); // (-1,1)
if ((cl.data[1] + cl.data[4] + cl.data[7]) != 0)
eta.y = static_cast<double>(-cl.data[1] + cl.data[2 * 3 + 1]) /
(cl.data[1] + cl.data[4] + cl.data[7]);
static_cast<double>(cl.data[1] + cl.data[4] + cl.data[7]);
return eta;
}
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
Eta2<T> calculate_cross_eta3(
const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
static_assert(ClusterSizeX > 2 && ClusterSizeY > 2,
"calculate_eta3 only defined for clusters larger than 2x2");
if constexpr (ClusterSizeX != 3 || ClusterSizeY != 3) {
auto reduced_cluster = reduce_cluster_to_3x3(cl);
return calculate_cross_eta3(reduced_cluster);
} else {
return calculate_cross_eta3(cl);
}
}
/**
* @brief calculates Eta3 for 3x3 cluster
* It calculates the eta by taking into account all pixels in the 3x3 cluster
*/
template <typename T, typename CoordType = uint16_t>
Eta2<T> calculate_eta3(const Cluster<T, 3, 3, CoordType> &cl) {
Eta2<T> eta{};
T photon_energy = cl.sum();
eta.sum = photon_energy;
// TODO: how do we handle potential arithmetic overflows? - T could be
// uint16
if (photon_energy != 0) {
std::array<T, 2> column_sums{
static_cast<T>(cl.data[0] + cl.data[3] + cl.data[6]),
static_cast<T>(cl.data[2] + cl.data[5] + cl.data[8])};
eta.x = static_cast<double>(-column_sums[0] + column_sums[1]) /
static_cast<double>(photon_energy);
std::array<T, 2> row_sums{
static_cast<T>(cl.data[0] + cl.data[1] + cl.data[2]),
static_cast<T>(cl.data[6] + cl.data[7] + cl.data[8])};
eta.y = static_cast<double>(-row_sums[0] + row_sums[1]) /
static_cast<double>(photon_energy);
}
return eta;
}
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
Eta2<T>
calculate_eta3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
static_assert(ClusterSizeX > 2 && ClusterSizeY > 2,
"calculate_eta3 only defined for clusters larger than 2x2");
if constexpr (ClusterSizeX != 3 || ClusterSizeY != 3) {
auto reduced_cluster = reduce_cluster_to_3x3(cl);
return calculate_eta3(reduced_cluster);
} else {
return calculate_eta3(cl);
}
}
} // namespace aare

1
include/aare/CircularFifo.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <chrono>

183
include/aare/Cluster.hpp Normal file → Executable file
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
/************************************************
* @file Cluster.hpp
@@ -8,6 +9,7 @@
#pragma once
#include "defs.hpp"
#include <algorithm>
#include <array>
#include <cstdint>
@@ -17,8 +19,12 @@
namespace aare {
// requires clause c++20 maybe update
/**
* @brief Cluster struct
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = int16_t>
typename CoordType = uint16_t>
struct Cluster {
static_assert(std::is_arithmetic_v<T>, "T needs to be an arithmetic type");
@@ -27,8 +33,11 @@ struct Cluster {
static_assert(ClusterSizeX > 0 && ClusterSizeY > 0,
"Cluster sizes must be bigger than zero");
/// @brief Cluster center x coordinate (in pixel coordinates)
CoordType x;
/// @brief Cluster center y coordinate (in pixel coordinates)
CoordType y;
/// @brief Cluster data stored in row-major order starting from top-left
std::array<T, ClusterSizeX * ClusterSizeY> data;
static constexpr uint8_t cluster_size_x = ClusterSizeX;
@@ -36,9 +45,18 @@ struct Cluster {
using value_type = T;
using coord_type = CoordType;
/**
* @brief Sum of all elements in the cluster
*/
T sum() const { return std::accumulate(data.begin(), data.end(), T{}); }
std::pair<T, int> max_sum_2x2() const {
// TODO: handle 1 dimensional clusters
/**
* @brief sum of 2x2 subcluster with highest energy
* @return photon energy of subcluster, 2x2 subcluster index relative to
* cluster center
*/
Sum_index_pair<T, corner> max_sum_2x2() const {
if constexpr (cluster_size_x == 3 && cluster_size_y == 3) {
std::array<T, 4> sum_2x2_subclusters;
@@ -49,31 +67,166 @@ struct Cluster {
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);
return Sum_index_pair<T, corner>{sum_2x2_subclusters[index],
corner{index}};
} else if constexpr (cluster_size_x == 2 && cluster_size_y == 2) {
return std::make_pair(data[0] + data[1] + data[2] + data[3], 0);
return Sum_index_pair<T, corner>{
data[0] + data[1] + data[2] + data[3], corner{0}};
} else {
constexpr size_t num_2x2_subclusters =
(ClusterSizeX - 1) * (ClusterSizeY - 1);
constexpr size_t cluster_center_index =
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
std::array<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] =
data[i * ClusterSizeX + j] +
data[i * ClusterSizeX + j + 1] +
data[(i + 1) * ClusterSizeX + j] +
data[(i + 1) * ClusterSizeX + j + 1];
std::array<T, 4> sum_2x2_subcluster{0};
// subcluster top left from center
sum_2x2_subcluster[0] =
data[cluster_center_index] + data[cluster_center_index - 1] +
data[cluster_center_index - ClusterSizeX] +
data[cluster_center_index - 1 - ClusterSizeX];
// subcluster top right from center
if (ClusterSizeX > 2) {
sum_2x2_subcluster[1] =
data[cluster_center_index] +
data[cluster_center_index + 1] +
data[cluster_center_index - ClusterSizeX] +
data[cluster_center_index - ClusterSizeX + 1];
}
// subcluster bottom left from center
if (ClusterSizeY > 2) {
sum_2x2_subcluster[2] =
data[cluster_center_index] +
data[cluster_center_index - 1] +
data[cluster_center_index + ClusterSizeX] +
data[cluster_center_index + ClusterSizeX - 1];
}
// subcluster bottom right from center
if (ClusterSizeX > 2 && ClusterSizeY > 2) {
sum_2x2_subcluster[3] =
data[cluster_center_index] +
data[cluster_center_index + 1] +
data[cluster_center_index + ClusterSizeX] +
data[cluster_center_index + ClusterSizeX + 1];
}
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);
return Sum_index_pair<T, corner>{sum_2x2_subcluster[index],
corner{index}};
}
}
};
/**
* @brief Reduce a cluster to a 2x2 cluster by selecting the 2x2 block with the
* highest sum.
* @param c Cluster to reduce
* @return reduced cluster
* @note The cluster is filled using row major ordering starting at the top-left
* (thus for a max subcluster in the top left cornern the photon hit is at
* the fourth position)
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = uint16_t>
Cluster<T, 2, 2, CoordType>
reduce_to_2x2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &c) {
static_assert(ClusterSizeX >= 2 && ClusterSizeY >= 2,
"Cluster sizes must be at least 2x2 for reduction to 2x2");
Cluster<T, 2, 2, CoordType> result{};
auto [sum, index] = c.max_sum_2x2();
constexpr int16_t cluster_center_index =
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
int16_t index_top_left_max_2x2_subcluster = cluster_center_index;
switch (index) {
case corner::cTopLeft:
index_top_left_max_2x2_subcluster -= (ClusterSizeX + 1);
break;
case corner::cTopRight:
index_top_left_max_2x2_subcluster -= ClusterSizeX;
break;
case corner::cBottomLeft:
index_top_left_max_2x2_subcluster -= 1;
break;
case corner::cBottomRight:
// no change needed
break;
}
result.x = c.x;
result.y = c.y;
result.data = {
c.data[index_top_left_max_2x2_subcluster],
c.data[index_top_left_max_2x2_subcluster + 1],
c.data[index_top_left_max_2x2_subcluster + ClusterSizeX],
c.data[index_top_left_max_2x2_subcluster + ClusterSizeX + 1]};
return result;
}
template <typename T>
Cluster<T, 2, 2, uint16_t> reduce_to_2x2(const Cluster<T, 3, 3, uint16_t> &c) {
Cluster<T, 2, 2, uint16_t> result{};
auto [s, i] = c.max_sum_2x2();
result.x = c.x;
result.y = c.y;
switch (i) {
case corner::cTopLeft:
result.data = {c.data[0], c.data[1], c.data[3], c.data[4]};
break;
case corner::cTopRight:
result.data = {c.data[1], c.data[2], c.data[4], c.data[5]};
break;
case corner::cBottomLeft:
result.data = {c.data[3], c.data[4], c.data[6], c.data[7]};
break;
case corner::cBottomRight:
result.data = {c.data[4], c.data[5], c.data[7], c.data[8]};
break;
}
return result;
}
/**
* @brief Reduce a cluster to a 3x3 cluster
* @param c Cluster to reduce
* @return reduced cluster
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = int16_t>
Cluster<T, 3, 3, CoordType>
reduce_to_3x3(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &c) {
static_assert(ClusterSizeX >= 3 && ClusterSizeY >= 3,
"Cluster sizes must be at least 3x3 for reduction to 3x3");
Cluster<T, 3, 3, CoordType> result{};
int16_t cluster_center_index =
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
result.x = c.x;
result.y = c.y;
result.data = {c.data[cluster_center_index - ClusterSizeX - 1],
c.data[cluster_center_index - ClusterSizeX],
c.data[cluster_center_index - ClusterSizeX + 1],
c.data[cluster_center_index - 1],
c.data[cluster_center_index],
c.data[cluster_center_index + 1],
c.data[cluster_center_index + ClusterSizeX - 1],
c.data[cluster_center_index + ClusterSizeX],
c.data[cluster_center_index + ClusterSizeX + 1]};
return result;
}
// Type Traits for is_cluster_type
template <typename T>
struct is_cluster : std::false_type {}; // Default case: Not a Cluster

2
include/aare/ClusterCollector.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <atomic>
#include <thread>
@@ -6,6 +7,7 @@
#include "aare/ClusterVector.hpp"
#include "aare/ProducerConsumerQueue.hpp"
#include "aare/BlockingQueue.hpp"
#include "aare/defs.hpp"
namespace aare {

26
include/aare/ClusterFile.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Cluster.hpp"
@@ -189,6 +190,16 @@ class ClusterFile {
}
}
/**
* @brief Return the current position in the file (bytes)
*/
int64_t tell() {
if (!fp) {
throw std::runtime_error(LOCATION + "File not opened");
}
return ftell(fp);
}
/** @brief Open the file in specific mode
*
*/
@@ -354,15 +365,20 @@ template <typename ClusterType, typename Enable>
ClusterVector<ClusterType>
ClusterFile<ClusterType, Enable>::read_frame_without_cut() {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
throw std::runtime_error(LOCATION + "File not opened for reading");
}
if (m_num_left) {
throw std::runtime_error(
"There are still photons left in the last frame");
LOCATION + "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");
if (feof(fp))
throw std::runtime_error(LOCATION + "Unexpected end of file");
else if (ferror(fp))
throw std::runtime_error(LOCATION + "Error reading from file");
throw std::runtime_error(LOCATION + "Unexpected error (not feof or ferror) when reading frame number");
}
int32_t n_clusters; // Saved as 32bit integer in the cluster file
@@ -438,8 +454,8 @@ bool ClusterFile<ClusterType, Enable>::is_selected(ClusterType &cl) {
if (m_noise_map) {
auto sum_1x1 = cl.data[cluster_center_index]; // central pixel
auto sum_2x2 = cl.max_sum_2x2().first; // highest sum of 2x2 subclusters
auto total_sum = cl.sum(); // sum of all pixels
auto sum_2x2 = cl.max_sum_2x2().sum; // highest sum of 2x2 subclusters
auto total_sum = cl.sum(); // sum of all pixels
auto noise =
(*m_noise_map)(cl.y, cl.x); // TODO! check if this is correct

4
include/aare/ClusterFileSink.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <atomic>
#include <filesystem>
@@ -11,7 +12,8 @@
namespace aare {
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
typename = std::enable_if_t<is_cluster_v<ClusterType>>,
typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
class ClusterFileSink {
// ProducerConsumerQueue<ClusterVector<ClusterType>> *m_source;
BlockingQueue<ClusterVector<ClusterType>> *m_source;

35
include/aare/ClusterFinder.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/ClusterFile.hpp"
#include "aare/ClusterVector.hpp"
@@ -10,17 +11,25 @@
namespace aare {
template <typename ClusterType,
typename = std::enable_if_t<is_cluster_v<ClusterType>>>
struct no_2x2_cluster {
constexpr static bool value =
ClusterType::cluster_size_x > 2 && ClusterType::cluster_size_y > 2;
};
template <typename ClusterType = Cluster<int32_t, 3, 3>,
typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
class ClusterFinder {
Shape<2> m_image_size;
const PEDESTAL_TYPE m_nSigma;
const PEDESTAL_TYPE c2;
const PEDESTAL_TYPE c3;
Pedestal<PEDESTAL_TYPE> m_pedestal;
ClusterVector<ClusterType> m_clusters;
const uint32_t ClusterSizeX;
const uint32_t ClusterSizeY;
Pedestal<PEDESTAL_TYPE> m_pedestal;
ClusterVector<ClusterType> m_clusters;
static const uint8_t SavedClusterSizeX = ClusterType::cluster_size_x;
static const uint8_t SavedClusterSizeY = ClusterType::cluster_size_y;
@@ -151,11 +160,21 @@ class ClusterFinder {
for (int ic = -dx2; ic < dx2 + has_center_pixel_y; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
iy + ir >= 0 && iy + ir < frame.shape(0)) {
CT tmp = static_cast<CT>(frame(iy + ir, ix + ic)) - static_cast<CT>(m_pedestal.mean(iy + ir, ix + ic));
cluster.data[i] = tmp; // Watch for out of bounds access
}
// If the cluster type is an integral type, then we need to round the value before storing it
if constexpr(std::is_integral<CT>::value) {
auto tmp = std::round(frame(iy + ir, ix + ic) - m_pedestal.mean(iy + ir, ix + ic));
cluster.data[i] = static_cast<CT>(tmp);
}
// On the other hand if it's a floating point type we can just static cast directly
else{
auto tmp = frame(iy + ir, ix + ic) - m_pedestal.mean(iy + ir, ix + ic);
cluster.data[i] = static_cast<CT>(tmp);
}
}
i++;
}
}

4
include/aare/ClusterFinderMT.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <atomic>
#include <cstdint>
@@ -44,7 +45,8 @@ static bool is_poison(const FrameWrapper& f) {
* @tparam CT type of the cluster data
*/
template <typename ClusterType = Cluster<int32_t, 3, 3>,
typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
typename = std::enable_if_t<no_2x2_cluster<ClusterType>::value>>
class ClusterFinderMT {
protected:

57
include/aare/ClusterVector.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Cluster.hpp" //TODO maybe store in seperate file !!!
#include <algorithm>
@@ -28,12 +29,11 @@ class ClusterVector; // Forward declaration
* needed.
* @tparam T data type of the pixels in the cluster
* @tparam CoordType data type of the x and y coordinates of the cluster
* (normally int16_t)
* (normally uint16_t)
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
{
class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
std::vector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> m_data{};
int32_t m_frame_number{0}; // TODO! Check frame number size and type
@@ -87,15 +87,14 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
/**
* @brief Sum the pixels in the 2x2 subcluster with the biggest pixel sum in
* each cluster
* @return std::vector<T> vector of sums for each cluster
* @return vector of sums index pairs for each cluster
*/
std::vector<T> sum_2x2() {
std::vector<T> sums_2x2(m_data.size());
std::vector<Sum_index_pair<T, corner>> sum_2x2() {
std::vector<Sum_index_pair<T, corner>> sums_2x2(m_data.size());
std::transform(m_data.begin(), m_data.end(), sums_2x2.begin(),
[](const ClusterType &cluster) {
return cluster.max_sum_2x2().first;
});
std::transform(
m_data.begin(), m_data.end(), sums_2x2.begin(),
[](const ClusterType &cluster) { return cluster.max_sum_2x2(); });
return sums_2x2;
}
@@ -173,4 +172,42 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
}
};
/**
* @brief Reduce a cluster to a 2x2 cluster by selecting the 2x2 block with the
* highest sum.
* @param cv Clustervector containing clusters to reduce
* @return Clustervector with reduced clusters
* @note The cluster is filled using row major ordering starting at the top-left
* (thus for a max subcluster in the top left cornern the photon hit is at
* the fourth position)
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
ClusterVector<Cluster<T, 2, 2, CoordType>> reduce_to_2x2(
const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
&cv) {
ClusterVector<Cluster<T, 2, 2, CoordType>> result;
for (const auto &c : cv) {
result.push_back(reduce_to_2x2(c));
}
return result;
}
/**
* @brief Reduce a cluster to a 3x3 cluster
* @param cv Clustervector containing clusters to reduce
* @return Clustervector with reduced clusters
*/
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType>
ClusterVector<Cluster<T, 3, 3, CoordType>> reduce_to_3x3(
const ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>
&cv) {
ClusterVector<Cluster<T, 3, 3, CoordType>> result;
for (const auto &c : cv) {
result.push_back(reduce_to_3x3(c));
}
return result;
}
} // namespace aare

1
include/aare/CtbRawFile.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/FileInterface.hpp"

1
include/aare/DetectorGeometry.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/RawMasterFile.hpp" //ROI refactor away
#include "aare/defs.hpp"

1
include/aare/Dtype.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <cstdint>
#include <map>

1
include/aare/File.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/FileInterface.hpp"
#include <memory>

19
include/aare/FileInterface.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Dtype.hpp"
#include "aare/Frame.hpp"
@@ -33,15 +34,15 @@ struct FileConfig {
DetectorType detector_type{DetectorType::Unknown};
int max_frames_per_file{};
size_t total_frames{};
std::string to_string() const {
return "{ dtype: " + dtype.to_string() +
", rows: " + std::to_string(rows) +
", cols: " + std::to_string(cols) +
", geometry: " + geometry.to_string() +
", detector_type: " + ToString(detector_type) +
", max_frames_per_file: " + std::to_string(max_frames_per_file) +
", total_frames: " + std::to_string(total_frames) + " }";
}
// std::string to_string() const {
// return "{ dtype: " + dtype.to_string() +
// ", rows: " + std::to_string(rows) +
// ", cols: " + std::to_string(cols) +
// ", geometry: " + geometry.to_string() +
// ", detector_type: " + ToString(detector_type) +
// ", max_frames_per_file: " + std::to_string(max_frames_per_file) +
// ", total_frames: " + std::to_string(total_frames) + " }";
// }
};
/**

1
include/aare/FilePtr.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <cstdio>
#include <filesystem>

1
include/aare/Fit.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <cmath>

3
include/aare/Frame.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Dtype.hpp"
#include "aare/NDArray.hpp"
@@ -105,7 +106,7 @@ class Frame {
* @tparam T type of the pixels
* @return NDView<T, 2>
*/
template <typename T> NDView<T, 2> view() {
template <typename T> NDView<T, 2> view() & {
std::array<ssize_t, 2> shape = {static_cast<ssize_t>(m_rows),
static_cast<ssize_t>(m_cols)};
T *data = reinterpret_cast<T *>(m_data);

1
include/aare/GainMap.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
/************************************************
* @file GainMap.hpp
* @short function to apply gain map of image size to a vector of clusters -

262
include/aare/Interpolator.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/CalculateEta.hpp"
@@ -16,8 +17,27 @@ struct Photon {
double energy;
};
struct Coordinate2D {
double x{};
double y{};
};
class Interpolator {
/**
* @brief
* marginal CDF of eta_x (if rosenblatt applied), conditional
* CDF of eta_x conditioned on eta_y
* value at (i, j, e): F(eta_x[i] |
*eta_y[j], energy[e])
*/
NDArray<double, 3> m_ietax;
/**
* @brief
* conditional CDF of eta_y conditioned on eta_x
* value at (i,j,e): F(eta_y[j] | eta_x[i], energy[e])
*/
NDArray<double, 3> m_ietay;
NDArray<double, 1> m_etabinsx;
@@ -25,103 +45,207 @@ class Interpolator {
NDArray<double, 1> m_energy_bins;
public:
/**
* @brief Constructor for the Interpolator class
* @param etacube joint distribution of etaX, etaY and photon energy (note
* first dimension is etaX, second etaY, third photon energy)
* @param xbins bin edges for etaX
* @param ybins bin edges for etaY
* @param ebins bin edges for photon energy
*/
Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
NDView<double, 1> ybins, NDView<double, 1> ebins);
/**
* @brief Constructor for the Interpolator class
* @param xbins bin edges for etaX
* @param ybins bin edges for etaY
* @param ebins bin edges for photon energy
*/
Interpolator(NDView<double, 1> xbins, NDView<double, 1> ybins,
NDView<double, 1> ebins);
/**
* @brief transforms the joint eta distribution of etaX and etaY to the two
* independant uniform distributions based on the Roseblatt transform for
* each energy level
* @param etacube joint distribution of etaX, etaY and photon energy (first
* dimension is etaX, second etaY, third photon energy)
*/
void rosenblatttransform(NDView<double, 3> etacube);
NDArray<double, 3> get_ietax() { return m_ietax; }
NDArray<double, 3> get_ietay() { return m_ietay; }
template <typename ClusterType,
/**
* @brief interpolates the cluster centers for all clusters to a better
* precision
* @tparam ClusterType Type of Clusters to interpolate
* @tparam Etafunction Function object that calculates desired eta default:
* calculate_eta2
* @return interpolated photons (photon positions are given as double but
* following row column format e.g. x=0, y=0 means top row and first column
* of frame) (An interpolated photon position of (1.5, 2.5) corresponds to
* an estimated photon hit at the pixel center of pixel (1,2))
*/
template <auto EtaFunction = calculate_eta2, typename ClusterType,
typename Eanble = std::enable_if_t<is_cluster_v<ClusterType>>>
std::vector<Photon> interpolate(const ClusterVector<ClusterType> &clusters);
std::vector<Photon>
interpolate(const ClusterVector<ClusterType> &clusters) const;
/**
* @brief transforms the eta values to uniform coordinates based on the CDF
* ieta_x and ieta_y
* @tparam eta Eta to transform
* @return uniform coordinates {x,y}
*/
template <typename T>
Coordinate2D transform_eta_values(const Eta2<T> &eta) const;
private:
/**
* @brief bilinear interpolation of the transformed eta values
* @param ix index of etaX bin
* @param iy index of etaY bin
* @param ie index of energy bin
* @return pair of interpolated transformed eta values (ietax, ietay)
*/
template <typename T>
std::pair<double, double>
bilinear_interpolation(const size_t ix, const size_t iy, const size_t ie,
const Eta2<T> &eta) const;
};
// TODO: generalize to support any clustertype!!! otherwise add std::enable_if_t
// to only take Cluster2x2 and Cluster3x3
template <typename ClusterType, typename Enable>
template <typename T>
std::pair<double, double>
Interpolator::bilinear_interpolation(const size_t ix, const size_t iy,
const size_t ie,
const Eta2<T> &eta) const {
auto next_index_y = static_cast<ssize_t>(iy + 1) >= m_ietax.shape(1)
? m_ietax.shape(1) - 1
: iy + 1;
auto next_index_x = static_cast<ssize_t>(ix + 1) >= m_ietax.shape(0)
? m_ietax.shape(0) - 1
: ix + 1;
// bilinear interpolation
double ietax_interp_left = linear_interpolation(
{m_etabinsy(iy), m_etabinsy(iy + 1)},
{m_ietax(ix, iy, ie), m_ietax(ix, next_index_y, ie)}, eta.y);
double ietax_interp_right =
linear_interpolation({m_etabinsy(iy), m_etabinsy(iy + 1)},
{m_ietax(next_index_x, iy, ie),
m_ietax(next_index_x, next_index_y, ie)},
eta.y);
// transformed photon position x between [0,1]
double ietax_interpolated =
linear_interpolation({m_etabinsx(ix), m_etabinsx(ix + 1)},
{ietax_interp_left, ietax_interp_right}, eta.x);
double ietay_interp_left = linear_interpolation(
{m_etabinsx(ix), m_etabinsx(ix + 1)},
{m_ietay(ix, iy, ie), m_ietay(next_index_x, iy, ie)}, eta.x);
double ietay_interp_right =
linear_interpolation({m_etabinsx(ix), m_etabinsx(ix + 1)},
{m_ietay(ix, next_index_y, ie),
m_ietay(next_index_x, next_index_y, ie)},
eta.x);
// transformed photon position y between [0,1]
double ietay_interpolated =
linear_interpolation({m_etabinsy(iy), m_etabinsy(iy + 1)},
{ietay_interp_left, ietay_interp_right}, eta.y);
return {ietax_interpolated, ietay_interpolated};
}
template <typename T>
Coordinate2D Interpolator::transform_eta_values(const Eta2<T> &eta) const {
// Finding the index of the last element that is smaller
// should work fine as long as we have many bins
auto ie = last_smaller(m_energy_bins, static_cast<double>(eta.sum));
auto ix = last_smaller(m_etabinsx, eta.x);
auto iy = last_smaller(m_etabinsy, eta.y);
if (static_cast<ssize_t>(ix) >= m_etabinsx.size() - 1 ||
static_cast<ssize_t>(iy) >= m_etabinsy.size() - 1 ||
static_cast<ssize_t>(ie) >= m_energy_bins.size() - 1)
throw std::runtime_error(
fmt::format("Eta values out of bounds of eta distribution: eta.x = "
"{:.4f}, eta.y = {:.4f}, energy = {:.4f}",
eta.x, eta.y, eta.sum));
// TODO: bilinear interpolation only works if all bins have a size > 1 -
// otherwise bilinear interpolation with zero values which skew the
// results
// TODO: maybe trim the bins at the edges with zero values beforehand
// auto [ietax_interpolated, ietay_interpolated] =
// bilinear_interpolation(ix, iy, ie, eta);
return Coordinate2D{m_ietax(ix, iy, ie), m_ietay(ix, iy, ie)};
}
template <auto EtaFunction, typename ClusterType, typename Enable>
std::vector<Photon>
Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) const {
std::vector<Photon> photons;
photons.reserve(clusters.size());
if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
for (const ClusterType &cluster : clusters) {
for (const ClusterType &cluster : clusters) {
auto eta = calculate_eta2(cluster);
auto eta = EtaFunction(cluster);
Photon photon;
photon.x = cluster.x;
photon.y = cluster.y;
photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
Photon photon;
photon.x = cluster.x;
photon.y = cluster.y;
photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
// auto ie = nearest_index(m_energy_bins, photon.energy)-1;
// auto ix = nearest_index(m_etabinsx, eta.x)-1;
// auto iy = nearest_index(m_etabinsy, eta.y)-1;
// Finding the index of the last element that is smaller
// should work fine as long as we have many bins
auto ie = last_smaller(m_energy_bins, photon.energy);
auto ix = last_smaller(m_etabinsx, eta.x);
auto iy = last_smaller(m_etabinsy, eta.y);
auto uniform_coordinates = transform_eta_values(eta);
// fmt::print("ex: {}, ix: {}, iy: {}\n", ie, ix, iy);
if (EtaFunction == &calculate_eta2<typename ClusterType::value_type,
ClusterType::cluster_size_x,
ClusterType::cluster_size_y,
typename ClusterType::coord_type> ||
EtaFunction ==
&calculate_full_eta2<typename ClusterType::value_type,
ClusterType::cluster_size_x,
ClusterType::cluster_size_y,
typename ClusterType::coord_type>) {
double dX{}, dY{};
double dX, dY;
// cBottomLeft = 0,
// cBottomRight = 1,
// cTopLeft = 2,
// cTopRight = 3
switch (static_cast<corner>(eta.c)) {
// TODO: could also chaneg the sign of the eta calculation
switch (eta.c) {
case corner::cTopLeft:
dX = -1.;
dY = 0;
dX = -1.0;
dY = -1.0;
break;
case corner::cTopRight:;
dX = 0;
dY = 0;
dX = 0.0;
dY = -1.0;
break;
case corner::cBottomLeft:
dX = -1.;
dY = -1.;
dX = -1.0;
dY = 0.0;
break;
case corner::cBottomRight:
dX = 0.;
dY = -1.;
dX = 0.0;
dY = 0.0;
break;
}
photon.x += m_ietax(ix, iy, ie) * 2 + dX;
photon.y += m_ietay(ix, iy, ie) * 2 + dY;
photons.push_back(photon);
}
} else if (clusters.cluster_size_x() == 2 ||
clusters.cluster_size_y() == 2) {
for (const ClusterType &cluster : clusters) {
auto eta = calculate_eta2(cluster);
Photon photon;
photon.x = cluster.x;
photon.y = cluster.y;
photon.energy = static_cast<decltype(photon.energy)>(eta.sum);
// Now do some actual interpolation.
// Find which energy bin the cluster is in
// auto ie = nearest_index(m_energy_bins, photon.energy)-1;
// auto ix = nearest_index(m_etabinsx, eta.x)-1;
// auto iy = nearest_index(m_etabinsy, eta.y)-1;
// Finding the index of the last element that is smaller
// should work fine as long as we have many bins
auto ie = last_smaller(m_energy_bins, photon.energy);
auto ix = last_smaller(m_etabinsx, eta.x);
auto iy = last_smaller(m_etabinsy, eta.y);
photon.x += m_ietax(ix, iy, ie) *
2; // eta goes between 0 and 1 but we could move the hit
// anywhere in the 2x2
photon.y += m_ietay(ix, iy, ie) * 2;
photons.push_back(photon);
photon.x = photon.x + 0.5 + uniform_coordinates.x +
dX; // use pixel center + 0.5
photon.y =
photon.y + 0.5 + uniform_coordinates.y +
dY; // eta2 calculates the ratio between bottom and sum of
// bottom and top shift by 1 add eta value correctly
} else {
photon.x += uniform_coordinates.x;
photon.y += uniform_coordinates.y;
}
} else {
throw std::runtime_error(
"Only 3x3 and 2x2 clusters are supported for interpolation");
photons.push_back(photon);
}
return photons;

1
include/aare/JungfrauDataFile.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <cstdint>
#include <filesystem>

687
include/aare/NDArray.hpp
View File

@@ -1,12 +1,10 @@
// SPDX-License-Identifier: MPL-2.0
//
// Container holding image data, or a time series of image data in contigious
// memory. Used for all data processing in Aare.
//
#pragma once
/*
Container holding image data, or a time series of image data in contigious
memory.
TODO! Add expression templates for operators
*/
#include "aare/ArrayExpr.hpp"
#include "aare/NDView.hpp"
@@ -25,12 +23,17 @@ template <typename T, ssize_t Ndim = 2>
class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
std::array<ssize_t, Ndim> shape_;
std::array<ssize_t, Ndim> strides_;
size_t size_{}; //TODO! do we need to store size when we have shape?
size_t size_{}; // TODO! do we need to store size when we have shape?
T *data_;
public:
///////////////////////////////////////////////////////////////////////////////
// Constructors
//
///////////////////////////////////////////////////////////////////////////////
/**
* @brief Default constructor. Will construct an empty NDArray.
* @brief Default constructor. Constructs an empty NDArray.
*
*/
NDArray() : shape_(), strides_(c_strides<Ndim>(shape_)), data_(nullptr) {};
@@ -43,8 +46,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
*/
explicit NDArray(std::array<ssize_t, Ndim> shape)
: shape_(shape), strides_(c_strides<Ndim>(shape_)),
size_(num_elements(shape_)),
data_(new T[size_]) {}
size_(num_elements(shape_)), data_(new T[size_]) {}
/**
* @brief Construct a new NDArray object with a shape and value.
@@ -56,6 +58,10 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
this->operator=(value);
}
// Allow NDArray of different type and dimension to be friend classes
// This is needed for the move constructor from NDArray<T,Ndim+1>
template <typename U, ssize_t Dim> friend class NDArray;
/**
* @brief Construct a new NDArray object from a NDView.
* @note The data is copied from the view to the NDArray.
@@ -66,44 +72,67 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
std::copy(v.begin(), v.end(), begin());
}
/**
* @brief Construct a new NDArray object from an std::array.
*/
template <size_t Size>
NDArray(const std::array<T, Size> &arr) : NDArray<T, 1>({Size}) {
std::copy(arr.begin(), arr.end(), begin());
}
// Move constructor
/**
* @brief Move construct a new NDArray object. Cheap since it just
* reassigns the pointer and copy size/strides.
*
* @param other
*/
NDArray(NDArray &&other) noexcept
: shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
size_(other.size_), data_(other.data_) {
other.reset(); // TODO! is this necessary?
other.reset(); // Needed to avoid double free
}
//Move constructor from an an array with Ndim + 1
/**
* @brief Move construct a new NDArray object from an array with Ndim + 1.
* Can be used to drop a dimension cheaply.
* @param other
*/
template <ssize_t M, typename = std::enable_if_t<(M == Ndim + 1)>>
NDArray(NDArray<T, M> &&other)
NDArray(NDArray<T, M> &&other)
: shape_(drop_first_dim(other.shape())),
strides_(c_strides<Ndim>(shape_)), size_(num_elements(shape_)),
data_(other.data()) {
// For now only allow move if the size matches, to avoid unreachable data
// if the use case arises we can remove this check
if(size() != other.size()) {
data_ = nullptr; // avoid double free, other will clean up the memory in it's destructor
throw std::runtime_error(LOCATION +
"Size mismatch in move constructor of NDArray<T, Ndim-1>");
}
// For now only allow move if the size matches, to avoid unreachable
// data if the use case arises we can remove this check
if (size() != other.size()) {
data_ = nullptr; // avoid double free, other will clean up the
// memory in it's destructor
throw std::runtime_error(
LOCATION +
"Size mismatch in move constructor of NDArray<T, Ndim-1>");
}
other.reset();
}
// Copy constructor
/**
* @brief Copy construct a new NDArray object from another NDArray.
*
* @param other
*/
NDArray(const NDArray &other)
: shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
size_(other.size_), data_(new T[size_]) {
std::copy(other.data_, other.data_ + size_, data_);
}
// Conversion operator from array expression to array
/**
* @brief Conversion from a ArrayExpr to an actual NDArray. Used when
* the expression is evaluated and data needed.
*
* @tparam E
* @param expr
*/
template <typename E>
NDArray(ArrayExpr<E, Ndim> &&expr) : NDArray(expr.shape()) {
for (size_t i = 0; i < size_; ++i) {
@@ -111,23 +140,129 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
}
}
/**
* @brief Destroy the NDArray object. Frees the allocated memory.
*
*/
~NDArray() { delete[] data_; }
auto begin() { return data_; }
auto end() { return data_ + size_; }
///////////////////////////////////////////////////////////////////////////////
// Iterators and indexing
//
///////////////////////////////////////////////////////////////////////////////
auto begin() const { return data_; }
auto end() const { return data_ + size_; }
auto *begin() { return data_; }
const auto *begin() const { return data_; }
auto *end() { return data_ + size_; }
const auto *end() const { return data_ + size_; }
/*
* @brief Access element at given multi-dimensional index.
* i.e. arr(i,j,k,...)
*
* @note The fast index is the last index. Please take care when iterating
* through the array.
*/
template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
return data_[element_offset(strides_, index...)];
}
/*
* @brief Access element at given multi-dimensional index (const version).
* i.e. arr(i,j,k,...)
*
* @note The fast index is the last index. Please take care when iterating
* through the array.
*/
template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, const T &>
operator()(Ix... index) const {
return data_[element_offset(strides_, index...)];
}
/*
@brief Index the array as it would be a 1D array. To get a certain
pixel in a multidimensional array use the (i,j,k,...) operator instead.
*/
T &operator()(ssize_t i) { return data_[i]; }
/*
@brief Index the array as it would be a 1D array. To get a certain
pixel in a multidimensional array use the (i,j,k,...) operator instead.
*/
const T &operator()(ssize_t i) const { return data_[i]; }
/*
@brief Index the array as it would be a 1D array. To get a certain
pixel in a multidimensional array use the (i,j,k,...) operator instead.
*/
T &operator[](ssize_t i) { return data_[i]; }
/*
@brief Index the array as it would be a 1D array. To get a certain
pixel in a multidimensional array use the (i,j,k,...) operator instead.
*/
const T &operator[](ssize_t i) const { return data_[i]; }
/* @brief Return a raw pointer to the data */
T *data() { return data_; }
/* @brief Return a const raw pointer to the data */
const T *data() const { return data_; }
/* @brief Return a byte pointer to the data. Useful for memcpy like
* operations */
std::byte *buffer() { return reinterpret_cast<std::byte *>(data_); }
/**
* @brief Return the total number of elements in the array as a signed
* integer
*/
ssize_t size() const { return static_cast<ssize_t>(size_); }
/** @brief Return the total number of bytes in the array */
size_t total_bytes() const { return size_ * sizeof(T); }
/** @brief Return the shape of the array */
Shape<Ndim> shape() const noexcept { return shape_; }
/** @brief Return the size of dimension i */
ssize_t shape(ssize_t i) const noexcept { return shape_[i]; }
/** @brief Return the strides of the array */
std::array<ssize_t, Ndim> strides() const noexcept { return strides_; }
/**
* @brief Return the bitdepth of the array. Useful for checking that
* detector data can fit in the array type.
*/
size_t bitdepth() const noexcept { return sizeof(T) * 8; }
/**
* @brief Return the number of bytes to step in each dimension when
* traversing the array.
*/
std::array<ssize_t, Ndim> byte_strides() const noexcept {
auto byte_strides = strides_;
for (auto &val : byte_strides)
val *= sizeof(T);
return byte_strides;
}
using value_type = T;
NDArray &operator=(NDArray &&other) noexcept; // Move assign
NDArray &operator=(const NDArray &other); // Copy assign
NDArray &operator+=(const NDArray &other);
NDArray &operator-=(const NDArray &other);
NDArray &operator*=(const NDArray &other);
///////////////////////////////////////////////////////////////////////////////
// Assignments
//
///////////////////////////////////////////////////////////////////////////////
// Write directly to the data array, or create a new one
/**
* @brief Copy to the NDArray from an std::array. If the size of the array
* is different we reallocate the data.
*
*/
template <size_t Size>
NDArray<T, 1> &operator=(const std::array<T, Size> &other) {
if (Size != size_) {
@@ -141,12 +276,94 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
return *this;
}
// NDArray& operator/=(const NDArray& other);
/**
* @brief Move assignment operator.
*/
NDArray &operator=(NDArray &&other) noexcept {
// TODO! Should we use swap?
if (this != &other) {
delete[] data_;
data_ = other.data_;
shape_ = other.shape_;
size_ = other.size_;
strides_ = other.strides_;
other.reset();
}
return *this;
}
/**
* @brief Copy assignment operator.
*/
NDArray &operator=(const NDArray &other) {
if (this != &other) {
delete[] data_;
shape_ = other.shape_;
strides_ = other.strides_;
size_ = other.size_;
data_ = new T[size_];
std::copy(other.data_, other.data_ + size_, data_);
}
return *this;
}
///////////////////////////////////////////////////////////////////////////////
// Math operators
//
///////////////////////////////////////////////////////////////////////////////
/**
* @brief Add elementwise from another NDArray.
*/
NDArray &operator+=(const NDArray &other) {
if (shape_ != other.shape_)
throw(std::runtime_error(
"Shape of NDArray must match for operator +="));
for (size_t i = 0; i < size_; ++i) {
data_[i] += other.data_[i];
}
return *this;
}
/**
* @brief Subtract elementwise with another NDArray.
*/
NDArray &operator-=(const NDArray &other) {
if (shape_ != other.shape_)
throw(std::runtime_error(
"Shape of NDArray must match for operator -="));
for (size_t i = 0; i < size_; ++i) {
data_[i] -= other.data_[i];
}
return *this;
}
/**
* @brief Multiply elementwise with another NDArray.
*/
NDArray &operator*=(const NDArray &other) {
if (shape_ != other.shape_)
throw(std::runtime_error(
"Shape of NDArray must match for operator *="));
for (size_t i = 0; i < size_; ++i) {
data_[i] *= other.data_[i];
}
return *this;
}
/**
* @brief Divide elementwise by another NDArray. Templated to allow division
* with different types.
*
* TODO! Why is this templated when the others are not?
*/
template <typename V> NDArray &operator/=(const NDArray<V, Ndim> &other) {
// check shape
if (shape_ == other.shape()) {
for (uint32_t i = 0; i < size_; ++i) {
for (size_t i = 0; i < size_; ++i) {
data_[i] /= other(i);
}
return *this;
@@ -154,67 +371,139 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
throw(std::runtime_error("Shape of NDArray must match"));
}
NDArray<bool, Ndim> operator>(const NDArray &other);
/**
* @brief Assign a scalar value to all elements in the NDArray.
*/
NDArray &operator=(const T &value) {
std::fill_n(data_, size_, value);
return *this;
}
bool operator==(const NDArray &other) const;
bool operator!=(const NDArray &other) const;
/**
* @brief Add a scalar value to all elements in the NDArray.
*/
NDArray &operator+=(const T &value) {
for (size_t i = 0; i < size_; ++i)
data_[i] += value;
return *this;
}
NDArray &operator=(const T & /*value*/);
NDArray &operator+=(const T & /*value*/);
NDArray operator+(const T & /*value*/);
NDArray &operator-=(const T & /*value*/);
NDArray operator-(const T & /*value*/);
NDArray &operator*=(const T & /*value*/);
NDArray operator*(const T & /*value*/);
NDArray &operator/=(const T & /*value*/);
NDArray operator/(const T & /*value*/);
/**
* @brief Subtract a scalar value to all elements in the NDArray.
*/
NDArray &operator-=(const T &value) {
for (size_t i = 0; i < size_; ++i)
data_[i] -= value;
return *this;
}
NDArray &operator&=(const T & /*mask*/);
/**
* @brief Multiply all elements in the NDArray with a scalar value
*/
NDArray &operator*=(const T &value) {
for (size_t i = 0; i < size_; ++i)
data_[i] *= value;
return *this;
}
/**
* @brief Divide all elements in the NDArray with a scalar value
*/
NDArray &operator/=(const T &value) {
for (size_t i = 0; i < size_; ++i)
data_[i] /= value;
return *this;
}
/**
* @brief Bitwise AND all elements in the NDArray with a scalar mask.
* Used for example to mask out gain bits for Jungfrau detectors.
*/
NDArray &operator&=(const T &mask) {
for (auto it = begin(); it != end(); ++it)
*it &= mask;
return *this;
}
/**
* @brief Operator + with a scalar value. Returns a new NDArray.
*
* TODO! Expression template version of this?
*/
NDArray operator+(const T &value) {
NDArray result = *this;
result += value;
return result;
}
/**
* @brief Operator - with a scalar value. Returns a new NDArray.
*
* TODO! Expression template version of this?
*/
NDArray operator-(const T &value) {
NDArray result = *this;
result -= value;
return result;
}
/**
* @brief Operator * with a scalar value. Returns a new NDArray.
*
* TODO! Expression template version of this?
*/
NDArray operator*(const T &value) {
NDArray result = *this;
result *= value;
return result;
}
/**
* @brief Operator / with a scalar value. Returns a new NDArray.
*
* TODO! Expression template version of this?
*/
NDArray operator/(const T &value) {
NDArray result = *this;
result /= value;
return result;
}
/**
* @brief Compare two NDArrays elementwise for equality.
*/
bool operator==(const NDArray &other) const {
if (shape_ != other.shape_)
return false;
for (size_t i = 0; i != size_; ++i)
if (data_[i] != other.data_[i])
return false;
return true;
}
/**
* @brief Compare two NDArrays elementwise for non-equality.
*/
bool operator!=(const NDArray &other) const { return !((*this) == other); }
/**
* @brief Compute the square root of all elements in the NDArray.
*/
void sqrt() {
for (int i = 0; i < size_; ++i) {
for (size_t i = 0; i < size_; ++i) {
data_[i] = std::sqrt(data_[i]);
}
}
NDArray &operator++(); // pre inc
template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
return data_[element_offset(strides_, index...)];
}
template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) const {
return data_[element_offset(strides_, index...)];
}
template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, T> value(Ix... index) {
return data_[element_offset(strides_, index...)];
}
// TODO! is int the right type for index?
T &operator()(ssize_t i) { return data_[i]; }
const T &operator()(ssize_t i) const { return data_[i]; }
T &operator[](ssize_t i) { return data_[i]; }
const T &operator[](ssize_t i) const { return data_[i]; }
T *data() { return data_; }
std::byte *buffer() { return reinterpret_cast<std::byte *>(data_); }
ssize_t size() const { return static_cast<ssize_t>(size_); }
size_t total_bytes() const { return size_ * sizeof(T); }
std::array<ssize_t, Ndim> shape() const noexcept { return shape_; }
ssize_t shape(ssize_t i) const noexcept { return shape_[i]; }
std::array<ssize_t, Ndim> strides() const noexcept { return strides_; }
size_t bitdepth() const noexcept { return sizeof(T) * 8; }
std::array<ssize_t, Ndim> byte_strides() const noexcept {
auto byte_strides = strides_;
for (auto &val : byte_strides)
val *= sizeof(T);
return byte_strides;
/*
* @brief Prefix increment operator. Increments all elements by 1.
*/
NDArray &operator++() {
for (size_t i = 0; i < size_; ++i)
data_[i] += T{1};
return *this;
}
/**
@@ -224,10 +513,12 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
*/
NDView<T, Ndim> view() const { return NDView<T, Ndim>{data_, shape_}; }
void Print();
void Print_all();
void Print_some();
private:
/**
* @brief Reset the NDArray to an empty state. Dropping the ownership of
* the data. Used internally for move operations to avoid double free or
* dangling pointers.
*/
void reset() {
data_ = nullptr;
size_ = 0;
@@ -236,167 +527,10 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
}
};
// Move assign
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &
NDArray<T, Ndim>::operator=(NDArray<T, Ndim> &&other) noexcept {
if (this != &other) {
delete[] data_;
data_ = other.data_;
shape_ = other.shape_;
size_ = other.size_;
strides_ = other.strides_;
other.reset();
}
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const NDArray<T, Ndim> &other) {
// check shape
if (shape_ == other.shape_) {
for (size_t i = 0; i < size_; ++i) {
data_[i] += other.data_[i];
}
return *this;
}
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const NDArray<T, Ndim> &other) {
// check shape
if (shape_ == other.shape_) {
for (uint32_t i = 0; i < size_; ++i) {
data_[i] -= other.data_[i];
}
return *this;
}
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const NDArray<T, Ndim> &other) {
// check shape
if (shape_ == other.shape_) {
for (uint32_t i = 0; i < size_; ++i) {
data_[i] *= other.data_[i];
}
return *this;
}
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator&=(const T &mask) {
for (auto it = begin(); it != end(); ++it)
*it &= mask;
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<bool, Ndim> NDArray<T, Ndim>::operator>(const NDArray &other) {
if (shape_ == other.shape_) {
NDArray<bool, Ndim> result{shape_};
for (int i = 0; i < size_; ++i) {
result(i) = (data_[i] > other.data_[i]);
}
return result;
}
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(const NDArray<T, Ndim> &other) {
if (this != &other) {
delete[] data_;
shape_ = other.shape_;
strides_ = other.strides_;
size_ = other.size_;
data_ = new T[size_];
std::copy(other.data_, other.data_ + size_, data_);
}
return *this;
}
template <typename T, ssize_t Ndim>
bool NDArray<T, Ndim>::operator==(const NDArray<T, Ndim> &other) const {
if (shape_ != other.shape_)
return false;
for (uint32_t i = 0; i != size_; ++i)
if (data_[i] != other.data_[i])
return false;
return true;
}
template <typename T, ssize_t Ndim>
bool NDArray<T, Ndim>::operator!=(const NDArray<T, Ndim> &other) const {
return !((*this) == other);
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator++() {
for (uint32_t i = 0; i < size_; ++i)
data_[i] += 1;
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(const T &value) {
std::fill_n(data_, size_, value);
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] += value;
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const T &value) {
NDArray result = *this;
result += value;
return result;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] -= value;
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator-(const T &value) {
NDArray result = *this;
result -= value;
return result;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator/=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] /= value;
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator/(const T &value) {
NDArray result = *this;
result /= value;
return result;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] *= value;
return *this;
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const T &value) {
NDArray result = *this;
result *= value;
return result;
}
///////////////////////////////////////////////////////////////////////////////
// Free functions closely related to NDArray
//
///////////////////////////////////////////////////////////////////////////////
template <typename T, ssize_t Ndim>
std::ostream &operator<<(std::ostream &os, const NDArray<T, Ndim> &arr) {
@@ -410,27 +544,9 @@ std::ostream &operator<<(std::ostream &os, const NDArray<T, Ndim> &arr) {
return os;
}
template <typename T, ssize_t Ndim> void NDArray<T, Ndim>::Print_all() {
for (auto row = 0; row < shape_[0]; ++row) {
for (auto col = 0; col < shape_[1]; ++col) {
std::cout << std::setw(3);
std::cout << (*this)(row, col) << " ";
}
std::cout << "\n";
}
}
template <typename T, ssize_t Ndim> void NDArray<T, Ndim>::Print_some() {
for (auto row = 0; row < 5; ++row) {
for (auto col = 0; col < 5; ++col) {
std::cout << std::setw(7);
std::cout << (*this)(row, col) << " ";
}
std::cout << "\n";
}
}
template <typename T, ssize_t Ndim>
void save(NDArray<T, Ndim> &img, std::string &pathname) {
[[deprecated("Saving of raw arrays without metadata is deprecated")]] void
save(NDArray<T, Ndim> &img, std::string &pathname) {
std::ofstream f;
f.open(pathname, std::ios::binary);
f.write(img.buffer(), img.size() * sizeof(T));
@@ -438,8 +554,9 @@ void save(NDArray<T, Ndim> &img, std::string &pathname) {
}
template <typename T, ssize_t Ndim>
NDArray<T, Ndim> load(const std::string &pathname,
std::array<ssize_t, Ndim> shape) {
[[deprecated(
"Loading of raw arrays without metadata is deprecated")]] NDArray<T, Ndim>
load(const std::string &pathname, std::array<ssize_t, Ndim> shape) {
NDArray<T, Ndim> img{shape};
std::ifstream f;
f.open(pathname, std::ios::binary);
@@ -448,6 +565,20 @@ NDArray<T, Ndim> load(const std::string &pathname,
return img;
}
/**
* @brief Free function to safely divide two NDArrays elementwise, handling
* division by zero. Uses static_cast to convert types as needed.
*
* @tparam RT Result type
* @tparam NT Numerator type
* @tparam DT Denominator type
* @tparam Ndim Number of dimensions
* @param numerator The numerator NDArray
* @param denominator The denominator NDArray
* @return NDArray<RT, Ndim> Resulting NDArray after safe division
* @throws std::runtime_error if the shapes of the numerator and denominator do
* not match
*/
template <typename RT, typename NT, typename DT, ssize_t Ndim>
NDArray<RT, Ndim> safe_divide(const NDArray<NT, Ndim> &numerator,
const NDArray<DT, Ndim> &denominator) {
@@ -467,4 +598,4 @@ NDArray<RT, Ndim> safe_divide(const NDArray<NT, Ndim> &numerator,
return result;
}
} // namespace aare
} // namespace aare

1
include/aare/NDView.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/ArrayExpr.hpp"
#include "aare/defs.hpp"

1
include/aare/NumpyFile.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Dtype.hpp"
#include "aare/FileInterface.hpp"

1
include/aare/NumpyHelpers.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <algorithm>

1
include/aare/Pedestal.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Frame.hpp"
#include "aare/NDArray.hpp"

1
include/aare/PixelMap.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/NDArray.hpp"

1
include/aare/ProducerConsumerQueue.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
/*
* Copyright (c) Meta Platforms, Inc. and affiliates.
*

1
include/aare/RawFile.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/DetectorGeometry.hpp"
#include "aare/FileInterface.hpp"

24
include/aare/RawMasterFile.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/defs.hpp"
#include <algorithm>
@@ -5,6 +6,7 @@
#include <fmt/format.h>
#include <fstream>
#include <optional>
#include <chrono>
#include <nlohmann/json.hpp>
using json = nlohmann::json;
@@ -42,14 +44,16 @@ class RawFileNameComponents {
class ScanParameters {
bool m_enabled = false;
std::string m_dac;
DACIndex m_dac{};
int m_start = 0;
int m_stop = 0;
int m_step = 0;
// TODO! add settleTime, requires string to time conversion
int64_t m_settleTime = 0; // [ns]
public:
ScanParameters(const std::string &par);
ScanParameters(const bool enabled, const DACIndex dac, const int start,
const int stop, const int step, const int64_t settleTime);
ScanParameters() = default;
ScanParameters(const ScanParameters &) = default;
ScanParameters &operator=(const ScanParameters &) = default;
@@ -57,8 +61,9 @@ class ScanParameters {
int start() const;
int stop() const;
int step() const;
const std::string &dac() const;
DACIndex dac() const;
bool enabled() const;
int64_t settleTime() const;
void increment_stop();
};
@@ -80,6 +85,9 @@ class RawMasterFile {
size_t m_bitdepth{};
uint8_t m_quad = 0;
std::optional<std::chrono::nanoseconds> m_exptime;
std::chrono::nanoseconds m_period{0};
xy m_geometry{};
xy m_udp_interfaces_per_module{1, 1};
@@ -99,11 +107,13 @@ class RawMasterFile {
std::optional<size_t> m_digital_samples;
std::optional<size_t> m_transceiver_samples;
std::optional<size_t> m_number_of_rows;
std::optional<uint8_t> m_counter_mask;
std::optional<ROI> m_roi;
public:
RawMasterFile(const std::filesystem::path &fpath);
RawMasterFile(std::istream &is, const std::string &fname); // for testing
std::filesystem::path data_fname(size_t mod_id, size_t file_id) const;
@@ -129,14 +139,18 @@ class RawMasterFile {
std::optional<size_t> digital_samples() const;
std::optional<size_t> transceiver_samples() const;
std::optional<size_t> number_of_rows() const;
std::optional<uint8_t> counter_mask() const;
std::optional<ROI> roi() const;
ScanParameters scan_parameters() const;
std::optional<std::chrono::nanoseconds> exptime() const { return m_exptime; }
std::chrono::nanoseconds period() const { return m_period; }
private:
void parse_json(const std::filesystem::path &fpath);
void parse_raw(const std::filesystem::path &fpath);
void parse_json(std::istream &is);
void parse_raw(std::istream &is);
void retrieve_geometry();
};

1
include/aare/RawSubFile.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Frame.hpp"
#include "aare/defs.hpp"

3
include/aare/VarClusterFinder.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <algorithm>
@@ -124,7 +125,7 @@ template <typename T> int VarClusterFinder<T>::check_neighbours(int i, int j) {
const auto row = i + di[k];
const auto col = j + dj[k];
if (row >= 0 && col >= 0 && row < shape_[0] && col < shape_[1]) {
auto tmp = labeled_.value(i + di[k], j + dj[k]);
auto tmp = labeled_(i + di[k], j + dj[k]);
if (tmp != 0)
neighbour_labels.push_back(tmp);
}

16
include/aare/algorithm.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <aare/NDArray.hpp>
@@ -109,4 +110,19 @@ template <typename Container> bool all_equal(const Container &c) {
return false;
}
/**
* linear interpolation
* @param bin_edge left and right bin edges
* @param bin_values function values at bin edges
* @param coord coordinate to interpolate at
* @return interpolated value at coord
*/
inline double linear_interpolation(const std::pair<double, double> &bin_edge,
const std::pair<double, double> &bin_values,
const double coord) {
const double bin_width = bin_edge.second - bin_edge.first;
return bin_values.first * (1 - (coord - bin_edge.first) / bin_width) +
bin_values.second * (coord - bin_edge.first) / bin_width;
}
} // namespace aare

1
include/aare/calibration.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/NDArray.hpp"

23
include/aare/decode.hpp
View File

@@ -1,10 +1,12 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/defs.hpp"
#include <aare/NDView.hpp>
#include <cstdint>
#include <vector>
namespace aare {
uint16_t adc_sar_05_decode64to16(uint64_t input);
uint16_t adc_sar_04_decode64to16(uint64_t input);
void adc_sar_05_decode64to16(NDView<uint64_t, 2> input,
@@ -12,6 +14,25 @@ void adc_sar_05_decode64to16(NDView<uint64_t, 2> input,
void adc_sar_04_decode64to16(NDView<uint64_t, 2> input,
NDView<uint16_t, 2> output);
/**
* @brief Called with a 32 bit unsigned integer, shift by offset
* and then return the lower 24 bits as an 32 bit integer
* @param input 32-ibt input value
* @param offset (should be in range 0-7 to allow for full 24 bits)
* @return uint32_t
*/
uint32_t mask32to24bits(uint32_t input, BitOffset offset={});
/**
* @brief Expand 24 bit values in a 8bit buffer to 32bit unsigned integers
* Used for detectors with 24bit counters in combination with CTB
*
* @param input View of the 24 bit data as uint8_t (no 24bit native data type exists)
* @param output Destination of the expanded data (32bit, unsigned)
* @param offset Offset within the first byte to where the data starts (0-7 bits)
*/
void expand24to32bit(NDView<uint8_t,1> input, NDView<uint32_t,1> output, BitOffset offset={});
/**
* @brief Apply custom weights to a 16-bit input value. Will sum up
* weights[i]**i for each bit i that is set in the input value.

156
include/aare/defs.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include "aare/Dtype.hpp"
@@ -215,22 +216,153 @@ enum class DetectorType {
Unknown
};
/**
* @brief Enum class to define the Digital to Analog converter
* The values are the same as in slsDetectorPackage
*/
enum DACIndex {
DAC_0,
DAC_1,
DAC_2,
DAC_3,
DAC_4,
DAC_5,
DAC_6,
DAC_7,
DAC_8,
DAC_9,
DAC_10,
DAC_11,
DAC_12,
DAC_13,
DAC_14,
DAC_15,
DAC_16,
DAC_17,
VSVP,
VTRIM,
VRPREAMP,
VRSHAPER,
VSVN,
VTGSTV,
VCMP_LL,
VCMP_LR,
VCAL,
VCMP_RL,
RXB_RB,
RXB_LB,
VCMP_RR,
VCP,
VCN,
VISHAPER,
VTHRESHOLD,
IO_DELAY,
VREF_DS,
VOUT_CM,
VIN_CM,
VREF_COMP,
VB_COMP,
VDD_PROT,
VIN_COM,
VREF_PRECH,
VB_PIXBUF,
VB_DS,
VREF_H_ADC,
VB_COMP_FE,
VB_COMP_ADC,
VCOM_CDS,
VREF_RSTORE,
VB_OPA_1ST,
VREF_COMP_FE,
VCOM_ADC1,
VREF_L_ADC,
VREF_CDS,
VB_CS,
VB_OPA_FD,
VCOM_ADC2,
VCASSH,
VTH2,
VRSHAPER_N,
VIPRE_OUT,
VTH3,
VTH1,
VICIN,
VCAS,
VCAL_N,
VIPRE,
VCAL_P,
VDCSH,
VBP_COLBUF,
VB_SDA,
VCASC_SFP,
VIPRE_CDS,
IBIAS_SFP,
ADC_VPP,
HIGH_VOLTAGE,
TEMPERATURE_ADC,
TEMPERATURE_FPGA,
TEMPERATURE_FPGAEXT,
TEMPERATURE_10GE,
TEMPERATURE_DCDC,
TEMPERATURE_SODL,
TEMPERATURE_SODR,
TEMPERATURE_FPGA2,
TEMPERATURE_FPGA3,
TRIMBIT_SCAN,
V_POWER_A = 100,
V_POWER_B = 101,
V_POWER_C = 102,
V_POWER_D = 103,
V_POWER_IO = 104,
V_POWER_CHIP = 105,
I_POWER_A = 106,
I_POWER_B = 107,
I_POWER_C = 108,
I_POWER_D = 109,
I_POWER_IO = 110,
V_LIMIT = 111,
SLOW_ADC0 = 1000,
SLOW_ADC1,
SLOW_ADC2,
SLOW_ADC3,
SLOW_ADC4,
SLOW_ADC5,
SLOW_ADC6,
SLOW_ADC7,
SLOW_ADC_TEMP
};
// helper pair class to easily expose in python
template <typename T1, typename T2> struct Sum_index_pair {
T1 sum;
T2 index;
};
enum class corner : int {
cTopLeft = 0,
cTopRight = 1,
cBottomLeft = 2,
cBottomRight = 3
};
enum class TimingMode { Auto, Trigger };
enum class FrameDiscardPolicy { NoDiscard, Discard, DiscardPartial };
template <class T> T StringTo(const std::string &arg) { return T(arg); }
template <class T> std::string ToString(T arg) { return T(arg); }
template <> DetectorType StringTo(const std::string & /*name*/);
template <> std::string ToString(DetectorType arg);
template <> TimingMode StringTo(const std::string & /*mode*/);
template <> FrameDiscardPolicy StringTo(const std::string & /*mode*/);
using DataTypeVariants = std::variant<uint16_t, uint32_t>;
constexpr uint16_t ADC_MASK = 0x3FFF; // used to mask out the gain bits in Jungfrau
constexpr uint16_t ADC_MASK =
0x3FFF; // used to mask out the gain bits in Jungfrau
class BitOffset{
uint8_t m_offset{};
public:
BitOffset() = default;
explicit BitOffset(uint32_t offset);
uint8_t value() const {return m_offset;}
bool operator==(const BitOffset& other) const;
bool operator<(const BitOffset& other) const;
};
} // namespace aare

3
include/aare/logger.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
/*Utility to log to console*/
@@ -105,7 +106,7 @@ class Logger {
}
std::ostringstream &Get() {
os << Color(m_level) << "- " << Timestamp() << " " << ToString(m_level)
os << Color(m_level) << "- " << Timestamp() << " " << Logger::ToString(m_level)
<< ": ";
return os;
}

1
include/aare/utils/ifstream_helpers.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <fstream>

1
include/aare/utils/par.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <thread>
#include <utility>

1
include/aare/utils/task.hpp
View File

@@ -1,3 +1,4 @@
// SPDX-License-Identifier: MPL-2.0
#pragma once
#include <utility>
#include <vector>