From 31ac10a8319b1becd4f16e5220d6fecc30d7058a Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Sun, 1 Jun 2025 17:38:01 +0200 Subject: [PATCH 1/2] ITS: simplify constants + mathutils Signed-off-by: Felix Schlepper --- .../GPU/ITStrackingGPU/TrackingKernels.h | 4 + .../GPU/ITStrackingGPU/VertexerTraitsGPU.h | 12 +-- .../ITS/tracking/GPU/cuda/TrackingKernels.cu | 30 ++---- .../tracking/GPU/cuda/VertexerTraitsGPU.cxx | 2 +- .../ITS/tracking/GPU/cuda/VertexingKernels.cu | 8 -- .../include/ITStracking/Configuration.h | 6 +- .../tracking/include/ITStracking/Constants.h | 102 ++---------------- .../include/ITStracking/IndexTableUtils.h | 4 +- .../tracking/include/ITStracking/MathUtils.h | 90 +++++++--------- .../include/ITStracking/TrackerTraits.h | 4 +- .../include/ITStracking/VertexerTraits.h | 2 - .../ITSMFT/ITS/tracking/src/Smoother.cxx | 4 +- .../ITSMFT/ITS/tracking/src/TimeFrame.cxx | 24 ++--- Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx | 8 +- .../ITSMFT/ITS/tracking/src/TrackerTraits.cxx | 55 ++++------ .../ITS/tracking/src/VertexerTraits.cxx | 12 +-- 16 files changed, 104 insertions(+), 263 deletions(-) diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h index 09c8c39725efa..b847aacd9bba5 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h +++ b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h @@ -22,7 +22,11 @@ class CellSeed; class ExternalAllocator; namespace gpu { + #ifdef GPUCA_GPUCODE // GPUg() global kernels must only when compiled by GPU compiler + +GPUdi() int4 getEmptyBinsRect() { return int4{0, 0, 0, 0}; } + GPUd() bool fitTrack(TrackITSExt& track, int start, int end, diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/VertexerTraitsGPU.h b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/VertexerTraitsGPU.h index a5c3709081a82..5b1d9194e1174 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/VertexerTraitsGPU.h +++ b/Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/VertexerTraitsGPU.h @@ -18,7 +18,6 @@ #define ITSTRACKINGGPU_VERTEXERTRAITSGPU_H_ #include -#include #include "ITStracking/VertexerTraits.h" #include "ITStracking/Configuration.h" @@ -29,13 +28,8 @@ #include "ITStrackingGPU/TimeFrameGPU.h" -namespace o2 +namespace o2::its { -namespace its -{ -class ROframe; - -using constants::its2::InversePhiBinSize; class VertexerTraitsGPU final : public VertexerTraits { @@ -63,6 +57,6 @@ inline void VertexerTraitsGPU::adoptTimeFrame(TimeFrame<7>* tf) noexcept mTimeFrame = static_cast*>(tf); } -} // namespace its -} // namespace o2 +} // namespace o2::its + #endif diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu index 18c89d39adda0..2191880374548 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu @@ -52,14 +52,8 @@ using namespace o2::track; namespace o2::its { -using namespace constants::its2; using Vertex = o2::dataformats::Vertex>; -GPUdii() float Sq(float v) -{ - return v * v; -} - namespace gpu { @@ -99,9 +93,9 @@ GPUd() const int4 getBinsRect(const Cluster& currentCluster, const int layerInde const float z1, const float z2, float maxdeltaz, float maxdeltaphi) { const float zRangeMin = o2::gpu::CAMath::Min(z1, z2) - maxdeltaz; - const float phiRangeMin = (maxdeltaphi > constants::math::Pi) ? 0.f : currentCluster.phi - maxdeltaphi; + const float phiRangeMin = (maxdeltaphi > o2::constants::math::PI) ? 0.f : currentCluster.phi - maxdeltaphi; const float zRangeMax = o2::gpu::CAMath::Max(z1, z2) + maxdeltaz; - const float phiRangeMax = (maxdeltaphi > constants::math::Pi) ? constants::math::TwoPi : currentCluster.phi + maxdeltaphi; + const float phiRangeMax = (maxdeltaphi > o2::constants::math::PI) ? o2::constants::math::TwoPI : currentCluster.phi + maxdeltaphi; if (zRangeMax < -utils.getLayerZ(layerIndex) || zRangeMin > utils.getLayerZ(layerIndex) || zRangeMin > zRangeMax) { @@ -129,7 +123,7 @@ GPUd() bool fitTrack(TrackITSExt& track, o2::base::PropagatorF::MatCorrType matCorrType) { for (int iLayer{start}; iLayer != end; iLayer += step) { - if (track.getClusterIndex(iLayer) == constants::its::UnusedIndex) { + if (track.getClusterIndex(iLayer) == constants::UnusedIndex) { continue; } const TrackingFrameInfo& trackingHit = tfInfos[iLayer][track.getClusterIndex(iLayer)]; @@ -316,7 +310,7 @@ GPUg() void fitTrackSeedsKernel( temporaryTrack.setChi2(0); int* clusters = seed.getClusters(); for (int iL{0}; iL < 7; ++iL) { - temporaryTrack.setExternalClusterIndex(iL, clusters[iL], clusters[iL] != constants::its::UnusedIndex); + temporaryTrack.setExternalClusterIndex(iL, clusters[iL], clusters[iL] != constants::UnusedIndex); } bool fitSuccess = fitTrack(temporaryTrack, // TrackITSExt& track, 0, // int lastLayer, @@ -422,8 +416,6 @@ GPUg() void computeLayerCellsKernel( const float cellDeltaTanLambdaSigma, const float nSigmaCut) { - constexpr float radl = 9.36f; // Radiation length of Si [cm]. - constexpr float rho = 2.33f; // Density of Si [g/cm^3]. constexpr float layerxX0[7] = {5.e-3f, 5.e-3f, 5.e-3f, 1.e-2f, 1.e-2f, 1.e-2f, 1.e-2f}; // Hardcoded here for the moment. for (int iCurrentTrackletIndex = blockIdx.x * blockDim.x + threadIdx.x; iCurrentTrackletIndex < nTrackletsCurrent; iCurrentTrackletIndex += blockDim.x * gridDim.x) { const Tracklet& currentTracklet = tracklets[layer][iCurrentTrackletIndex]; @@ -462,7 +454,7 @@ GPUg() void computeLayerCellsKernel( break; } - if (!track.correctForMaterial(layerxX0[layer + iC], layerxX0[layer] * radl * rho, true)) { + if (!track.correctForMaterial(layerxX0[layer + iC], layerxX0[layer] * constants::Radl * constants::Rho, true)) { break; } @@ -548,12 +540,12 @@ GPUg() void computeLayerTrackletsMultiROFKernel( if (primaryVertex.isFlagSet(2) && iteration != 3) { continue; } - const float resolution = o2::gpu::CAMath::Sqrt(Sq(resolutionPV) / primaryVertex.getNContributors() + Sq(positionResolution)); + const float resolution = o2::gpu::CAMath::Sqrt(math_utils::Sq(resolutionPV) / primaryVertex.getNContributors() + math_utils::Sq(positionResolution)); const float tanLambda{(currentCluster.zCoordinate - primaryVertex.getZ()) * inverseR0}; const float zAtRmin{tanLambda * (minR - currentCluster.radius) + currentCluster.zCoordinate}; const float zAtRmax{tanLambda * (maxR - currentCluster.radius) + currentCluster.zCoordinate}; - const float sqInverseDeltaZ0{1.f / (Sq(currentCluster.zCoordinate - primaryVertex.getZ()) + 2.e-8f)}; /// protecting from overflows adding the detector resolution - const float sigmaZ{o2::gpu::CAMath::Sqrt(Sq(resolution) * Sq(tanLambda) * ((Sq(inverseR0) + sqInverseDeltaZ0) * Sq(meanDeltaR) + 1.f) + Sq(meanDeltaR * MSAngle))}; + const float sqInverseDeltaZ0{1.f / (math_utils::Sq(currentCluster.zCoordinate - primaryVertex.getZ()) + 2.e-8f)}; /// protecting from overflows adding the detector resolution + const float sigmaZ{o2::gpu::CAMath::Sqrt(math_utils::Sq(resolution) * math_utils::Sq(tanLambda) * ((math_utils::Sq(inverseR0) + sqInverseDeltaZ0) * math_utils::Sq(meanDeltaR) + 1.f) + math_utils::Sq(meanDeltaR * MSAngle))}; const int4 selectedBinsRect{getBinsRect(currentCluster, layerIndex + 1, *utils, zAtRmin, zAtRmax, sigmaZ * NSigmaCut, phiCut)}; if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) { continue; @@ -587,7 +579,7 @@ GPUg() void computeLayerTrackletsMultiROFKernel( const float deltaPhi{o2::gpu::CAMath::Abs(currentCluster.phi - nextCluster.phi)}; const float deltaZ{o2::gpu::CAMath::Abs(tanLambda * (nextCluster.radius - currentCluster.radius) + currentCluster.zCoordinate - nextCluster.zCoordinate)}; const int nextSortedIndex{ROFClusters[layerIndex + 1][rof1] + nextClusterIndex}; - if (deltaZ / sigmaZ < NSigmaCut && (deltaPhi < phiCut || o2::gpu::CAMath::Abs(deltaPhi - constants::math::TwoPi) < phiCut)) { + if (deltaZ / sigmaZ < NSigmaCut && (deltaPhi < phiCut || o2::gpu::CAMath::Abs(deltaPhi - o2::constants::math::TwoPI) < phiCut)) { if constexpr (initRun) { trackletsLUT[layerIndex][currentSortedIndex]++; // we need l0 as well for usual exclusive sums. } else { @@ -634,8 +626,6 @@ GPUg() void processNeighboursKernel(const int layer, const o2::base::Propagator* propagator, const o2::base::PropagatorF::MatCorrType matCorrType) { - constexpr float radl = 9.36f; // Radiation length of Si [cm]. - constexpr float rho = 2.33f; // Density of Si [g/cm^3]. constexpr float layerxX0[7] = {5.e-3f, 5.e-3f, 5.e-3f, 1.e-2f, 1.e-2f, 1.e-2f, 1.e-2f}; // Hardcoded here for the moment. for (unsigned int iCurrentCell = blockIdx.x * blockDim.x + threadIdx.x; iCurrentCell < nCurrentCells; iCurrentCell += blockDim.x * gridDim.x) { int foundSeeds{0}; @@ -678,7 +668,7 @@ GPUg() void processNeighboursKernel(const int layer, } if (matCorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) { - if (!seed.correctForMaterial(layerxX0[layer - 1], layerxX0[layer - 1] * radl * rho, true)) { + if (!seed.correctForMaterial(layerxX0[layer - 1], layerxX0[layer - 1] * constants::Radl * constants::Rho, true)) { continue; } } diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/VertexerTraitsGPU.cxx b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/VertexerTraitsGPU.cxx index 2a6debe8f652e..90d654a26a43d 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/VertexerTraitsGPU.cxx +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/VertexerTraitsGPU.cxx @@ -40,7 +40,7 @@ void VertexerTraitsGPU::updateVertexingParameters(const std::vector(std::ceil(mIndexTableUtils.getNphiBins() * par.phiCut / constants::math::TwoPi)); + par.phiSpan = static_cast(std::ceil(mIndexTableUtils.getNphiBins() * par.phiCut / o2::constants::math::TwoPI)); par.zSpan = static_cast(std::ceil(par.zCut * mIndexTableUtils.getInverseZCoordinate(0))); } } diff --git a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/VertexingKernels.cu b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/VertexingKernels.cu index acbd77585df37..126e799efce5d 100644 --- a/Detectors/ITSMFT/ITS/tracking/GPU/cuda/VertexingKernels.cu +++ b/Detectors/ITSMFT/ITS/tracking/GPU/cuda/VertexingKernels.cu @@ -20,10 +20,7 @@ namespace o2 { namespace its { -using constants::its::VertexerHistogramVolume; -using constants::math::TwoPi; using math_utils::getNormalizedPhi; -using namespace constants::its2; namespace gpu { @@ -58,11 +55,6 @@ void trackletFinderHandler(const Cluster* clustersNextLayer, // 0 2 maxTrackletsPerCluster); // const unsigned int maxTrackletsPerCluster = 1e2 } /* -GPUd() float smallestAngleDifference(float a, float b) -{ - float diff = fmod(b - a + constants::math::Pi, constants::math::TwoPi) - constants::math::Pi; - return (diff < -constants::math::Pi) ? diff + constants::math::TwoPi : ((diff > constants::math::Pi) ? diff - constants::math::TwoPi : diff); -} GPUd() const int4 getBinsRect(const Cluster& currentCluster, const int layerIndex, const float z1, float maxdeltaz, float maxdeltaphi) diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h index d3f7597ae314b..14edd0b81e049 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h @@ -58,9 +58,9 @@ class Configuration : public Param }; struct TrackingParameters { - int CellMinimumLevel() { return MinTrackLength - constants::its::ClustersPerCell + 1; } - int CellsPerRoad() const { return NLayers - 2; } - int TrackletsPerRoad() const { return NLayers - 1; } + int CellMinimumLevel() const noexcept { return MinTrackLength - constants::ClustersPerCell + 1; } + int CellsPerRoad() const noexcept { return NLayers - 2; } + int TrackletsPerRoad() const noexcept { return NLayers - 1; } std::string asString() const; int NLayers = 7; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Constants.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Constants.h index c29ad2e01c588..9b932352ee8e4 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Constants.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Constants.h @@ -17,112 +17,24 @@ #define TRACKINGITSU_INCLUDE_CONSTANTS_H_ #include "ITStracking/Definitions.h" -#include "CommonConstants/MathConstants.h" -#include "GPUCommonMath.h" -#include "GPUCommonDef.h" - -#ifndef GPUCA_GPUCODE_DEVICE -#include -#include -#include -#endif - -namespace o2 -{ -namespace its -{ - -namespace constants +namespace o2::its::constants { constexpr float MB = 1024.f * 1024.f; constexpr float GB = 1024.f * 1024.f * 1024.f; constexpr bool DoTimeBenchmarks = true; constexpr bool SaveTimeBenchmarks = false; -namespace math -{ -constexpr float Pi{3.14159265359f}; -constexpr float TwoPi{2.0f * Pi}; -constexpr float FloatMinThreshold{1e-20f}; -} // namespace math - -namespace its -{ -constexpr int LayersNumberVertexer{3}; constexpr int ClustersPerCell{3}; constexpr int UnusedIndex{-1}; constexpr float Resolution{0.0005f}; - -GPUhdi() constexpr std::array VertexerHistogramVolume() -{ - return std::array{{1.98, 1.98, 40.f}}; -} -} // namespace its - -namespace its2 -{ -constexpr int LayersNumber{7}; -constexpr int TrackletsPerRoad{LayersNumber - 1}; -constexpr int CellsPerRoad{LayersNumber - 2}; - -GPUhdi() constexpr std::array LayersZCoordinate() -{ - constexpr double s = 1.; // safety margin - return std::array{16.333f + s, 16.333f + s, 16.333f + s, 42.140f + s, 42.140f + s, 73.745f + s, 73.745f + s}; -} - -GPUhdi() constexpr std::array LayersRCoordinate() +constexpr float Radl = 9.36f; // Radiation length of Si [cm] +constexpr float Rho = 2.33f; // Density of Si [g/cm^3] +namespace its // to be removed { - return std::array{{2.33959f, 3.14076f, 3.91924f, 19.6213f, 24.5597f, 34.388f, 39.3329f}}; -} - -constexpr int ZBins{256}; -constexpr int PhiBins{128}; -constexpr float InversePhiBinSize{PhiBins / constants::math::TwoPi}; -GPUhdi() constexpr std::array InverseZBinSize() -{ - constexpr auto zSize = LayersZCoordinate(); - return std::array{0.5f * ZBins / (zSize[0]), 0.5f * ZBins / (zSize[1]), 0.5f * ZBins / (zSize[2]), - 0.5f * ZBins / (zSize[3]), 0.5f * ZBins / (zSize[4]), 0.5f * ZBins / (zSize[5]), - 0.5f * ZBins / (zSize[6])}; -} - -GPUhdi() constexpr float getInverseZCoordinate(const int layerIndex) -{ - return 0.5f * ZBins / LayersZCoordinate()[layerIndex]; -} - -GPUhdi() int getZBinIndex(const int layerIndex, const float zCoordinate) -{ - return (zCoordinate + LayersZCoordinate()[layerIndex]) * - InverseZBinSize()[layerIndex]; -} - -GPUhdi() int getPhiBinIndex(const float currentPhi) -{ - return (currentPhi * InversePhiBinSize); -} - -GPUhdi() int getBinIndex(const int zIndex, const int phiIndex) -{ - return o2::gpu::GPUCommonMath::Min(phiIndex * ZBins + zIndex, - ZBins * PhiBins - 1); -} - -GPUhdi() constexpr int4 getEmptyBinsRect() { return int4{0, 0, 0, 0}; } - -} // namespace its2 - -namespace pdgcodes -{ -constexpr int PionCode{211}; -} -} // namespace constants -#ifndef GPUCA_GPUCODE_DEVICE -typedef std::vector> index_table_t; -#endif +constexpr int UnusedIndex{-1}; +constexpr float Resolution{0.0005f}; } // namespace its -} // namespace o2 +} // namespace o2::its::constants #endif /* TRACKINGITSU_INCLUDE_CONSTANTS_H_ */ diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/IndexTableUtils.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/IndexTableUtils.h index ed4027f77f360..61072cb2410b7 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/IndexTableUtils.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/IndexTableUtils.h @@ -16,9 +16,9 @@ #ifndef TRACKINGITSU_INCLUDE_INDEXTABLEUTILS_H_ #define TRACKINGITSU_INCLUDE_INDEXTABLEUTILS_H_ -#include "ITStracking/Constants.h" #include "ITStracking/Configuration.h" #include "ITStracking/Definitions.h" +#include "CommonConstants/MathConstants.h" #include "GPUCommonMath.h" #include "GPUCommonDef.h" @@ -55,7 +55,7 @@ class IndexTableUtils template inline void IndexTableUtils::setTrackingParameters(const T& params) { - mInversePhiBinSize = params.PhiBins / constants::math::TwoPi; + mInversePhiBinSize = params.PhiBins / o2::constants::math::TwoPI; mNzBins = params.ZBins; mNphiBins = params.PhiBins; for (int iLayer{0}; iLayer < params.LayerZ.size(); ++iLayer) { diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h index 9093609144283..bc99822d4b2d2 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h @@ -13,75 +13,48 @@ /// \brief /// -#ifndef TRACKINGITSU_INCLUDE_CAUTILS_H_ -#define TRACKINGITSU_INCLUDE_CAUTILS_H_ - -#ifndef GPUCA_GPUCODE_DEVICE -#include -#include -#include -#include -#endif +#ifndef O2_ITS_TRACKING_MATHUTILS_H_ +#define O2_ITS_TRACKING_MATHUTILS_H_ +#include "CommonConstants/MathConstants.h" #include "MathUtils/Utils.h" -#include "ITStracking/Constants.h" #include "GPUCommonMath.h" #include "GPUCommonDef.h" -namespace o2 -{ -namespace its -{ - -namespace math_utils +namespace o2::its::math_utils { -GPUhdni() float computePhi(const float, const float); -GPUhdni() float hypot(const float, const float); -GPUhdni() constexpr float getNormalizedPhi(const float); -GPUhdni() constexpr float3 crossProduct(const float3&, const float3&); -GPUhdni() float computeCurvature(float x1, float y1, float x2, float y2, float x3, float y3); -GPUhdni() float computeCurvatureCentreX(float x1, float y1, float x2, float y2, float x3, float y3); -GPUhdni() float computeTanDipAngle(float x1, float y1, float x2, float y2, float z1, float z2); - -} // namespace math_utils -GPUhdi() float math_utils::computePhi(const float x, const float y) +GPUhdi() float computePhi(float x, float y) { - //return o2::gpu::CAMath::ATan2(-yCoordinate, -xCoordinate) + constants::math::Pi; - return o2::math_utils::fastATan2(-y, -x) + constants::math::Pi; + return o2::math_utils::fastATan2(-y, -x) + o2::constants::math::PI; } -GPUhdi() float math_utils::hypot(const float x, const float y) +GPUhdi() constexpr float hypot(float x, float y) { - return o2::gpu::CAMath::Sqrt(x * x + y * y); + return o2::gpu::CAMath::Hypot(x, y); } -GPUhdi() constexpr float math_utils::getNormalizedPhi(const float phi) +GPUhdi() constexpr float getNormalizedPhi(float phi) { - return (phi < 0) ? phi + constants::math::TwoPi : (phi > constants::math::TwoPi) ? phi - constants::math::TwoPi - : phi; + phi -= o2::constants::math::TwoPI * o2::gpu::CAMath::Floor(phi * (1.f / o2::constants::math::TwoPI)); + return phi; } -GPUhdi() constexpr float3 math_utils::crossProduct(const float3& firstVector, const float3& secondVector) -{ - - return float3{(firstVector.y * secondVector.z) - (firstVector.z * secondVector.y), - (firstVector.z * secondVector.x) - (firstVector.x * secondVector.z), - (firstVector.x * secondVector.y) - (firstVector.y * secondVector.x)}; -} - -GPUhdi() float math_utils::computeCurvature(float x1, float y1, float x2, float y2, float x3, float y3) +GPUhdi() float computeCurvature(float x1, float y1, float x2, float y2, float x3, float y3) { const float d = (x2 - x1) * (y3 - y2) - (x3 - x2) * (y2 - y1); + if (o2::gpu::CAMath::Abs(d) < o2::constants::math::Almost0) { + return 0.f; + } const float a = 0.5f * ((y3 - y2) * (y2 * y2 - y1 * y1 + x2 * x2 - x1 * x1) - (y2 - y1) * (y3 * y3 - y2 * y2 + x3 * x3 - x2 * x2)); const float b = 0.5f * ((x2 - x1) * (y3 * y3 - y2 * y2 + x3 * x3 - x2 * x2) - (x3 - x2) * (y2 * y2 - y1 * y1 + x2 * x2 - x1 * x1)); - const float den2 = (d * x1 - a) * (d * x1 - a) + (d * y1 - b) * (d * y1 - b); - return den2 > 0.f ? -1.f * d / o2::gpu::CAMath::Sqrt(den2) : 0.f; + const float den = o2::gpu::CAMath::Hypot(d * x1 - a, d * y1 - b); + return -d / den; } -GPUhdi() float math_utils::computeCurvatureCentreX(float x1, float y1, float x2, float y2, float x3, float y3) +GPUhdi() float computeCurvatureCentreX(float x1, float y1, float x2, float y2, float x3, float y3) { float dx21 = x2 - x1, dx32 = x3 - x2; if (dx21 == 0.f || dx32 == 0.f) { // add small offset @@ -90,15 +63,30 @@ GPUhdi() float math_utils::computeCurvatureCentreX(float x1, float y1, float x2, dx32 = x3 - x2; } float k1 = (y2 - y1) / dx21, k2 = (y3 - y2) / dx32; - return (k1 != k2) ? 0.5f * (k1 * k2 * (y1 - y3) + k2 * (x1 + x2) - k1 * (x2 + x3)) / (k2 - k1) : 1e5; + return (k1 != k2) ? 0.5f * (k1 * k2 * (y1 - y3) + k2 * (x1 + x2) - k1 * (x2 + x3)) / (k2 - k1) : 1e5f; } -GPUhdi() float math_utils::computeTanDipAngle(float x1, float y1, float x2, float y2, float z1, float z2) +GPUhdi() float computeTanDipAngle(float x1, float y1, float x2, float y2, float z1, float z2) { - return (z1 - z2) / o2::gpu::CAMath::Sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)); + return (z1 - z2) / o2::gpu::CAMath::Hypot(x1 - x2, y1 - y2); } -} // namespace its -} // namespace o2 +GPUhdi() float smallestAngleDifference(float a, float b) +{ + return o2::gpu::CAMath::Remainderf(b - a, o2::constants::math::TwoPI); +} + +GPUhdi() float Sq(float v) +{ + return v * v; +} -#endif /* TRACKINGITSU_INCLUDE_CAUTILS_H_ */ +GPUhdi() float MSangle(float mass, float p, float xX0) +{ + float beta = p / o2::gpu::CAMath::Hypot(mass, p); + return 0.0136f * o2::gpu::CAMath::Sqrt(xX0) * (1.f + 0.038f * o2::gpu::CAMath::Log(xX0)) / (beta * p); +} + +} // namespace o2::its::math_utils + +#endif diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h index 5f4e40b92ba82..36956a5206277 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackerTraits.h @@ -113,9 +113,9 @@ template inline const int4 TrackerTraits::getBinsRect(const int layerIndex, float phi, float maxdeltaphi, float z1, float z2, float maxdeltaz) const noexcept { const float zRangeMin = o2::gpu::GPUCommonMath::Min(z1, z2) - maxdeltaz; - const float phiRangeMin = (maxdeltaphi > constants::math::Pi) ? 0.f : phi - maxdeltaphi; + const float phiRangeMin = (maxdeltaphi > o2::constants::math::PI) ? 0.f : phi - maxdeltaphi; const float zRangeMax = o2::gpu::GPUCommonMath::Max(z1, z2) + maxdeltaz; - const float phiRangeMax = (maxdeltaphi > constants::math::Pi) ? constants::math::TwoPi : phi + maxdeltaphi; + const float phiRangeMax = (maxdeltaphi > o2::constants::math::PI) ? o2::constants::math::TwoPI : phi + maxdeltaphi; if (zRangeMax < -mTrkParams[0].LayerZ[layerIndex] || zRangeMin > mTrkParams[0].LayerZ[layerIndex] || zRangeMin > zRangeMax) { diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexerTraits.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexerTraits.h index e4ecced6d67fb..e1e1d44e8ead9 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexerTraits.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/VertexerTraits.h @@ -41,8 +41,6 @@ class MCCompLabel; namespace its { -class ROframe; -using constants::its::LayersNumberVertexer; enum class TrackletMode { Layer0Layer1 = 0, diff --git a/Detectors/ITSMFT/ITS/tracking/src/Smoother.cxx b/Detectors/ITSMFT/ITS/tracking/src/Smoother.cxx index 6259e20a02dfb..9bc65161c3cbb 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Smoother.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Smoother.cxx @@ -50,7 +50,7 @@ Smoother::Smoother(TrackITSExt& track, size_t smoothingLayer, const ROframe& ////////////////////// // Outward propagation for (size_t iLayer{0}; iLayer < mLayerToSmooth; ++iLayer) { - if (mOutwardsTrack.getClusterIndex(iLayer) == constants::its::UnusedIndex) { // Shorter tracks + if (mOutwardsTrack.getClusterIndex(iLayer) == constants::UnusedIndex) { // Shorter tracks continue; } const TrackingFrameInfo& tF = event.getTrackingFrameInfoOnLayer(iLayer).at(mOutwardsTrack.getClusterIndex(iLayer)); @@ -78,7 +78,7 @@ Smoother::Smoother(TrackITSExt& track, size_t smoothingLayer, const ROframe& ///////////////////// // Inward propagation for (size_t iLayer{D - 1}; iLayer > mLayerToSmooth; --iLayer) { - if (mInwardsTrack.getClusterIndex(iLayer) == constants::its::UnusedIndex) { // Shorter tracks + if (mInwardsTrack.getClusterIndex(iLayer) == constants::UnusedIndex) { // Shorter tracks continue; } const TrackingFrameInfo& tF = event.getTrackingFrameInfoOnLayer(iLayer).at(mInwardsTrack.getClusterIndex(iLayer)); diff --git a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx index 34d8967c6a5bb..80dbae42fc387 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TimeFrame.cxx @@ -34,18 +34,6 @@ struct ClusterHelper { int bin; int ind; }; - -inline float MSangle(float mass, float p, float xX0) -{ - float beta = p / o2::gpu::CAMath::Hypot(mass, p); - return 0.0136f * o2::gpu::CAMath::Sqrt(xX0) * (1.f + 0.038f * o2::gpu::CAMath::Log(xX0)) / (beta * p); -} - -inline float Sq(float v) -{ - return v * v; -} - } // namespace namespace o2::its @@ -323,7 +311,7 @@ void TimeFrame::initialise(const int iteration, const TrackingParameter for (unsigned int iLayer{0}; iLayer < std::min((int)mClusters.size(), maxLayers); ++iLayer) { clearResizeBoundedVector(mClusters[iLayer], mUnsortedClusters[iLayer].size(), mMemoryPool.get()); clearResizeBoundedVector(mUsedClusters[iLayer], mUnsortedClusters[iLayer].size(), mMemoryPool.get()); - mPositionResolution[iLayer] = o2::gpu::CAMath::Sqrt(0.5 * (trkParam.SystErrorZ2[iLayer] + trkParam.SystErrorY2[iLayer]) + trkParam.LayerResolution[iLayer] * trkParam.LayerResolution[iLayer]); + mPositionResolution[iLayer] = o2::gpu::CAMath::Sqrt(0.5f * (trkParam.SystErrorZ2[iLayer] + trkParam.SystErrorY2[iLayer]) + trkParam.LayerResolution[iLayer] * trkParam.LayerResolution[iLayer]); } clearResizeBoundedArray(mIndexTables, mNrof * (trkParam.ZBins * trkParam.PhiBins + 1), mMemoryPool.get()); clearResizeBoundedVector(mLines, mNrof, mMemoryPool.get()); @@ -364,18 +352,18 @@ void TimeFrame::initialise(const int iteration, const TrackingParameter float oneOverR{0.001f * 0.3f * std::abs(mBz) / trkParam.TrackletMinPt}; for (unsigned int iLayer{0}; iLayer < mClusters.size(); ++iLayer) { - mMSangles[iLayer] = MSangle(0.14f, trkParam.TrackletMinPt, trkParam.LayerxX0[iLayer]); + mMSangles[iLayer] = math_utils::MSangle(0.14f, trkParam.TrackletMinPt, trkParam.LayerxX0[iLayer]); mPositionResolution[iLayer] = o2::gpu::CAMath::Sqrt(0.5f * (trkParam.SystErrorZ2[iLayer] + trkParam.SystErrorY2[iLayer]) + trkParam.LayerResolution[iLayer] * trkParam.LayerResolution[iLayer]); if (iLayer < mClusters.size() - 1) { const float& r1 = trkParam.LayerRadii[iLayer]; const float& r2 = trkParam.LayerRadii[iLayer + 1]; const float res1 = o2::gpu::CAMath::Hypot(trkParam.PVres, mPositionResolution[iLayer]); const float res2 = o2::gpu::CAMath::Hypot(trkParam.PVres, mPositionResolution[iLayer + 1]); - const float cosTheta1half = o2::gpu::CAMath::Sqrt(1.f - Sq(0.5f * r1 * oneOverR)); - const float cosTheta2half = o2::gpu::CAMath::Sqrt(1.f - Sq(0.5f * r2 * oneOverR)); + const float cosTheta1half = o2::gpu::CAMath::Sqrt(1.f - math_utils::Sq(0.5f * r1 * oneOverR)); + const float cosTheta2half = o2::gpu::CAMath::Sqrt(1.f - math_utils::Sq(0.5f * r2 * oneOverR)); float x = r2 * cosTheta1half - r1 * cosTheta2half; - float delta = o2::gpu::CAMath::Sqrt(1.f / (1.f - 0.25f * Sq(x * oneOverR)) * (Sq(0.25f * r1 * r2 * Sq(oneOverR) / cosTheta2half + cosTheta1half) * Sq(res1) + Sq(0.25f * r1 * r2 * Sq(oneOverR) / cosTheta1half + cosTheta2half) * Sq(res2))); - mPhiCuts[iLayer] = std::min(o2::gpu::CAMath::ASin(0.5f * x * oneOverR) + 2.f * mMSangles[iLayer] + delta, constants::math::Pi * 0.5f); + float delta = o2::gpu::CAMath::Sqrt(1.f / (1.f - 0.25f * math_utils::Sq(x * oneOverR)) * (math_utils::Sq(0.25f * r1 * r2 * math_utils::Sq(oneOverR) / cosTheta2half + cosTheta1half) * math_utils::Sq(res1) + math_utils::Sq(0.25f * r1 * r2 * math_utils::Sq(oneOverR) / cosTheta1half + cosTheta2half) * math_utils::Sq(res2))); + mPhiCuts[iLayer] = std::min(o2::gpu::CAMath::ASin(0.5f * x * oneOverR) + 2.f * mMSangles[iLayer] + delta, o2::constants::math::PI * 0.5f); } } diff --git a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx index 75c265dbdb703..f4da1a86818bb 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/Tracker.cxx @@ -178,7 +178,7 @@ void Tracker::computeRoadsMClabels() for (int iCell{0}; iCell < mTrkParams[0].CellsPerRoad(); ++iCell) { const int currentCellIndex{currentRoad[iCell]}; - if (currentCellIndex == constants::its::UnusedIndex) { + if (currentCellIndex == constants::UnusedIndex) { if (isFirstRoadCell) { continue; } else { @@ -270,7 +270,7 @@ void Tracker::computeTracksMClabels() for (int iCluster = 0; iCluster < TrackITSExt::MaxClusters; ++iCluster) { const int index = track.getClusterIndex(iCluster); - if (index == constants::its::UnusedIndex) { + if (index == constants::UnusedIndex) { continue; } auto labels = mTimeFrame->getClusterLabels(iCluster, index); @@ -300,7 +300,7 @@ void Tracker::computeTracksMClabels() // set fake clusters pattern for (int ic{TrackITSExt::MaxClusters}; ic--;) { auto clid = track.getClusterIndex(ic); - if (clid != constants::its::UnusedIndex) { + if (clid != constants::UnusedIndex) { auto labelsSpan = mTimeFrame->getClusterLabels(ic, clid); for (const auto& currentLabel : labelsSpan) { if (currentLabel == maxOccurrencesValue) { @@ -325,7 +325,7 @@ void Tracker::rectifyClusterIndices() for (auto& track : mTimeFrame->getTracks(iROF)) { for (int iCluster = 0; iCluster < TrackITSExt::MaxClusters; ++iCluster) { const int index = track.getClusterIndex(iCluster); - if (index != constants::its::UnusedIndex) { + if (index != constants::UnusedIndex) { track.setExternalClusterIndex(iCluster, mTimeFrame->getClusterExternalIndex(iCluster, index)); } } diff --git a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx index 766dc25cd6d8e..36636069137f3 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/TrackerTraits.cxx @@ -22,7 +22,7 @@ #include #endif -#include +#include #include #include "CommonConstants/MathConstants.h" @@ -38,18 +38,10 @@ using o2::base::PropagatorF; -namespace -{ -inline float Sq(float q) -{ - return q * q; -} -} // namespace - namespace o2::its { -constexpr int debugLevel{0}; +static constexpr int debugLevel{0}; template void TrackerTraits::computeLayerTracklets(const int iteration, int iROFslice, int iVertex) @@ -105,15 +97,15 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iROF if (primaryVertex.isFlagSet(2) && iteration != 3) { continue; } - const float resolution = o2::gpu::CAMath::Sqrt(Sq(mTrkParams[iteration].PVres) / primaryVertex.getNContributors() + Sq(mTimeFrame->getPositionResolution(iLayer))); + const float resolution = o2::gpu::CAMath::Sqrt(math_utils::Sq(mTrkParams[iteration].PVres) / primaryVertex.getNContributors() + math_utils::Sq(mTimeFrame->getPositionResolution(iLayer))); const float tanLambda{(currentCluster.zCoordinate - primaryVertex.getZ()) * inverseR0}; const float zAtRmin{tanLambda * (mTimeFrame->getMinR(iLayer + 1) - currentCluster.radius) + currentCluster.zCoordinate}; const float zAtRmax{tanLambda * (mTimeFrame->getMaxR(iLayer + 1) - currentCluster.radius) + currentCluster.zCoordinate}; - const float sqInverseDeltaZ0{1.f / (Sq(currentCluster.zCoordinate - primaryVertex.getZ()) + 2.e-8f)}; /// protecting from overflows adding the detector resolution - const float sigmaZ{o2::gpu::CAMath::Sqrt(Sq(resolution) * Sq(tanLambda) * ((Sq(inverseR0) + sqInverseDeltaZ0) * Sq(meanDeltaR) + 1.f) + Sq(meanDeltaR * mTimeFrame->getMSangle(iLayer)))}; + const float sqInverseDeltaZ0{1.f / (math_utils::Sq(currentCluster.zCoordinate - primaryVertex.getZ()) + 2.e-8f)}; /// protecting from overflows adding the detector resolution + const float sigmaZ{o2::gpu::CAMath::Sqrt(math_utils::Sq(resolution) * math_utils::Sq(tanLambda) * ((math_utils::Sq(inverseR0) + sqInverseDeltaZ0) * math_utils::Sq(meanDeltaR) + 1.f) + math_utils::Sq(meanDeltaR * mTimeFrame->getMSangle(iLayer)))}; const int4 selectedBinsRect{getBinsRect(currentCluster, iLayer + 1, zAtRmin, zAtRmax, sigmaZ * mTrkParams[iteration].NSigmaCut, mTimeFrame->getPhiCut(iLayer))}; if (selectedBinsRect.x == 0 && selectedBinsRect.y == 0 && selectedBinsRect.z == 0 && selectedBinsRect.w == 0) { @@ -181,7 +173,7 @@ void TrackerTraits::computeLayerTracklets(const int iteration, int iROF if (deltaZ / sigmaZ < mTrkParams[iteration].NSigmaCut && (deltaPhi < mTimeFrame->getPhiCut(iLayer) || - o2::gpu::GPUCommonMath::Abs(deltaPhi - constants::math::TwoPi) < mTimeFrame->getPhiCut(iLayer))) { + o2::gpu::GPUCommonMath::Abs(deltaPhi - o2::constants::math::TwoPI) < mTimeFrame->getPhiCut(iLayer))) { if (iLayer > 0) { mTimeFrame->getTrackletsLookupTable()[iLayer - 1][currentSortedIndex]++; } @@ -273,9 +265,6 @@ void TrackerTraits::computeLayerCells(const int iteration) std::ofstream off(std::format("cells{}.txt", iter++)); #endif - constexpr float radl = 9.36f; // Radiation length of Si [cm] - constexpr float rho = 2.33f; // Density of Si [g/cm^3] - for (int iLayer = 0; iLayer < mTrkParams[iteration].CellsPerRoad(); ++iLayer) { deepVectorClear(mTimeFrame->getCells()[iLayer]); if (iLayer > 0) { @@ -359,7 +348,7 @@ void TrackerTraits::computeLayerCells(const int iteration) break; } - if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer + iC], mTrkParams[0].LayerxX0[iLayer] * radl * rho, true)) { + if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer + iC], mTrkParams[0].LayerxX0[iLayer] * constants::Radl * constants::Rho, true)) { break; } @@ -441,7 +430,7 @@ void TrackerTraits::computeLayerCells(const int iteration) break; } - if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer + iC], mTrkParams[0].LayerxX0[iLayer] * radl * rho, true)) { + if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer + iC], mTrkParams[0].LayerxX0[iLayer] * constants::Radl * constants::Rho, true)) { break; } @@ -682,9 +671,7 @@ void TrackerTraits::processNeighbours(int iLayer, int iLevel, const bou } if (mCorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) { - float radl = 9.36f; // Radiation length of Si [cm] - float rho = 2.33f; // Density of Si [g/cm^3] - if (!seed.correctForMaterial(mTrkParams[0].LayerxX0[iLayer - 1], mTrkParams[0].LayerxX0[iLayer - 1] * radl * rho, true)) { + if (!seed.correctForMaterial(mTrkParams[0].LayerxX0[iLayer - 1], mTrkParams[0].LayerxX0[iLayer - 1] * constants::Radl * constants::Rho, true)) { continue; } } @@ -745,9 +732,7 @@ void TrackerTraits::processNeighbours(int iLayer, int iLevel, const bou } if (mCorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) { - float radl = 9.36f; // Radiation length of Si [cm] - float rho = 2.33f; // Density of Si [g/cm^3] - if (!seed.correctForMaterial(mTrkParams[0].LayerxX0[iLayer - 1], mTrkParams[0].LayerxX0[iLayer - 1] * radl * rho, true)) { + if (!seed.correctForMaterial(mTrkParams[0].LayerxX0[iLayer - 1], mTrkParams[0].LayerxX0[iLayer - 1] * constants::Radl * constants::Rho, true)) { continue; } } @@ -838,7 +823,7 @@ void TrackerTraits::findRoads(const int iteration) temporaryTrack.resetCovariance(); temporaryTrack.setChi2(0); for (int iL{0}; iL < 7; ++iL) { - temporaryTrack.setExternalClusterIndex(iL, seed.getCluster(iL), seed.getCluster(iL) != constants::its::UnusedIndex); + temporaryTrack.setExternalClusterIndex(iL, seed.getCluster(iL), seed.getCluster(iL) != constants::UnusedIndex); } bool fitSuccess = fitTrack(temporaryTrack, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF); @@ -874,7 +859,7 @@ void TrackerTraits::findRoads(const int iteration) trk.resetCovariance(); trk.setChi2(0); for (int iL{0}; iL < 7; ++iL) { - trk.setExternalClusterIndex(iL, seed.getCluster(iL), seed.getCluster(iL) != constants::its::UnusedIndex); + trk.setExternalClusterIndex(iL, seed.getCluster(iL), seed.getCluster(iL) != constants::UnusedIndex); } bool fitSuccess = fitTrack(trk, 0, mTrkParams[0].NLayers, 1, mTrkParams[0].MaxChi2ClusterAttachment, mTrkParams[0].MaxChi2NDF); @@ -898,7 +883,7 @@ void TrackerTraits::findRoads(const int iteration) int nShared = 0; bool isFirstShared{false}; for (int iLayer{0}; iLayer < mTrkParams[0].NLayers; ++iLayer) { - if (track.getClusterIndex(iLayer) == constants::its::UnusedIndex) { + if (track.getClusterIndex(iLayer) == constants::UnusedIndex) { continue; } nShared += int(mTimeFrame->isClusterUsed(iLayer, track.getClusterIndex(iLayer))); @@ -911,7 +896,7 @@ void TrackerTraits::findRoads(const int iteration) std::array rofs{INT_MAX, INT_MAX, INT_MAX}; for (int iLayer{0}; iLayer < mTrkParams[0].NLayers; ++iLayer) { - if (track.getClusterIndex(iLayer) == constants::its::UnusedIndex) { + if (track.getClusterIndex(iLayer) == constants::UnusedIndex) { continue; } mTimeFrame->markUsedCluster(iLayer, track.getClusterIndex(iLayer)); @@ -977,7 +962,7 @@ void TrackerTraits::extendTracks(const int iteration) track.setPattern(pattern); /// Make sure that the newly attached clusters get marked as used for (int iLayer{0}; iLayer < mTrkParams[iteration].NLayers; ++iLayer) { - if (track.getClusterIndex(iLayer) == constants::its::UnusedIndex) { + if (track.getClusterIndex(iLayer) == constants::UnusedIndex) { continue; } mTimeFrame->markUsedCluster(iLayer, track.getClusterIndex(iLayer)); @@ -1082,7 +1067,7 @@ bool TrackerTraits::fitTrack(TrackITSExt& track, int start, int end, in auto propInstance = o2::base::Propagator::Instance(); for (int iLayer{start}; iLayer != end; iLayer += step) { - if (track.getClusterIndex(iLayer) == constants::its::UnusedIndex) { + if (track.getClusterIndex(iLayer) == constants::UnusedIndex) { continue; } const TrackingFrameInfo& trackingHit = mTimeFrame->getTrackingFrameInfoOnLayer(iLayer)[track.getClusterIndex(iLayer)]; @@ -1096,9 +1081,7 @@ bool TrackerTraits::fitTrack(TrackITSExt& track, int start, int end, in } if (mCorrType == o2::base::PropagatorF::MatCorrType::USEMatCorrNONE) { - constexpr float radl = 9.36f; // Radiation length of Si [cm] - constexpr float rho = 2.33f; // Density of Si [g/cm^3] - if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer], mTrkParams[0].LayerxX0[iLayer] * radl * rho, true)) { + if (!track.correctForMaterial(mTrkParams[0].LayerxX0[iLayer], mTrkParams[0].LayerxX0[iLayer] * constants::Radl * constants::Rho, true)) { continue; } } @@ -1143,9 +1126,7 @@ bool TrackerTraits::trackFollowing(TrackITSExt* track, int rof, bool ou } if (mTrkParams[iteration].CorrType == PropagatorF::MatCorrType::USEMatCorrNONE) { // account for material affects if propagator does not - constexpr float radl = 9.36f; // Radiation length of Si [cm] - constexpr float rho = 2.33f; // Density of Si [g/cm^3] - if (!hypoParam.correctForMaterial(mTrkParams[iteration].LayerxX0[iLayer], mTrkParams[iteration].LayerxX0[iLayer] * radl * rho, true)) { + if (!hypoParam.correctForMaterial(mTrkParams[iteration].LayerxX0[iLayer], mTrkParams[iteration].LayerxX0[iLayer] * constants::Radl * constants::Rho, true)) { continue; } } diff --git a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx index fe1619efaa192..37b650c05bd61 100644 --- a/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx +++ b/Detectors/ITSMFT/ITS/tracking/src/VertexerTraits.cxx @@ -31,12 +31,6 @@ using namespace o2::its; -float smallestAngleDifference(float a, float b) -{ - float diff = fmod(b - a + constants::math::Pi, constants::math::TwoPi) - constants::math::Pi; - return (diff < -constants::math::Pi) ? diff + constants::math::TwoPi : ((diff > constants::math::Pi) ? diff - constants::math::TwoPi : diff); -} - template void trackleterKernelHost( const gsl::span& clustersNextLayer, // 0 2 @@ -75,7 +69,7 @@ void trackleterKernelHost( continue; } const Cluster& nextCluster{clustersNextLayer[iNextLayerClusterIndex]}; - if (o2::gpu::GPUCommonMath::Abs(smallestAngleDifference(currentCluster.phi, nextCluster.phi)) < phiCut) { + if (o2::gpu::GPUCommonMath::Abs(math_utils::smallestAngleDifference(currentCluster.phi, nextCluster.phi)) < phiCut) { if (storedTracklets < maxTrackletsPerCluster) { if constexpr (!EvalRun) { if constexpr (Mode == TrackletMode::Layer0Layer1) { @@ -128,7 +122,7 @@ void trackletSelectionKernelHost( continue; } const float deltaTanLambda{o2::gpu::GPUCommonMath::Abs(tracklet01.tanLambda - tracklet12.tanLambda)}; - const float deltaPhi{o2::gpu::GPUCommonMath::Abs(smallestAngleDifference(tracklet01.phi, tracklet12.phi))}; + const float deltaPhi{o2::gpu::GPUCommonMath::Abs(math_utils::smallestAngleDifference(tracklet01.phi, tracklet12.phi))}; if (!usedTracklets[iTracklet01] && deltaTanLambda < tanLambdaCut && deltaPhi < phiCut && validTracklets != maxTracklets) { usedClusters0[tracklet01.firstClusterIndex] = true; usedClusters2[tracklet12.secondClusterIndex] = true; @@ -171,7 +165,7 @@ void VertexerTraits::updateVertexingParameters(const std::vector(std::ceil(mIndexTableUtils.getNphiBins() * par.phiCut / constants::math::TwoPi)); + par.phiSpan = static_cast(std::ceil(mIndexTableUtils.getNphiBins() * par.phiCut / o2::constants::math::TwoPI)); par.zSpan = static_cast(std::ceil(par.zCut * mIndexTableUtils.getInverseZCoordinate(0))); } setNThreads(vrtPar[0].nThreads); From 5d3f3f1a85e0c6239c3f3c8182dad60f5d60066e Mon Sep 17 00:00:00 2001 From: Felix Schlepper Date: Tue, 10 Jun 2025 19:06:45 +0200 Subject: [PATCH 2/2] ITS: mathutils protect against degeneracy Signed-off-by: Felix Schlepper --- .../tracking/include/ITStracking/Constants.h | 1 + .../tracking/include/ITStracking/MathUtils.h | 25 +++++++++++++++---- 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Constants.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Constants.h index 9b932352ee8e4..48cc45e44cf1c 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Constants.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/Constants.h @@ -25,6 +25,7 @@ constexpr float GB = 1024.f * 1024.f * 1024.f; constexpr bool DoTimeBenchmarks = true; constexpr bool SaveTimeBenchmarks = false; +constexpr float Tolerance{1e-12}; // numerical tolerance constexpr int ClustersPerCell{3}; constexpr int UnusedIndex{-1}; constexpr float Resolution{0.0005f}; diff --git a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h index bc99822d4b2d2..c5c1e4a8ce220 100644 --- a/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h +++ b/Detectors/ITSMFT/ITS/tracking/include/ITStracking/MathUtils.h @@ -17,6 +17,7 @@ #define O2_ITS_TRACKING_MATHUTILS_H_ #include "CommonConstants/MathConstants.h" +#include "ITStracking/Constants.h" #include "MathUtils/Utils.h" #include "GPUCommonMath.h" #include "GPUCommonDef.h" @@ -42,8 +43,9 @@ GPUhdi() constexpr float getNormalizedPhi(float phi) GPUhdi() float computeCurvature(float x1, float y1, float x2, float y2, float x3, float y3) { + // in case the triangle is degenerate we return infinite curvature. const float d = (x2 - x1) * (y3 - y2) - (x3 - x2) * (y2 - y1); - if (o2::gpu::CAMath::Abs(d) < o2::constants::math::Almost0) { + if (o2::gpu::CAMath::Abs(d) < o2::its::constants::Tolerance) { return 0.f; } const float a = @@ -51,24 +53,37 @@ GPUhdi() float computeCurvature(float x1, float y1, float x2, float y2, float x3 const float b = 0.5f * ((x2 - x1) * (y3 * y3 - y2 * y2 + x3 * x3 - x2 * x2) - (x3 - x2) * (y2 * y2 - y1 * y1 + x2 * x2 - x1 * x1)); const float den = o2::gpu::CAMath::Hypot(d * x1 - a, d * y1 - b); + if (den < o2::its::constants::Tolerance) { + return 0.f; + } return -d / den; } GPUhdi() float computeCurvatureCentreX(float x1, float y1, float x2, float y2, float x3, float y3) { + // in case the triangle is degenerate we return set the centre to infinity. float dx21 = x2 - x1, dx32 = x3 - x2; - if (dx21 == 0.f || dx32 == 0.f) { // add small offset + if (o2::gpu::CAMath::Abs(dx21) < o2::its::constants::Tolerance || + o2::gpu::CAMath::Abs(dx32) < o2::its::constants::Tolerance) { // add small offset x2 += 1e-4; dx21 = x2 - x1; dx32 = x3 - x2; } - float k1 = (y2 - y1) / dx21, k2 = (y3 - y2) / dx32; - return (k1 != k2) ? 0.5f * (k1 * k2 * (y1 - y3) + k2 * (x1 + x2) - k1 * (x2 + x3)) / (k2 - k1) : 1e5f; + const float k1 = (y2 - y1) / dx21, k2 = (y3 - y2) / dx32; + if (o2::gpu::CAMath::Abs(k2 - k1) < o2::its::constants::Tolerance) { + return o2::constants::math::VeryBig; + } + return 0.5f * (k1 * k2 * (y1 - y3) + k2 * (x1 + x2) - k1 * (x2 + x3)) / (k2 - k1); } GPUhdi() float computeTanDipAngle(float x1, float y1, float x2, float y2, float z1, float z2) { - return (z1 - z2) / o2::gpu::CAMath::Hypot(x1 - x2, y1 - y2); + // in case the points vertically align we go to pos/neg inifinity. + const float d = o2::gpu::CAMath::Hypot(x1 - x2, y1 - y2); + if (o2::gpu::CAMath::Abs(d) < o2::its::constants::Tolerance) { + return ((z1 > z2) ? -1.f : 1.f) * o2::constants::math::VeryBig; + } + return (z1 - z2) / d; } GPUhdi() float smallestAngleDifference(float a, float b)