diff --git a/Core/include/Acts/Seeding2/GbtsTrackingFilter.hpp b/Core/include/Acts/Seeding2/GbtsTrackingFilter.hpp index 8ccffe99b1a..f30b5ca7473 100644 --- a/Core/include/Acts/Seeding2/GbtsTrackingFilter.hpp +++ b/Core/include/Acts/Seeding2/GbtsTrackingFilter.hpp @@ -8,6 +8,7 @@ #pragma once +#include "Acts/Definitions/Units.hpp" #include "Acts/Seeding2/GbtsDataStorage.hpp" #include "Acts/Seeding2/GbtsGeometry.hpp" @@ -70,9 +71,9 @@ class GbtsTrackingFilter final { float radLen = 0.025; /// Measurement uncertainty in x direction. - float sigmaX = 0.08; + float sigmaX = 0.08 * Acts::UnitConstants::mm; /// Measurement uncertainty in y direction. - float sigmaY = 0.25; + float sigmaY = 0.25 * Acts::UnitConstants::mm; /// Measurement weight in x direction. float weightX = 0.5; @@ -88,9 +89,9 @@ class GbtsTrackingFilter final { float addHit = 14.0; /// Maximum track curvature. - float maxCurvature = 1e-3f; + float maxCurvature = 1e-3f / Acts::UnitConstants::mm; /// Maximum longitudinal impact parameter. - float maxZ0 = 170.0; + float maxZ0 = 170.0 * Acts::UnitConstants::mm; }; /// State for the tracking filter, containing edge states and a global diff --git a/Core/include/Acts/Seeding2/GraphBasedTrackSeeder.hpp b/Core/include/Acts/Seeding2/GraphBasedTrackSeeder.hpp index 4edba30d9b5..6cf16df65ef 100644 --- a/Core/include/Acts/Seeding2/GraphBasedTrackSeeder.hpp +++ b/Core/include/Acts/Seeding2/GraphBasedTrackSeeder.hpp @@ -51,10 +51,16 @@ class GraphBasedTrackSeeder { bool matchBeforeCreate = false; /// Use legacy tuning parameters. bool useOldTunings = false; + /// description here + bool validateTriplets = true; + /// description here + bool useAdaptiveCuts = true; /// Tau ratio cut threshold. float tauRatioCut = 0.007; /// Tau ratio precut threshold. float tauRatioPrecut = 0.009f; + /// description here + float tauRatioCorr = 0.006; /// Eta bin width override (0 uses default from connection file). // specify non-zero to override eta bin width from connection file (default // 0.2 in createLinkingScheme.py) @@ -66,9 +72,6 @@ class GraphBasedTrackSeeder { float minPt = 1.0f * UnitConstants::GeV; // graph building options - /// Transverse momentum coefficient (~0.3*B/2 - assumes nominal field of - /// 2*T). - double ptCoeff = 0.29997 * 1.9972 / 2.0; /// Use eta binning from geometry structure. bool useEtaBinning = true; /// Apply RZ cuts on doublets. @@ -76,17 +79,23 @@ class GraphBasedTrackSeeder { /// Maximum number of Gbts edges/doublets. std::uint32_t nMaxEdges = 2000000; /// Minimum delta radius between layers. - float minDeltaRadius = 2.0; + float minDeltaRadius = 2.0 * Acts::UnitConstants::mm; + /// Maximum d0 impact parameter when validating edge connection triplet + float d0Max = 3.0 * UnitConstants::mm; // Seed extraction options /// Minimum eta for edge masking. float edgeMaskMinEta = 1.5; /// Threshold for hit sharing between seeds. float hitShareThreshold = 0.49; - + /// Max seed eta value considered for splitting. + float maxSeedSplitEta = 0.6; + /// Max allowed curvature for seed self consistency check. + // Units of inverse meters + float maxInvRadDiff = 0.7e-2 / UnitConstants::m; // GbtsDataStorage options /// Maximum endcap cluster width. - float maxEndcapClusterWidth = 0.35; + float maxEndcapClusterWidth = 0.35 * Acts::UnitConstants::mm; }; /// Derived configuration struct that contains calculated parameters based on @@ -100,27 +109,56 @@ class GraphBasedTrackSeeder { float phiSliceWidth = std::numeric_limits::quiet_NaN(); }; - /// Seed metadata produced by the GBTS algorithm. - struct SeedProperties { + /// Optional inputs for variables passed in + /// or derived during runtime. + struct Options { + /// Constructor. + /// @param bFieldInZ_ the magnetic field in z + explicit Options(float bFieldInZ_); + + /// Magnetic field in z + /// units of GeV/(e*mm). + float bFieldInZ{}; + + /// Transverse momentum coefficient (~0.3*B/2 - assumes nominal field of + /// 2*T). + double ptCoeff{}; + }; + /// candidate seed metadata produced by the GBTS algorithm. + struct SeedCandidateProperties { /// Constructor. /// @param quality Seed quality score /// @param clone Clone flag - /// @param sps Space point indices - SeedProperties(float quality, std::int32_t clone, - std::vector sps) - : seedQuality(quality), isClone(clone), spacePoints(std::move(sps)) {} + /// @param sps Vector of pointers to actual space points + /// @param splitFlag used to flag if seed needs to be split in two + SeedCandidateProperties(float quality, std::int32_t clone, + std::vector sps, + std::uint32_t splitFlag) + : seedQuality(quality), + isClone(clone), + spacePoints(std::move(sps)), + forSeedSplitting(splitFlag) {} /// Seed quality score. float seedQuality{}; /// Clone flag. std::int32_t isClone{}; /// Space point indices. - std::vector spacePoints; + std::vector spacePoints; + /// Flag for seed splitting. + std::uint32_t forSeedSplitting{}; + }; + + /// Output seed metadata + struct OutputSeedProperties { + /// Constructor. + OutputSeedProperties(float Quality, std::vector sps) + : seedQuality(Quality), spacePoints(std::move(sps)) {} - /// Comparison operator. - /// @param o Other seed properties to compare - /// @return True if this is less than other - auto operator<=>(const SeedProperties& o) const = default; + /// Quality of seed. + float seedQuality{}; + /// Index of spacepoints in seed. + std::vector spacePoints; }; /// Sliding window in phi used to define range used for edge creation @@ -154,7 +192,8 @@ class GraphBasedTrackSeeder { SeedContainer2 createSeeds(const SpacePointContainer2& spacePoints, const GbtsRoiDescriptor& roi, std::uint32_t maxLayers, - const GbtsTrackingFilter& filter) const; + const GbtsTrackingFilter& filter, + Options options) const; private: DerivedConfig m_cfg; @@ -187,7 +226,7 @@ class GraphBasedTrackSeeder { /// @return Pair of edge count and maximum level std::pair buildTheGraph( const GbtsRoiDescriptor& roi, GbtsNodeStorage& nodeStorage, - std::vector& edgeStorage) const; + std::vector& edgeStorage, Options options) const; /// Run connected component analysis on the graph. /// @param nEdges Number of edges in the graph @@ -206,7 +245,7 @@ class GraphBasedTrackSeeder { void extractSeedsFromTheGraph(std::uint32_t maxLevel, std::uint32_t nEdges, std::int32_t nHits, std::vector& edgeStorage, - std::vector& vSeedCandidates, + std::vector& vOutputSeeds, const GbtsTrackingFilter& filter) const; /// Check to see if z0 of segment is within the expected z range of the @@ -218,6 +257,12 @@ class GraphBasedTrackSeeder { /// @return Whether segment is within beamspot range bool checkZ0BitMask(std::uint16_t z0BitMask, float z0, float minZ0, float z0HistoCoeff) const; + + float estimateCurvature(const std::array& nodes) const; + + bool validateTriplet(const std::array candidateTriplet, + float tripletMinPt, float tauRatio, float tauRatioCut, + Options options) const; }; } // namespace Acts::Experimental diff --git a/Core/src/Seeding2/GraphBasedTrackSeeder.cpp b/Core/src/Seeding2/GraphBasedTrackSeeder.cpp index 5fb1083c192..df473cbd190 100644 --- a/Core/src/Seeding2/GraphBasedTrackSeeder.cpp +++ b/Core/src/Seeding2/GraphBasedTrackSeeder.cpp @@ -27,6 +27,11 @@ GraphBasedTrackSeeder::DerivedConfig::DerivedConfig(const Config& config) phiSliceWidth = 2 * std::numbers::pi_v / config.nMaxPhiSlice; } +GraphBasedTrackSeeder::Options::Options(float bFieldInZ_) + : bFieldInZ(bFieldInZ_) { + ptCoeff = 0.5f * bFieldInZ * Acts::UnitConstants::m; +} + GraphBasedTrackSeeder::GraphBasedTrackSeeder( const DerivedConfig& config, std::shared_ptr geometry, std::unique_ptr logger) @@ -38,7 +43,8 @@ GraphBasedTrackSeeder::GraphBasedTrackSeeder( SeedContainer2 GraphBasedTrackSeeder::createSeeds( const SpacePointContainer2& spacePoints, const GbtsRoiDescriptor& roi, - const std::uint32_t maxLayers, const GbtsTrackingFilter& filter) const { + const std::uint32_t maxLayers, const GbtsTrackingFilter& filter, + Options options) const { GbtsNodeStorage nodeStorage(m_geometry, m_mlLut); SeedContainer2 SeedContainer; @@ -75,7 +81,7 @@ SeedContainer2 GraphBasedTrackSeeder::createSeeds( std::vector edgeStorage; std::pair graphStats = - buildTheGraph(roi, nodeStorage, edgeStorage); + buildTheGraph(roi, nodeStorage, edgeStorage, options); ACTS_DEBUG("Created graph with " << graphStats.first << " edges and " << graphStats.second << " edge links"); @@ -88,29 +94,17 @@ SeedContainer2 GraphBasedTrackSeeder::createSeeds( ACTS_DEBUG("Reached Level " << maxLevel << " after GNN iterations"); - std::vector vSeedCandidates; + std::vector vOutputSeeds; extractSeedsFromTheGraph(maxLevel, graphStats.first, spacePoints.size(), - edgeStorage, vSeedCandidates, filter); + edgeStorage, vOutputSeeds, filter); - if (vSeedCandidates.empty()) { + if (vOutputSeeds.empty()) { ACTS_WARNING("No Seed Candidates"); } - - for (const auto& seed : vSeedCandidates) { - if (seed.isClone != 0) { - continue; - } - - // add to seed container: - std::vector Sp_Indexes{}; - Sp_Indexes.reserve(seed.spacePoints.size()); - - for (const auto& sp_idx : seed.spacePoints) { - Sp_Indexes.emplace_back(sp_idx); - } - + // add to seed container: + for (const auto& seed : vOutputSeeds) { auto newSeed = SeedContainer.createSeed(); - newSeed.assignSpacePointIndices(Sp_Indexes); + newSeed.assignSpacePointIndices(seed.spacePoints); newSeed.quality() = seed.seedQuality; } @@ -191,7 +185,7 @@ std::vector> GraphBasedTrackSeeder::createNodes( std::pair GraphBasedTrackSeeder::buildTheGraph( const GbtsRoiDescriptor& roi, GbtsNodeStorage& nodeStorage, - std::vector& edgeStorage) const { + std::vector& edgeStorage, Options options) const { // phi cut for triplets const float cutDPhiMax = m_cfg.lrtMode ? 0.07f : 0.012f; // curv cut for triplets @@ -214,9 +208,9 @@ std::pair GraphBasedTrackSeeder::buildTheGraph( const float tripletPtMin = 0.8f * m_cfg.minPt; // to re-scale original tunings done for the 900 MeV pT cut - const float ptScale = (0.9f * UnitConstants::GeV) / m_cfg.minPt; + const float ptScale = 0.9f / m_cfg.minPt; - const float maxCurv = m_cfg.ptCoeff / tripletPtMin; + const float maxCurv = options.ptCoeff / tripletPtMin; float maxKappaHighEta = m_cfg.lrtMode ? 1.0f * maxCurv : std::sqrt(0.8f) * maxCurv; @@ -258,6 +252,8 @@ std::pair GraphBasedTrackSeeder::buildTheGraph( const std::uint32_t layerId1 = B1.layerId; + const bool isBarrel1 = (layerId1 / 10000) == 8; + // prepare a sliding window for each bin2 in the group std::vector phiSlidingWindow; @@ -326,6 +322,10 @@ std::pair GraphBasedTrackSeeder::buildTheGraph( const GbtsEtaBin& B2 = *slw.bin; + const std::uint32_t lk2 = B2.layerId; + + const bool isBarrel2 = (lk2 / 10000) == 8; + const float deltaPhi = slw.deltaPhi; // sliding window phi1 +/- deltaPhi @@ -480,9 +480,32 @@ std::pair GraphBasedTrackSeeder::buildTheGraph( continue; } - const float tauRatio = pS->p[0] * uat2 - 1.0f; + const std::uint32_t lk3 = + m_geometry->layerIdByIndex(pS->n2->layer); - if (std::abs(tauRatio) > cutTauRatioMax) { // bad match + const bool isBarrel3 = (lk3 / 10000) == 8; + + const float absTauRatio = std::abs(pS->p[0] * uat2 - 1.0f); + float addTauRatioCorr = 0; + + if (m_cfg.useAdaptiveCuts) { + if (isBarrel1 && isBarrel2 && isBarrel3) { + const bool noGap = + ((lk3 - lk2) == 1000) && ((lk2 - layerId1) == 1000); + + // assume more scattering due to the layer in between + if (!noGap) { + addTauRatioCorr = m_cfg.tauRatioCorr; + } + } else { + bool mixedTriplet = isBarrel1 && isBarrel2 && !isBarrel3; + if (mixedTriplet) { + addTauRatioCorr = m_cfg.tauRatioCorr; + } + } + } + // bad match + if (absTauRatio > cutTauRatioMax + addTauRatioCorr) { continue; } @@ -504,6 +527,20 @@ std::pair GraphBasedTrackSeeder::buildTheGraph( continue; } + // final check: cuts on pT and d0 + if (m_cfg.validateTriplets) { + // Pixel barrel + if (isBarrel1 && isBarrel2 && isBarrel3) { + const std::array candidateTriplet = { + B1.vn[n1Idx], B2.vn[n2Idx], pS->n2}; + + if (!validateTriplet(candidateTriplet, tripletPtMin, + absTauRatio, cutTauRatioMax, options)) { + continue; + } + } + } + pS->vNei[pS->nNei] = outEdgeIdx; ++pS->nNei; @@ -560,7 +597,7 @@ std::int32_t GraphBasedTrackSeeder::runCCA( std::uint32_t iter = 0; - std::vector v_old; + std::vector vOld; for (std::uint32_t edgeIndex = 0; edgeIndex < nEdges; ++edgeIndex) { GbtsEdge* pS = &(edgeStorage[edgeIndex]); @@ -570,17 +607,17 @@ std::int32_t GraphBasedTrackSeeder::runCCA( // TODO: increment level for segments as they already have at least one // neighbour - v_old.push_back(pS); + vOld.push_back(pS); } - std::vector v_new; - v_new.reserve(v_old.size()); + std::vector vNew; + vNew.reserve(vOld.size()); // generate proposals for (; iter < maxIter; iter++) { - v_new.clear(); + vNew.clear(); - for (GbtsEdge* pS : v_old) { + for (GbtsEdge* pS : vOld) { std::int32_t nextLevel = pS->level; for (std::uint32_t nIdx = 0; nIdx < pS->nNei; ++nIdx) { @@ -590,7 +627,7 @@ std::int32_t GraphBasedTrackSeeder::runCCA( if (pS->level == pN->level) { nextLevel = pS->level + 1; - v_new.push_back(pS); + vNew.push_back(pS); break; } } @@ -603,7 +640,7 @@ std::int32_t GraphBasedTrackSeeder::runCCA( std::uint32_t nChanges = 0; - for (auto pS : v_new) { + for (auto pS : vNew) { if (pS->next != pS->level) { nChanges++; pS->level = pS->next; @@ -617,8 +654,8 @@ std::int32_t GraphBasedTrackSeeder::runCCA( break; } - v_old.swap(v_new); - v_new.clear(); + vOld.swap(vNew); + vNew.clear(); } return maxLevel; @@ -627,10 +664,8 @@ std::int32_t GraphBasedTrackSeeder::runCCA( void GraphBasedTrackSeeder::extractSeedsFromTheGraph( std::uint32_t maxLevel, std::uint32_t nEdges, std::int32_t nHits, std::vector& edgeStorage, - std::vector& vSeedCandidates, + std::vector& vOutputSeeds, const GbtsTrackingFilter& filter) const { - vSeedCandidates.clear(); - // a triplet + 2 confirmation std::uint32_t minLevel = 3; @@ -666,8 +701,16 @@ void GraphBasedTrackSeeder::extractSeedsFromTheGraph( // backtracking + std::vector vSeedCandidates; + vSeedCandidates.reserve(vSeeds.size()); + std::vector> vArgSort; + + vArgSort.reserve(vSeeds.size()); + + std::uint32_t seedCounter = 0; + GbtsTrackingFilter::State filterState{}; for (GbtsEdge* pS : vSeeds) { @@ -706,33 +749,86 @@ void GraphBasedTrackSeeder::extractSeedsFromTheGraph( continue; } - std::vector vSpIdx; - vSpIdx.resize(vN.size()); + const std::uint32_t origSeedSize = vN.size(); + + const float origSeedQuality = -rs.j / origSeedSize; + + std::uint32_t seedSplitFlag = (seedEta < m_cfg.maxSeedSplitEta) && + (origSeedSize > 3) && + (origSeedSize <= 5) + ? 1 + : 0; + + // split the seed by dropping spacepoints + if (seedSplitFlag != 0) { + // 2. "drop-outs" and the original seed candidate + std::array, 3> triplets{}; + + // triplet parameter estimate + std::array invRads{}; + + triplets[0] = {vN[0], vN[origSeedSize / 2], vN[origSeedSize - 1]}; + + // all but the first one + const std::vector dropOut1(vN.begin() + 1, vN.end()); + + triplets[1] = {dropOut1[0], dropOut1[(origSeedSize - 1) / 2], + dropOut1[origSeedSize - 2]}; + + std::vector dopOut2; - for (std::uint32_t k = 0; k < vN.size(); ++k) { - vSpIdx[k] = vN[k]->idx; + dopOut2.reserve(origSeedSize - 1); + + for (std::uint32_t k = 0; k < origSeedSize; k++) { + if (k == origSeedSize / 2) { + continue; // drop the middle SP in the original seed + } + + dopOut2.emplace_back(vN[k]); + } + + triplets[2] = {dopOut2[0], dopOut2[(origSeedSize - 1) / 2], + dopOut2[origSeedSize - 2]}; + + for (std::uint32_t k = 0; k < invRads.size(); k++) { + invRads[k] = estimateCurvature(triplets[k]); + } + + const std::array diffs = {std::abs(invRads[1] - invRads[0]), + std::abs(invRads[2] - invRads[0]), + std::abs(invRads[2] - invRads[1])}; + + const bool confirmed = diffs[0] < m_cfg.maxInvRadDiff && + diffs[1] < m_cfg.maxInvRadDiff && + diffs[2] < m_cfg.maxInvRadDiff; + + if (confirmed) { + seedSplitFlag = 0; // reset the flag + } } - vSeedCandidates.emplace_back(-rs.j / vN.size(), 0, std::move(vSpIdx)); + vSeedCandidates.emplace_back(origSeedQuality, 0, vN, seedSplitFlag); + + vArgSort.emplace_back(origSeedQuality, seedCounter); + + ++seedCounter; } // clone removal code goes below ... - std::ranges::sort(vSeedCandidates, - [](const SeedProperties& s1, const SeedProperties& s2) { - return s1 < s2; - }); + std::ranges::sort(vArgSort); std::vector h2t(nHits + 1, 0); // hit to track associations std::uint32_t trackId = 0; - for (const auto& seed : vSeedCandidates) { + for (const auto& args : vArgSort) { + const auto& seed = vSeedCandidates[args.second]; ++trackId; // loop over space points indices for (const auto& h : seed.spacePoints) { - const std::uint32_t hitId = h + 1; + const std::uint32_t hitId = h->idx + 1; const std::uint32_t tid = h2t[hitId]; @@ -744,15 +840,21 @@ void GraphBasedTrackSeeder::extractSeedsFromTheGraph( } } - for (std::uint32_t trackIdx = 0; trackIdx < vSeedCandidates.size(); - ++trackIdx) { - const std::uint32_t nTotal = vSeedCandidates[trackIdx].spacePoints.size(); + std::uint32_t trackIdx = 0; + + for (const auto& args : vArgSort) { + const auto& seed = vSeedCandidates[args.second].spacePoints; + + const std::uint32_t nTotal = seed.size(); + std::uint32_t nOther = 0; trackId = trackIdx + 1; - for (const auto& h : vSeedCandidates[trackIdx].spacePoints) { - const std::uint32_t hitId = h + 1; + ++trackIdx; + + for (const auto& h : seed) { + const std::uint32_t hitId = h->idx + 1; const std::uint32_t tid = h2t[hitId]; @@ -764,7 +866,59 @@ void GraphBasedTrackSeeder::extractSeedsFromTheGraph( if (nOther > m_cfg.hitShareThreshold * nTotal) { // reject - vSeedCandidates[trackIdx].isClone = -1; + vSeedCandidates[args.second].isClone = -1; // reject + } + } + vOutputSeeds.reserve(vSeedCandidates.size()); + + // drop the clones and split seeds if need be + + for (const auto& args : vArgSort) { + const auto& seed = vSeedCandidates[args.second]; + + if (seed.isClone != 0) { + continue; // identified as a clone of a better candidate + } + + const auto& vN = seed.spacePoints; + + if (seed.forSeedSplitting == 0) { + // add seed to output + + std::vector vSpIdx; + + vSpIdx.resize(vN.size()); + + for (std::uint32_t k = 0; k < vSpIdx.size(); k++) { + vSpIdx[k] = vN[k]->idx; + } + + vOutputSeeds.emplace_back(seed.seedQuality, vSpIdx); + + continue; + } + + // seed split into "drop-out" seeds + + const std::uint32_t seedSize = vN.size(); + + const std::array indices2drop = { + 0, seedSize / 2ul}; // the first and the middle + + for (const auto& skipIdx : indices2drop) { + std::vector newSeed; + + newSeed.reserve(seedSize - 1); + + for (std::uint32_t k = 0; k < seedSize; k++) { + if (k == skipIdx) { + continue; + } + + newSeed.emplace_back(vN[k]->idx); + } + + vOutputSeeds.emplace_back(seed.seedQuality, newSeed); } } } @@ -810,4 +964,122 @@ bool GraphBasedTrackSeeder::checkZ0BitMask(const std::uint16_t z0BitMask, return false; } +float GraphBasedTrackSeeder::estimateCurvature( + const std::array& nodes) const { + // conformal mapping with the center at the last spacepoint + + std::array u{}; + std::array v{}; + + const float x0 = nodes[2]->x; + const float y0 = nodes[2]->y; + + const float r0 = nodes[2]->r; + + const float cosA = x0 / r0; + + const float sinA = y0 / r0; + + for (std::uint32_t k = 0; k < 2; k++) { + const float dx = nodes[k]->x - x0; + + const float dy = nodes[k]->y - y0; + + const float r2Inv = 1.0 / (dx * dx + dy * dy); + + const float xn = dx * cosA + dy * sinA; + + const float yn = -dx * sinA + dy * cosA; + + u[k] = xn * r2Inv; + v[k] = yn * r2Inv; + } + + const float du = u[0] - u[1]; + + if (du == 0.0) { + return 0.0; + } + + const float A = (v[0] - v[1]) / du; + + const float B = v[1] - A * u[1]; + + // curavture in units of 1/mm + return B / std::sqrt(1 + A * A); +} + +bool GraphBasedTrackSeeder::validateTriplet( + const std::array candidateTriplet, + const float tripletMinPt, const float tauRatio, const float tauRatioCut, + Options options) const { + // conformal mapping with the center at the middle spacepoint + + std::array u{}; + std::array v{}; + + const float x0 = candidateTriplet[1]->x; + const float y0 = candidateTriplet[1]->y; + + const float r0 = candidateTriplet[1]->r; + + const float cosA = x0 / r0; + + const float sinA = y0 / r0; + + for (std::uint32_t k = 0; k < 2; k++) { + std::uint32_t spIdx = (k == 1) ? 2 : k; + + const float dx = candidateTriplet[spIdx]->x - x0; + + const float dy = candidateTriplet[spIdx]->y - y0; + + const float r2Inv = 1.0f / (dx * dx + dy * dy); + + const float xn = dx * cosA + dy * sinA; + + const float yn = -dx * sinA + dy * cosA; + + u[k] = xn * r2Inv; + v[k] = yn * r2Inv; + } + + const float du = u[0] - u[1]; + + if (du == 0.0) { + return false; + } + + const float A = (v[0] - v[1]) / du; + + const float B = v[1] - A * u[1]; + + const float d0 = r0 * (B * r0 - A); + + if (std::abs(d0) > m_cfg.d0Max) { + return false; + } + + if (B != 0.0) { // straight-line track is OK + + const float R = std::sqrt(1 + A * A) / B; + + // 1T magnetic field used + const float pT = std::abs(options.bFieldInZ * R / 2); + + if (pT < tripletMinPt) { + return false; + } + + if (pT > 5 * tripletMinPt) { // relatively high-pT track + + if (tauRatio > 0.9f * tauRatioCut) { + return false; + } + } + } + + return true; +} + } // namespace Acts::Experimental