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

changed eta interpolation to take into account photon center

Browse Source
This commit is contained in:
mazzol_a
2025-10-06 17:25:05 +02:00
parent c0ee17275e
commit 77a9788891
4 changed files with 148 additions and 83 deletions
Download Patch File Download Diff File Expand all files Collapse all files

150
include/aare/CalculateEta.hpp
View File

@@ -7,10 +7,10 @@
namespace aare { namespace aare {
enum class corner : int { enum class corner : int {
cBottomLeft = 0, cTopLeft = 0,
cBottomRight = 1, cTopRight = 1,
cTopLeft = 2, cBottomLeft = 2,
cTopRight = 3 cBottomRight = 3
}; };
enum class pixel : int { enum class pixel : int {
@@ -58,90 +58,126 @@ template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType> typename CoordType>
Eta2<T> Eta2<T>
calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) { calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
Eta2<T> eta{};
auto max_sum = cl.max_sum_2x2(); assert(ClusterSizeX > 1 && ClusterSizeY > 1);
eta.sum = max_sum.first; Eta2<T> eta{};
auto c = max_sum.second;
size_t cluster_center_index = size_t cluster_center_index =
(ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX; (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
size_t index_bottom_left_max_2x2_subcluster = auto max_sum = cl.max_sum_2x2();
(int(c / (ClusterSizeX - 1))) * ClusterSizeX + c % (ClusterSizeX - 1); eta.sum = max_sum.first;
int c = max_sum.second;
// calculate direction of gradient // subcluster top right from center
switch (static_cast<corner>(c)) {
// check that cluster center is in max subcluster case (corner::cTopLeft):
if (cluster_center_index != index_bottom_left_max_2x2_subcluster && if ((cl.data[cluster_center_index - 1] +
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) 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]);
if ((cl.data[cluster_center_index - ClusterSizeX] +
cl.data[cluster_center_index]) != 0)
eta.y = static_cast<double>(
cl.data[cluster_center_index - ClusterSizeX]) /
static_cast<double>(
cl.data[cluster_center_index - ClusterSizeX] +
cl.data[cluster_center_index]);
eta.x = static_cast<double>(cl.data[cluster_center_index + 1]) / // dx = 0
static_cast<double>((cl.data[cluster_center_index + 1] + // dy = 0
cl.data[cluster_center_index])); break;
} else { case (corner::cTopRight):
if ((cl.data[cluster_center_index] + if (cl.data[cluster_center_index] + cl.data[cluster_center_index + 1] !=
cl.data[cluster_center_index - 1]) != 0) 0)
eta.x = static_cast<double>(cl.data[cluster_center_index]) / eta.x = static_cast<double>(cl.data[cluster_center_index]) /
static_cast<double>((cl.data[cluster_center_index - 1] + static_cast<double>(cl.data[cluster_center_index] +
cl.data[cluster_center_index])); cl.data[cluster_center_index + 1]);
} if ((cl.data[cluster_center_index - ClusterSizeX] +
if ((cluster_center_index - index_bottom_left_max_2x2_subcluster) / cl.data[cluster_center_index]) != 0)
ClusterSizeX < eta.y = static_cast<double>(
1) { cl.data[cluster_center_index - ClusterSizeX]) /
assert(cluster_center_index + ClusterSizeX < static_cast<double>(
ClusterSizeX * ClusterSizeY); // suppress warning cl.data[cluster_center_index - ClusterSizeX] +
cl.data[cluster_center_index]);
// dx = 1
// dy = 0
break;
case (corner::cBottomLeft):
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]);
if ((cl.data[cluster_center_index] + if ((cl.data[cluster_center_index] +
cl.data[cluster_center_index + ClusterSizeX]) != 0) cl.data[cluster_center_index + ClusterSizeX]) != 0)
eta.y = static_cast<double>(
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]) / eta.y = static_cast<double>(cl.data[cluster_center_index]) /
static_cast<double>( static_cast<double>(
(cl.data[cluster_center_index] + cl.data[cluster_center_index] +
cl.data[cluster_center_index - ClusterSizeX])); cl.data[cluster_center_index + ClusterSizeX]);
// dx = 0
// dy = 1
break;
case (corner::cBottomRight):
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] +
cl.data[cluster_center_index + 1]);
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]);
// dx = 1
// dy = 1
break;
} }
eta.c = c; // TODO only supported for 2x2 and 3x3 clusters -> at least no eta.c = c;
// underyling enum class
return eta; return eta;
} }
// TODO! Look up eta2 calculation - photon center should be top right corner // TODO! Look up eta2 calculation - photon center should be bottom right corner
template <typename T> template <typename T>
Eta2<T> calculate_eta2(const Cluster<T, 2, 2, int16_t> &cl) { Eta2<T> calculate_eta2(const Cluster<T, 2, 2, int16_t> &cl) {
Eta2<T> eta{}; Eta2<T> eta{};
if ((cl.data[0] + cl.data[1]) != 0) if ((cl.data[0] + cl.data[1]) != 0)
eta.x = static_cast<double>(cl.data[1]) / eta.x = static_cast<double>(cl.data[2]) /
(cl.data[0] + cl.data[1]); // between (0,1) the closer to zero (cl.data[2] + cl.data[3]); // between (0,1) the closer to zero
// left value probably larger // left value probably larger
if ((cl.data[0] + cl.data[2]) != 0) if ((cl.data[0] + cl.data[2]) != 0)
eta.y = static_cast<double>(cl.data[2]) / eta.y = static_cast<double>(cl.data[1]) /
(cl.data[0] + cl.data[2]); // between (0,1) the closer to zero (cl.data[1] + cl.data[3]); // between (0,1) the closer to zero
// bottom value probably larger // bottom value probably larger
eta.sum = cl.sum(); eta.sum = cl.sum();
return eta; return eta;
} }
// TODO generalize
template <typename T>
Eta2<T> calculate_eta2(const Cluster<T, 1, 2, int16_t> &cl) {
Eta2<T> eta{};
eta.x = 0;
eta.y = static_cast<double>(cl.data[0]) / cl.data[1];
eta.sum = cl.sum();
}
template <typename T>
Eta2<T> calculate_eta2(const Cluster<T, 2, 1, int16_t> &cl) {
Eta2<T> eta{};
eta.x = static_cast<double>(cl.data[0]) / cl.data[1];
eta.y = 0;
eta.sum = cl.sum();
}
// calculates Eta3 for 3x3 cluster based on code from analyze_cluster // calculates Eta3 for 3x3 cluster based on code from analyze_cluster
// TODO only supported for 3x3 Clusters // 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_eta3(const Cluster<T, 3, 3> &cl) {

