Skip to content
Snippets Groups Projects
Commit 10a8de1d authored by Simon Schuett's avatar Simon Schuett Committed by Volker Friese
Browse files

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%.
parent 6aa036a5
No related branches found
No related tags found
1 merge request!2044online: Fix float precision issues in time calculations in STS hit finder
Pipeline #33257 passed
...@@ -172,14 +172,13 @@ XPU_D void sts::Hitfinder::FindClustersParallel(FindClusters::context& ctx) cons ...@@ -172,14 +172,13 @@ XPU_D void sts::Hitfinder::FindClustersParallel(FindClusters::context& ctx) cons
const u32 nextChannelEnd = (channel + 2 < nChannels) ? channelOffsets[channel + 2] : nDigis; const u32 nextChannelEnd = (channel + 2 < nChannels) ? channelOffsets[channel + 2] : nDigis;
float firstPossibleTime = float(myDigi.GetTimeU32()) - deltaT;
// Binary search for the first possible digi // Binary search for the first possible digi
u32 start = channelOffsets[channel + 1]; u32 start = channelOffsets[channel + 1];
u32 end = nextChannelEnd; u32 end = nextChannelEnd;
while (start < end) { while (start < end) {
u32 mid = (start + end) / 2; u32 mid = (start + end) / 2;
if (float(digis[mid].GetTimeU32()) < firstPossibleTime) { real timeDiffDigis = GetTimeDiff(myDigi, digis[mid]);
if (deltaT < timeDiffDigis) {
start = mid + 1; start = mid + 1;
} }
else { else {
...@@ -194,15 +193,11 @@ XPU_D void sts::Hitfinder::FindClustersParallel(FindClusters::context& ctx) cons ...@@ -194,15 +193,11 @@ XPU_D void sts::Hitfinder::FindClustersParallel(FindClusters::context& ctx) cons
// This should never happen? // This should never happen?
if (myDigi.GetChannel() >= otherDigi.GetChannel()) continue; 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; 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 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? // 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 ...@@ -607,8 +602,7 @@ XPU_D void sts::Hitfinder::FindHits(FindHits::context& ctx) const
float maxSigmaBoth = 4.f * xpu::sqrt(maxTerrF * maxTerrF + maxTerrB * maxTerrB); float maxSigmaBoth = 4.f * xpu::sqrt(maxTerrF * maxTerrF + maxTerrB * maxTerrB);
int startB = 0; int startB = 0;
float firstPossibleTime = float(clusterIdxF[iClusterF].fTime) - maxSigmaBoth;
// Use binary search to find the first cluster on back side that can be matched // Use binary search to find the first cluster on back side that can be matched
// with the current cluster on front side // with the current cluster on front side
...@@ -616,7 +610,8 @@ XPU_D void sts::Hitfinder::FindHits(FindHits::context& ctx) const ...@@ -616,7 +610,8 @@ XPU_D void sts::Hitfinder::FindHits(FindHits::context& ctx) const
int end = nClustersB; int end = nClustersB;
while (end - start > 1) { while (end - start > 1) {
int mid = (start + end) / 2; int mid = (start + end) / 2;
if (float(clusterIdxB[mid].fTime) < firstPossibleTime) { auto timeDiffCluster = GetTimeDiff(clusterIdxF[iClusterF], clusterIdxB[mid]);
if (maxSigmaBoth < timeDiffCluster) {
start = mid; start = mid;
} }
else { else {
...@@ -642,7 +637,7 @@ XPU_D void sts::Hitfinder::FindHits(FindHits::context& ctx) const ...@@ -642,7 +637,7 @@ XPU_D void sts::Hitfinder::FindHits(FindHits::context& ctx) const
ClusterIdx clsIdxB = clusterIdxB[iClusterB]; ClusterIdx clsIdxB = clusterIdxB[iClusterB];
sts::Cluster clsDataB = clusterDataB[clsIdxB.fIdx]; sts::Cluster clsDataB = clusterDataB[clsIdxB.fIdx];
float timeDiff = float(clsIdxF.fTime) - float(clsIdxB.fTime); auto timeDiff = GetTimeDiff(clsIdxF, clsIdxB);
if (timeDiff > 0 && timeDiff > maxSigmaBoth) { if (timeDiff > 0 && timeDiff > maxSigmaBoth) {
startB++; startB++;
...@@ -674,7 +669,7 @@ XPU_D void sts::Hitfinder::FindHits(FindHits::context& ctx) const ...@@ -674,7 +669,7 @@ XPU_D void sts::Hitfinder::FindHits(FindHits::context& ctx) const
timeCut = params.timeCutClusterSig * xpu::sqrt(eF * eF + eB * eB); 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; continue;
} }
......
...@@ -387,6 +387,15 @@ namespace cbm::algo::sts ...@@ -387,6 +387,15 @@ namespace cbm::algo::sts
const ClusterIdx& idxF, const Cluster& clsF, const ClusterIdx& idxB, const sts::Cluster& clsB, const ClusterIdx& idxF, const Cluster& clsF, const ClusterIdx& idxB, const sts::Cluster& clsB,
float du, float dv) const; 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 XPU_D void SaveMaxError(float errorValue, int iModule) const
{ {
float* maxError = &maxClusterTimeErrorByModuleSide[iModule]; float* maxError = &maxClusterTimeErrorByModuleSide[iModule];
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment