Skip to content

ITS: simplify constants + mathutils #14383

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mconcas merged 2 commits into from
Jun 11, 2025
Merged

ITS: simplify constants + mathutils #14383

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
31ac10a
ITS: simplify constants + mathutils
f3sch Jun 1, 2025
5d3f3f1
ITS: mathutils protect against degeneracy
f3sch Jun 10, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/TrackingKernels.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
12 changes: 3 additions & 9 deletions Detectors/ITSMFT/ITS/tracking/GPU/ITStrackingGPU/VertexerTraitsGPU.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,6 @@
#define ITSTRACKINGGPU_VERTEXERTRAITSGPU_H_

#include <vector>
#include <array>

#include "ITStracking/VertexerTraits.h"
#include "ITStracking/Configuration.h"
Expand All @@ -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
{
Expand Down Expand Up @@ -63,6 +57,6 @@ inline void VertexerTraitsGPU::adoptTimeFrame(TimeFrame<7>* tf) noexcept
mTimeFrame = static_cast<TimeFrame<7>*>(tf);
}

} // namespace its
} // namespace o2
} // namespace o2::its

#endif
30 changes: 10 additions & 20 deletions Detectors/ITSMFT/ITS/tracking/GPU/cuda/TrackingKernels.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,14 +52,8 @@ using namespace o2::track;

namespace o2::its
{
using namespace constants::its2;
using Vertex = o2::dataformats::Vertex<o2::dataformats::TimeStamp<int>>;

GPUdii() float Sq(float v)
{
return v * v;
}

namespace gpu
{

Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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)];
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -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;
}

Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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 {
Expand Down Expand Up @@ -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};
Expand Down Expand Up @@ -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;
}
}
Expand Down
2 changes: 1 addition & 1 deletion Detectors/ITSMFT/ITS/tracking/GPU/cuda/VertexerTraitsGPU.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ void VertexerTraitsGPU::updateVertexingParameters(const std::vector<VertexingPar
mTfGPUParams = tfPar;
mIndexTableUtils.setTrackingParameters(vrtPar[0]);
for (auto& par : mVrtParams) {
par.phiSpan = static_cast<int>(std::ceil(mIndexTableUtils.getNphiBins() * par.phiCut / constants::math::TwoPi));
par.phiSpan = static_cast<int>(std::ceil(mIndexTableUtils.getNphiBins() * par.phiCut / o2::constants::math::TwoPI));
par.zSpan = static_cast<int>(std::ceil(par.zCut * mIndexTableUtils.getInverseZCoordinate(0)));
}
}
Expand Down
8 changes: 0 additions & 8 deletions Detectors/ITSMFT/ITS/tracking/GPU/cuda/VertexingKernels.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
{
Expand Down Expand Up @@ -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)
Expand Down
6 changes: 3 additions & 3 deletions Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
103 changes: 8 additions & 95 deletions Detectors/ITSMFT/ITS/tracking/include/ITStracking/Constants.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,112 +17,25 @@
#define TRACKINGITSU_INCLUDE_CONSTANTS_H_

#include "ITStracking/Definitions.h"
#include "CommonConstants/MathConstants.h"

#include "GPUCommonMath.h"
#include "GPUCommonDef.h"

#ifndef GPUCA_GPUCODE_DEVICE
#include <climits>
#include <vector>
#include <array>
#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 float Tolerance{1e-12}; // numerical tolerance
constexpr int ClustersPerCell{3};
constexpr int UnusedIndex{-1};
constexpr float Resolution{0.0005f};

GPUhdi() constexpr std::array<float, 3> VertexerHistogramVolume()
{
return std::array<float, 3>{{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<float, LayersNumber> LayersZCoordinate()
{
constexpr double s = 1.; // safety margin
return std::array<float, LayersNumber>{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<float, LayersNumber> 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<float, LayersNumber>{{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<float, LayersNumber> InverseZBinSize()
{
constexpr auto zSize = LayersZCoordinate();
return std::array<float, LayersNumber>{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<std::vector<int>> 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_ */
4 changes: 2 additions & 2 deletions Detectors/ITSMFT/ITS/tracking/include/ITStracking/IndexTableUtils.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down Expand Up @@ -55,7 +55,7 @@ class IndexTableUtils
template <class T>
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) {
Expand Down
Loading