51
include/aare/Cluster.hpp
View File

@@ -8,7 +8,6 @@
#pragma once #pragma once
#include "logger.hpp"
#include <algorithm> #include <algorithm>
#include <array> #include <array>
#include <cstdint> #include <cstdint>
@@ -19,7 +18,7 @@ namespace aare {
// requires clause c++20 maybe update // requires clause c++20 maybe update
template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY, template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
typename CoordType = int16_t> typename CoordType = uint16_t>
struct Cluster { struct Cluster {
static_assert(std::is_arithmetic_v<T>, "T needs to be an arithmetic type"); static_assert(std::is_arithmetic_v<T>, "T needs to be an arithmetic type");
@@ -39,6 +38,13 @@ struct Cluster {
T sum() const { return std::accumulate(data.begin(), data.end(), T{}); } T sum() const { return std::accumulate(data.begin(), data.end(), T{}); }
// TODO: handle 1 dimensional clusters
// TODO: change int to corner
/**
* @brief sum of 2x2 subcluster with highest energy
* @return photon energy of subcluster, 2x2 subcluster index relative to
* cluster center
*/
std::pair<T, int> max_sum_2x2() const { std::pair<T, int> max_sum_2x2() const {
if constexpr (cluster_size_x == 3 && cluster_size_y == 3) { if constexpr (cluster_size_x == 3 && cluster_size_y == 3) {
@@ -54,17 +60,38 @@ struct Cluster {
} else if constexpr (cluster_size_x == 2 && cluster_size_y == 2) { } 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 std::make_pair(data[0] + data[1] + data[2] + data[3], 0);
} else { } else {
constexpr size_t num_2x2_subclusters = constexpr size_t cluster_center_index =
(ClusterSizeX - 1) * (ClusterSizeY - 1); (ClusterSizeX / 2) + (ClusterSizeY / 2) * ClusterSizeX;
std::array<T, num_2x2_subclusters> sum_2x2_subcluster; std::array<T, 4> sum_2x2_subcluster{0};
for (size_t i = 0; i < ClusterSizeY - 1; ++i) { // subcluster top left from center
for (size_t j = 0; j < ClusterSizeX - 1; ++j) sum_2x2_subcluster[0] =
sum_2x2_subcluster[i * (ClusterSizeX - 1) + j] = data[cluster_center_index] + data[cluster_center_index - 1] +
data[i * ClusterSizeX + j] + data[cluster_center_index - ClusterSizeX] +
data[i * ClusterSizeX + j + 1] + data[cluster_center_index - 1 - ClusterSizeX];
data[(i + 1) * ClusterSizeX + j] + // subcluster top right from center
data[(i + 1) * ClusterSizeX + j + 1]; 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(), int index = std::max_element(sum_2x2_subcluster.begin(),

4
include/aare/ClusterFinder.hpp
View File

@@ -136,7 +136,7 @@ class ClusterFinder {
// don't have a photon // don't have a photon
int i = 0; int i = 0;
for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) { for (int ir = -dy; ir < dy + has_center_pixel_y; ir++) {
for (int ic = -dx; ic < dx + has_center_pixel_y; ic++) { for (int ic = -dx; ic < dx + has_center_pixel_x; ic++) {
if (ix + ic >= 0 && ix + ic < frame.shape(1) && if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
iy + ir >= 0 && iy + ir < frame.shape(0)) { iy + ir >= 0 && iy + ir < frame.shape(0)) {
CT tmp = CT tmp =
@@ -144,7 +144,7 @@ class ClusterFinder {
static_cast<CT>( static_cast<CT>(
m_pedestal.mean(iy + ir, ix + ic)); m_pedestal.mean(iy + ir, ix + ic));
cluster.data[i] = cluster.data[i] =
tmp; // Watch for out of bounds access tmp; // Watch for out of bounds access
} }
i++; i++;
} }

26
include/aare/Interpolator.hpp
View File

@@ -69,26 +69,27 @@ Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
// cBottomRight = 1, // cBottomRight = 1,
// cTopLeft = 2, // cTopLeft = 2,
// cTopRight = 3 // cTopRight = 3
// TODO: could also chaneg the sign of the eta calculation
switch (static_cast<corner>(eta.c)) { switch (static_cast<corner>(eta.c)) {
case corner::cTopLeft: case corner::cTopLeft:
dX = -1.; dX = 0.0;
dY = 0; dY = 0.0;
break; break;
case corner::cTopRight:; case corner::cTopRight:;
dX = 0; dX = 1.0;
dY = 0; dY = 0.0;
break; break;
case corner::cBottomLeft: case corner::cBottomLeft:
dX = -1.; dX = 0.0;
dY = -1.; dY = 1.0;
break; break;
case corner::cBottomRight: case corner::cBottomRight:
dX = 0.; dX = 1.0;
dY = -1.; dY = 1.0;
break; break;
} }
photon.x += m_ietax(ix, iy, ie) * 2 + dX; photon.x -= m_ietax(ix, iy, ie) * 2 - dX;
photon.y += m_ietay(ix, iy, ie) * 2 + dY; photon.y -= m_ietay(ix, iy, ie) * 2 - dY;
photons.push_back(photon); photons.push_back(photon);
} }
} else if (clusters.cluster_size_x() == 2 || } else if (clusters.cluster_size_x() == 2 ||
@@ -112,10 +113,11 @@ Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
auto ix = last_smaller(m_etabinsx, eta.x); auto ix = last_smaller(m_etabinsx, eta.x);
auto iy = last_smaller(m_etabinsy, eta.y); auto iy = last_smaller(m_etabinsy, eta.y);
photon.x += m_ietax(ix, iy, ie) * // TODO: why 2?
photon.x -= m_ietax(ix, iy, ie) *
2; // eta goes between 0 and 1 but we could move the hit 2; // eta goes between 0 and 1 but we could move the hit
// anywhere in the 2x2 // anywhere in the 2x2
photon.y += m_ietay(ix, iy, ie) * 2; photon.y -= m_ietay(ix, iy, ie) * 2;
photons.push_back(photon); photons.push_back(photon);
} }