From 10a8de1d0995d613deec9e5f5a4d99a38177986d Mon Sep 17 00:00:00 2001
From: Simon Schuett <s6004392@stud.uni-frankfurt.de>
Date: Tue, 4 Mar 2025 16:32:59 +0100
Subject: [PATCH] online: Fix float precision issues in time calculations STS
hit finder
Fix cluster and hit finding in the STS detector when using 32-bit floating-point
precision. The issue occurs when calculating time differences between digis or
clusters with large timestamps - 32-bit floats lack sufficient precision for
representing nanosecond timestamps in the mCBM setup.
Add helper methods GetTimeDiff() that first convert timestamps to int before
calculating the difference, preserving precision and sign. Replace direct
floating-point timestamp subtractions throughout the code with calls to these
helper methods.
This resolves hit count discrepancies between 32-bit and 64-bit float versions,
reducing the difference from ~2.6% to ~0.04%.
---
algo/detectors/sts/Hitfinder.cxx | 23 +++++++++--------------
algo/detectors/sts/Hitfinder.h | 9 +++++++++
2 files changed, 18 insertions(+), 14 deletions(-)
diff --git a/algo/detectors/sts/Hitfinder.cxx b/algo/detectors/sts/Hitfinder.cxx
index cede2b7fdc..3ddf2df2d6 100644
--- a/algo/detectors/sts/Hitfinder.cxx
+++ b/algo/detectors/sts/Hitfinder.cxx
@@ -172,14 +172,13 @@ XPU_D void sts::Hitfinder::FindClustersParallel(FindClusters::context& ctx) cons
const u32 nextChannelEnd = (channel + 2 < nChannels) ? channelOffsets[channel + 2] : nDigis;
- float firstPossibleTime = float(myDigi.GetTimeU32()) - deltaT;
-
// Binary search for the first possible digi
u32 start = channelOffsets[channel + 1];
u32 end = nextChannelEnd;
while (start < end) {
u32 mid = (start + end) / 2;
- if (float(digis[mid].GetTimeU32()) < firstPossibleTime) {
+ real timeDiffDigis = GetTimeDiff(myDigi, digis[mid]);
+ if (deltaT < timeDiffDigis) {
start = mid + 1;
}
else {
@@ -194,15 +193,11 @@ XPU_D void sts::Hitfinder::FindClustersParallel(FindClusters::context& ctx) cons
// This should never happen?
if (myDigi.GetChannel() >= otherDigi.GetChannel()) continue;
- if (float(myDigi.GetTimeU32()) + deltaT < float(otherDigi.GetTimeU32())) {
+ real timeDiffDigis = xpu::abs(GetTimeDiff(myDigi, otherDigi));
+ if (deltaT < timeDiffDigis) {
break;
}
- // This check is not necessary? We already found the first possible digi via binary search...
- if (float(myDigi.GetTimeU32()) - deltaT > float(otherDigi.GetTimeU32())) {
- continue;
- }
-
// How to handle noise? Digis in same channel within a small time window of the previous digi are discarded by old cluster finder
// How to handle if multiple digis try to connect to the same digi?
@@ -607,8 +602,7 @@ XPU_D void sts::Hitfinder::FindHits(FindHits::context& ctx) const
float maxSigmaBoth = 4.f * xpu::sqrt(maxTerrF * maxTerrF + maxTerrB * maxTerrB);
- int startB = 0;
- float firstPossibleTime = float(clusterIdxF[iClusterF].fTime) - maxSigmaBoth;
+ int startB = 0;
// Use binary search to find the first cluster on back side that can be matched
// with the current cluster on front side
@@ -616,7 +610,8 @@ XPU_D void sts::Hitfinder::FindHits(FindHits::context& ctx) const
int end = nClustersB;
while (end - start > 1) {
int mid = (start + end) / 2;
- if (float(clusterIdxB[mid].fTime) < firstPossibleTime) {
+ auto timeDiffCluster = GetTimeDiff(clusterIdxF[iClusterF], clusterIdxB[mid]);
+ if (maxSigmaBoth < timeDiffCluster) {
start = mid;
}
else {
@@ -642,7 +637,7 @@ XPU_D void sts::Hitfinder::FindHits(FindHits::context& ctx) const
ClusterIdx clsIdxB = clusterIdxB[iClusterB];
sts::Cluster clsDataB = clusterDataB[clsIdxB.fIdx];
- float timeDiff = float(clsIdxF.fTime) - float(clsIdxB.fTime);
+ auto timeDiff = GetTimeDiff(clsIdxF, clsIdxB);
if (timeDiff > 0 && timeDiff > maxSigmaBoth) {
startB++;
@@ -674,7 +669,7 @@ XPU_D void sts::Hitfinder::FindHits(FindHits::context& ctx) const
timeCut = params.timeCutClusterSig * xpu::sqrt(eF * eF + eB * eB);
}
- if (xpu::abs(float(clsIdxF.fTime) - float(clsIdxB.fTime)) > timeCut) {
+ if (xpu::abs(timeDiff) > timeCut) {
continue;
}
diff --git a/algo/detectors/sts/Hitfinder.h b/algo/detectors/sts/Hitfinder.h
index c96c552253..c0accf7499 100644
--- a/algo/detectors/sts/Hitfinder.h
+++ b/algo/detectors/sts/Hitfinder.h
@@ -387,6 +387,15 @@ namespace cbm::algo::sts
const ClusterIdx& idxF, const Cluster& clsF, const ClusterIdx& idxB, const sts::Cluster& clsB,
float du, float dv) const;
+ XPU_D float GetTimeDiff(const CbmStsDigi& d1, const CbmStsDigi& d2) const
+ {
+ // Preserve sign of difference by first casting to int.
+ // Can't cast to float immediately as not enough precision with 32 bit floats for large timestamps in mCBM setup.
+ return int(d1.GetTimeU32()) - int(d2.GetTimeU32());
+ }
+
+ XPU_D float GetTimeDiff(const ClusterIdx& d1, const ClusterIdx& d2) const { return int(d1.fTime) - int(d2.fTime); }
+
XPU_D void SaveMaxError(float errorValue, int iModule) const
{
float* maxError = &maxClusterTimeErrorByModuleSide[iModule];
--
GitLab