From f8dbeeaed420bb9f960ccbda01d2a235ea740fbd Mon Sep 17 00:00:00 2001 From: decianm Date: Sat, 26 Apr 2025 18:48:55 +0200 Subject: [PATCH 01/13] try to understand differences --- Pr/PrAlgorithms/src/PrDownTrack.h | 8 ++ Pr/PrAlgorithms/src/PrLongLivedTracking.cpp | 99 +++++++++++++++++++-- 2 files changed, 100 insertions(+), 7 deletions(-) diff --git a/Pr/PrAlgorithms/src/PrDownTrack.h b/Pr/PrAlgorithms/src/PrDownTrack.h index 67a5b180bac..5e99b0e7d7f 100644 --- a/Pr/PrAlgorithms/src/PrDownTrack.h +++ b/Pr/PrAlgorithms/src/PrDownTrack.h @@ -88,6 +88,14 @@ namespace Downstream { : // lhs.lhcbID() < rhs.lhcbID() ); }; + + inline constexpr auto IncreaseByAbsProj = []( const Hit& lhs, const Hit& rhs ) { + return ( std::abs(lhs.projection) < std::abs(rhs.projection) ? true : // + std::abs(rhs.projection) < std::abs(lhs.projection) ? false + : // + lhs.lhcbID() < rhs.lhcbID() ); + }; + } // namespace Downstream /** @class PrDownTrack PrDownTrack.h diff --git a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp index 3361d140b2d..70dc70479de 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -582,6 +582,82 @@ private: } } while ( again ); } + //========================================================================= + // Fit and remove the worst hit, as long as over tolerance + // Perform a chi2 fit to the track and remove outliers + //========================================================================= + template + void simplyFit( PrDownTrack& track ) const { + + //== Fit, using the magnet point as constraint. + float mat[6], rhs[3]; + mat[0] = track.weightXMag(); + mat[1] = 0.; + mat[2] = 0.; + mat[3] = 0.; + mat[4] = 0.; + mat[5] = 0.; + rhs[0] = track.dxMagnet() * track.weightXMag(); + rhs[1] = 0.; + rhs[2] = 0.; + + //for ( std::int8_t iHit = 0; iHit < nHits; ++iHit ) { + for ( const auto& hit: track.hits()) { + const float dz = 0.001 * ( hit.z - track.zMagnet() ); + const float w = hit.weight(); + const float t = hit.sin(); + + mat[0] += w; + mat[1] += w * dz; + mat[2] += w * dz * dz; + mat[3] += w * t; + mat[4] += w * dz * t; + mat[5] += w * t * t; + + if constexpr(useProjections){ + const float dist = hit.projection; + rhs[0] += w * dist; + rhs[1] += w * dist * dz; + rhs[2] += w * dist * t; + }else{ + const float dist = track.distance( hit ); + rhs[0] += w * dist; + rhs[1] += w * dist * dz; + rhs[2] += w * dist * t; + } + + } + + // -- solve the equation and update the parameters of the track + CholeskyDecomp decomp( mat ); + if ( !decomp ) { + track.setChi2( 1e42 ); + return; + } else { + decomp.Solve( rhs ); + } + + const float dx = rhs[0]; + const float dsl = 0.001 * rhs[1]; + const float dy = rhs[2]; + + track.updateX( dx, dsl ); + track.updateY( dy ); + + //== Remove worst hit and retry, if too far. + float chi2 = track.initialChi2(); + + for ( const auto& hit: track.hits()) { + const float dist = track.distance( hit ); + chi2 += dist * dist * hit.weight(); + // -- only flag this hit as removable if it is not alone in a plane or there are 4 planes that fired + } + + if( track.hits().size() > 2 ) chi2 /= ( track.hits().size() - 2 ); + track.setChi2( chi2 ); + + } + //========================================================================= // Collect the hits in the other x layer @@ -645,7 +721,7 @@ private: //========================================================================= // Add the U hits. //========================================================================= - void addUHits( const PrDownTrack& track, Downstream::Hits& preSelHits, Downstream::Hits& uHitsTemp, + void addUHits( PrDownTrack& track, const Downstream::Hits& preSelHits, Downstream::Hits& uHitsTemp, PrDownTracks& goodXUTracks ) const { goodXUTracks.clear(); @@ -660,20 +736,29 @@ private: // -- first select all hits, and then // -- accept until over a tolerance - for ( auto& hit : preSelHits ) { - const float dist = std::abs( track.distance( hit ) ); - if ( dist > tol ) continue; - hit.projection = dist; + for ( const auto& hit : preSelHits ) { + const float dist = track.distance( hit ); + if ( std::abs(dist) > tol ) continue; + //hit.projection = dist; uHitsTemp.push_back( hit ); + uHitsTemp.back().projection = dist; + } + + std::sort( uHitsTemp.begin(), uHitsTemp.end(), Downstream::IncreaseByAbsProj ); + + for(auto& hit : track.hits()){ + hit.projection = track.distance( hit ); } - std::sort( uHitsTemp.begin(), uHitsTemp.end(), Downstream::IncreaseByProj ); for ( auto& hit : uHitsTemp ) { auto& greatTrack = goodXUTracks.emplace_back( track ); + + greatTrack.hits().push_back( hit ); - fitAndRemove( greatTrack ); + //fitAndRemove( greatTrack ); + simplyFit( greatTrack ); if ( !greatTrack.isYCompatible( m_yTol ) ) { goodXUTracks.pop_back(); -- GitLab From 783e5b844e78c9435b6984a67482f641ee30048e Mon Sep 17 00:00:00 2001 From: decianm Date: Sun, 27 Apr 2025 13:44:14 +0200 Subject: [PATCH 02/13] clean up to avoid different number of tracks --- Pr/PrAlgorithms/src/PrDownTrack.h | 10 +-- Pr/PrAlgorithms/src/PrLongLivedTracking.cpp | 78 +++++++++++---------- 2 files changed, 46 insertions(+), 42 deletions(-) diff --git a/Pr/PrAlgorithms/src/PrDownTrack.h b/Pr/PrAlgorithms/src/PrDownTrack.h index 5e99b0e7d7f..dac6c40f420 100644 --- a/Pr/PrAlgorithms/src/PrDownTrack.h +++ b/Pr/PrAlgorithms/src/PrDownTrack.h @@ -74,7 +74,7 @@ namespace Downstream { std::max( mH.get().cast(), mH.get().cast() ); return ( ( ( yMin - tol ) <= y ) && ( y <= ( yMax + tol ) ) ); } - [[nodiscard]] auto xAt( const float y ) const { + [[nodiscard]] inline auto xAt( const float y ) const { const auto mH = getScalarHit(); return mH.get().cast() + y * mH.get().cast(); } @@ -115,7 +115,7 @@ namespace Downstream { class PrDownTrack final { public: - using Hits = boost::container::small_vector>; + using Hits = boost::container::small_vector>; // using Hits = boost::container::static_vector; // Until we can put a bound on the number of hits, use a small_vector @@ -194,12 +194,12 @@ public: void setChi2( float chi2 ) noexcept { m_chi2 = chi2; } // functions - float xAtZ( float z ) const noexcept { + inline float xAtZ( float z ) const noexcept { const float curvature = 1.6e-5 * ( stateTx() - m_slopeX ); return xMagnet() + ( z - zMagnet() ) * m_slopeX + curvature * ( z - m_zUT ) * ( z - m_zUT ); } - float yAtZ( float z ) const noexcept { return yMagnet() + m_displY + ( z - zMagnet() ) * slopeY(); } + inline float yAtZ( float z ) const noexcept { return yMagnet() + m_displY + ( z - zMagnet() ) * slopeY(); } void updateX( float dx, float dsl ) noexcept { m_displX += dx; @@ -221,7 +221,7 @@ public: return sinTrack * std::abs( momentum() ); } - float distance( const Downstream::Hit& hit ) const noexcept { return hit.xAt( yAtZ( hit.z ) ) - xAtZ( hit.z ); } + inline float distance( const Downstream::Hit& hit ) const noexcept { return hit.xAt( yAtZ( hit.z ) ) - xAtZ( hit.z ); } bool isYCompatible( const float tol ) const noexcept { return std ::all_of( m_hits.begin(), m_hits.end(), diff --git a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp index 70dc70479de..9621f86f55d 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -373,13 +373,16 @@ private: // -- Take all xTracks into account whose chi2 is close to the best // if ( xTrack.chi2() - m_maxChi2DistXTracks >= goodXTracks[0].chi2() ) break; + + addUHits( xTrack, preSelHits[plane_indx[2]], uHitsTemp, goodXUTracks ); // -- Loop over good xu tracks for ( PrDownTrack& xuTrack : goodXUTracks ) { addVHits( xuTrack, preSelHits[plane_indx[3]] ); if ( xuTrack.hits().size() < 3 ) continue; - fitAndRemove( xuTrack ); + simplyFit( xuTrack ); + //fitAndRemove( xuTrack ); if ( xuTrack.chi2() > m_maxChi2 || !xuTrack.isYCompatible( m_yTol ) ) continue; trackCandidates.push_back( std::move( xuTrack ) ); } // Loop over good xu tracks @@ -461,8 +464,8 @@ private: if ( xAtYEq0 < lowerBoundX || !( yMin - yTol <= yPredLay && yPredLay <= yMax + yTol ) ) continue; auto xx = mH.get().cast() + y * dxDy; - if ( xPredTol < pos - xx ) continue; // go from -x to +x - if ( xPredTol < xx - pos ) continue; // can break if we go out of the right bound + if ( xPredTol < std::abs(pos - xx) ) continue; // go from -x to +x + if ( xPredTol < xx - pos ) break; // can break if we go out of the right bound preSelHits[layer].emplace_back( &hitHandler, itHit, xx, zLayer, fabs( xx - pos ) ); } @@ -478,7 +481,6 @@ private: // Fit and remove the worst hit, as long as over tolerance // Perform a chi2 fit to the track and remove outliers //========================================================================= - template void fitAndRemove( PrDownTrack& track ) const { if ( track.hits().size() < 2 ) { return; // no fit if single point only ! @@ -499,7 +501,7 @@ private: rhs[1] = 0.; rhs[2] = 0.; - std::array differentPlanes = { 0, 0, 0, 0 }; + std::array differentPlanes = { 0, 0, 0, 0 }; for ( auto& hit : track.hits() ) { const float dz = 0.001 * ( hit.z - track.zMagnet() ); @@ -521,9 +523,9 @@ private: differentPlanes[hit.planeCode()]++; } - const unsigned int nDoF = + const std::uint8_t nDoF = std::count_if( differentPlanes.begin(), differentPlanes.end(), []( const int a ) { return a > 0; } ); - const int nbUV = differentPlanes[1] + differentPlanes[2]; + const std::uint8_t nbUV = differentPlanes[1] + differentPlanes[2]; // -- solve the equation and update the parameters of the track CholeskyDecomp decomp( mat ); @@ -550,20 +552,19 @@ private: for ( auto itH = track.hits().begin(); itH != track.hits().end(); ++itH ) { Downstream::Hit& hit = *itH; - if ( !onlyFit ) { - const float yTrackAtZ = track.yAtZ( hit.z ); - if ( !hit.isYCompatible( yTrackAtZ, m_yTol ) ) { - track.hits().erase( itH ); - if ( 2 < track.hits().size() ) again = true; - break; - } + const float yTrackAtZ = track.yAtZ( hit.z ); + if ( !hit.isYCompatible( yTrackAtZ, m_yTol ) ) { + track.hits().erase( itH ); + if ( 2 < track.hits().size() ) again = true; + break; } + const float dist = std::abs( track.distance( hit ) ); hit.projection = dist; chi2 += dist * dist * hit.weight(); // -- only flag this hit as removable if it is not alone in a plane or there are 4 planes that fired - if ( !onlyFit && maxDist < dist && ( 1 < differentPlanes[hit.planeCode()] || nDoF == track.hits().size() ) ) { + if ( maxDist < dist && ( 1 < differentPlanes[hit.planeCode()] || nDoF == track.hits().size() ) ) { maxDist = dist; worst = itH; } @@ -574,8 +575,6 @@ private: if ( track.hits().size() > 2 ) chi2 /= ( track.hits().size() - 2 ); track.setChi2( chi2 ); - if ( onlyFit ) { return; } - if ( m_maxChi2 < chi2 && track.hits().size() > 3 && maxDist > 0 ) { track.hits().erase( worst ); again = true; @@ -644,11 +643,11 @@ private: track.updateX( dx, dsl ); track.updateY( dy ); - //== Remove worst hit and retry, if too far. float chi2 = track.initialChi2(); - for ( const auto& hit: track.hits()) { + for ( auto& hit: track.hits()) { const float dist = track.distance( hit ); + hit.projection = dist;// this is needed later on for fitting with projections chi2 += dist * dist * hit.weight(); // -- only flag this hit as removable if it is not alone in a plane or there are 4 planes that fired } @@ -750,8 +749,7 @@ private: hit.projection = track.distance( hit ); } - - for ( auto& hit : uHitsTemp ) { + for ( const auto& hit : uHitsTemp ) { auto& greatTrack = goodXUTracks.emplace_back( track ); @@ -760,16 +758,20 @@ private: //fitAndRemove( greatTrack ); simplyFit( greatTrack ); - if ( !greatTrack.isYCompatible( m_yTol ) ) { - goodXUTracks.pop_back(); - continue; - } + // -- it's sorted if ( greatTrack.chi2() > minChi2 ) { goodXUTracks.pop_back(); break; } + + if ( !greatTrack.isYCompatible( m_yTol ) ) { + goodXUTracks.pop_back(); + continue; + } + + if ( goodXUTracks.size() >= m_maxXUTracks ) { break; } } @@ -780,23 +782,23 @@ private: //========================================================================= // Add the V hits. Take the one which has the best chi2 //========================================================================= - void addVHits( PrDownTrack& track, Downstream::Hits& preSelHits ) const { + void addVHits( PrDownTrack& track, const Downstream::Hits& preSelHits ) const { if ( preSelHits.empty() ) { return; } - auto p = std::abs( track.momentum() ); - float tol = ( track.hits().size() == 2 ) ? m_tolUOffset + m_tolUConst / p : m_tolVOffset + m_tolVConst / p; + const auto p = std::abs( track.momentum() ); + const float tol = ( track.hits().size() == 2 ) ? m_tolUOffset + m_tolUConst / p : m_tolVOffset + m_tolVConst / p; // Find the closest hit to the track // It's enough to look at the closest hit only, since absolute majority of the tracks // with not a closest hit get killed later anyways // This also doesn't affect efficiency, as the addOverlapRegions function probes all close enough hits as well float bestDist = 1e9f; - Downstream::Hit* bestHit = nullptr; + const Downstream::Hit* bestHit = nullptr; for ( auto& hit : preSelHits ) { - const float dist = std::abs( track.distance( hit ) ); - if ( dist > tol ) continue; - hit.projection = dist; - if ( dist < bestDist ) { + const float dist = track.distance( hit ); + if ( std::abs(dist) > tol ) continue; + //hit.projection = dist; + if ( std::abs(dist) < std::abs(bestDist) ) { bestDist = dist; bestHit = &hit; } @@ -805,7 +807,9 @@ private: // Fit the track if ( bestHit ) { track.hits().push_back( *bestHit ); - fitAndRemove( track ); + track.hits().back().projection = bestDist; + simplyFit( track ); + //fitAndRemove( track ); } track.sortFinalHits(); } @@ -815,7 +819,7 @@ private: //============================================================================= // Fit the projection in the zx plane, one hit in each x layer //============================================================================= - void fitXProjection( const PrDownTrack& track, Downstream::Hit& firstHit, Downstream::Hits& matchingXHits, + void fitXProjection( const PrDownTrack& track, const Downstream::Hit& firstHit, const Downstream::Hits& matchingXHits, PrDownTracks& goodXTracks, const bool skip_second_layer ) const { goodXTracks.clear(); @@ -879,7 +883,7 @@ private: // Maybe we could make this smarter and do it for every track and add the 'second best' // this, such that we do not need to loop over them again //============================================================================= - void addOverlapRegions( PrDownTrack& track, std::array& preSelHits ) const { + void addOverlapRegions( PrDownTrack& track, const std::array& preSelHits ) const { bool hitAdded = false; const int trackNbHits = track.hits().size(); @@ -905,7 +909,7 @@ private: if ( hitAdded ) { track.sortFinalHits(); - fitAndRemove( track ); + fitAndRemove( track ); } } -- GitLab From e603009cd5183552a11dc181488be450df147d17 Mon Sep 17 00:00:00 2001 From: decianm Date: Wed, 14 May 2025 15:17:53 +0200 Subject: [PATCH 03/13] Use floats throughout in PrDownTrack, move easier to compute statements earlier in OverlapRegions, check for beampipe earlier, write simplified distance function --- Pr/PrAlgorithms/src/PrDownTrack.h | 29 +++++++++++++---- Pr/PrAlgorithms/src/PrLongLivedTracking.cpp | 35 ++++++++++++++------- 2 files changed, 47 insertions(+), 17 deletions(-) diff --git a/Pr/PrAlgorithms/src/PrDownTrack.h b/Pr/PrAlgorithms/src/PrDownTrack.h index dac6c40f420..20fd190a6d6 100644 --- a/Pr/PrAlgorithms/src/PrDownTrack.h +++ b/Pr/PrAlgorithms/src/PrDownTrack.h @@ -121,16 +121,23 @@ public: PrDownTrack( Gaudi::TrackVector stateVector, float stateZ, float zUT, LHCb::span magnetParams, LHCb::span yParams, LHCb::span momPar, float magnetScale ) - : m_stateVector( stateVector ), m_stateZ( stateZ ), m_zUT( zUT ) { + : m_stateZ( stateZ ), m_zUT( zUT ) { + + m_stateVector[0] = (float)stateVector[0]; + m_stateVector[1] = (float)stateVector[1]; + m_stateVector[2] = (float)stateVector[2]; + m_stateVector[3] = (float)stateVector[3]; + m_stateVector[4] = (float)stateVector[4]; + const auto tx2 = stateTx() * stateTx(); const auto ty2 = stateTy() * stateTy(); m_momentumParam = ( momPar[0] + momPar[1] * tx2 + momPar[2] * ty2 ) * magnetScale; // -- See PrFitKsParams to see how these coefficients are derived. const float zMagnet = magnetParams[0] + magnetParams[1] * ty2 + magnetParams[2] * tx2 + - magnetParams[3] * std::abs( stateQoP() ) + /// this is where the old one stopped. - magnetParams[4] * std::abs( stateX() ) + magnetParams[5] * std::abs( stateY() ) + - magnetParams[6] * std::abs( stateTy() ); + magnetParams[3] * std::abs( stateQoP() ) + /// this is where the old one stopped. + magnetParams[4] * std::abs( stateX() ) + magnetParams[5] * std::abs( stateY() ) + + magnetParams[6] * std::abs( stateTy() ); const float dz = zMagnet - stateZ; const float xMagnet = stateX() + dz * stateTx(); @@ -199,6 +206,12 @@ public: return xMagnet() + ( z - zMagnet() ) * m_slopeX + curvature * ( z - m_zUT ) * ( z - m_zUT ); } + // -- this ignores the change in curvature, to be used for small distances + inline float xAtZFast( float z, float xStart, float zStart ) const noexcept { + return xStart + ( z - zStart ) * m_slopeX; + } + + inline float yAtZ( float z ) const noexcept { return yMagnet() + m_displY + ( z - zMagnet() ) * slopeY(); } void updateX( float dx, float dsl ) noexcept { @@ -223,6 +236,8 @@ public: inline float distance( const Downstream::Hit& hit ) const noexcept { return hit.xAt( yAtZ( hit.z ) ) - xAtZ( hit.z ); } + inline float distanceFast( const Downstream::Hit& hit, float zStart, float xStart ) const noexcept { return hit.xAt( yAtZ( hit.z ) ) - xAtZFast( hit.z, zStart, xStart ); } + bool isYCompatible( const float tol ) const noexcept { return std ::all_of( m_hits.begin(), m_hits.end(), [&]( const auto& hit ) { return hit.isYCompatible( yAtZ( hit.z ), tol ); } ); @@ -235,9 +250,11 @@ public: } private: - Gaudi::TrackVector m_stateVector; + + // -- using the float versions everywhere + Gaudi::TrackVectorF m_stateVector; float m_stateZ{ 0 }; - Gaudi::XYZPoint m_magnet; + Gaudi::XYZPointF m_magnet; float m_momentumParam{ 0 }; float m_zUT{ 0 }; diff --git a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp index 9621f86f55d..3ecd841d446 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -197,6 +197,8 @@ public: for ( PrDownTrack& track : trackCandidates ) { if ( has4LayerTrack && track.hits().size() < 4 ) continue; + + addOverlapRegions( track, preSelHits ); if ( track.chi2() > m_maxChi2 || insideBeampipe( track ) || track.hits().size() < m_MinNumUTHits || @@ -373,8 +375,6 @@ private: // -- Take all xTracks into account whose chi2 is close to the best // if ( xTrack.chi2() - m_maxChi2DistXTracks >= goodXTracks[0].chi2() ) break; - - addUHits( xTrack, preSelHits[plane_indx[2]], uHitsTemp, goodXUTracks ); // -- Loop over good xu tracks @@ -383,7 +383,7 @@ private: if ( xuTrack.hits().size() < 3 ) continue; simplyFit( xuTrack ); //fitAndRemove( xuTrack ); - if ( xuTrack.chi2() > m_maxChi2 || !xuTrack.isYCompatible( m_yTol ) ) continue; + if ( xuTrack.chi2() > m_maxChi2 || !xuTrack.isYCompatible( m_yTol ) || insideBeampipe( xuTrack ) ) continue; trackCandidates.push_back( std::move( xuTrack ) ); } // Loop over good xu tracks } // Loop over good x tracks @@ -735,8 +735,12 @@ private: // -- first select all hits, and then // -- accept until over a tolerance + const float xStart = track.xAtZ(preSelHits[0].z); + const float zStart = preSelHits[0].z; + for ( const auto& hit : preSelHits ) { - const float dist = track.distance( hit ); + //const float dist = track.distance( hit ); + const float dist = track.distanceFast( hit, xStart, zStart ); if ( std::abs(dist) > tol ) continue; //hit.projection = dist; uHitsTemp.push_back( hit ); @@ -794,8 +798,13 @@ private: // This also doesn't affect efficiency, as the addOverlapRegions function probes all close enough hits as well float bestDist = 1e9f; const Downstream::Hit* bestHit = nullptr; + + const float xStart = track.xAtZ(preSelHits[0].z); + const float zStart = preSelHits[0].z; + for ( auto& hit : preSelHits ) { - const float dist = track.distance( hit ); + //const float dist = track.distance( hit ); + const float dist = track.distanceFast( hit, xStart, zStart ); if ( std::abs(dist) > tol ) continue; //hit.projection = dist; if ( std::abs(dist) < std::abs(bestDist) ) { @@ -834,9 +843,14 @@ private: goodXTracks.pop_back(); continue; } + + + tr.hits().push_back( firstHit ); tr.hits().push_back( hit ); + + if ( goodXTracks.size() >= 3 ) { break; } } @@ -893,16 +907,15 @@ private: if ( j++ >= trackNbHits ) break; for ( auto& hit : preSelHits[trackHit.planeCode()] ) { + // -- the displacement in z between overlap modules is larger than 1mm + if ( std::abs( hit.z - trackHit.z ) < 1.0 ) continue; if ( m_overlapTol > std::abs( track.distance( hit ) ) ) { float yTrack = track.yAtZ( hit.z ); if ( !hit.isYCompatible( yTrack, m_yTol ) ) continue; - - // -- the displacement in z between overlap modules is larger than 1mm - if ( std::abs( hit.z - trackHit.z ) >= 1.0 ) { - track.hits().push_back( hit ); - hitAdded = true; - } + + track.hits().push_back( hit ); + hitAdded = true; } } } -- GitLab From 87a98b66466bb1c9136b9dda6f423bdfbe1d96cd Mon Sep 17 00:00:00 2001 From: Gitlab CI Date: Wed, 14 May 2025 13:36:23 +0000 Subject: [PATCH 04/13] pre-commit fixes patch generated by https://gitlab.cern.ch/lhcb/Rec/-/jobs/55792616 --- Pr/PrAlgorithms/src/PrDownTrack.h | 26 +++---- Pr/PrAlgorithms/src/PrLongLivedTracking.cpp | 79 ++++++++------------- 2 files changed, 45 insertions(+), 60 deletions(-) diff --git a/Pr/PrAlgorithms/src/PrDownTrack.h b/Pr/PrAlgorithms/src/PrDownTrack.h index 20fd190a6d6..55fb8524d23 100644 --- a/Pr/PrAlgorithms/src/PrDownTrack.h +++ b/Pr/PrAlgorithms/src/PrDownTrack.h @@ -90,9 +90,9 @@ namespace Downstream { }; inline constexpr auto IncreaseByAbsProj = []( const Hit& lhs, const Hit& rhs ) { - return ( std::abs(lhs.projection) < std::abs(rhs.projection) ? true : // - std::abs(rhs.projection) < std::abs(lhs.projection) ? false - : // + return ( std::abs( lhs.projection ) < std::abs( rhs.projection ) ? true : // + std::abs( rhs.projection ) < std::abs( lhs.projection ) ? false + : // lhs.lhcbID() < rhs.lhcbID() ); }; @@ -121,7 +121,7 @@ public: PrDownTrack( Gaudi::TrackVector stateVector, float stateZ, float zUT, LHCb::span magnetParams, LHCb::span yParams, LHCb::span momPar, float magnetScale ) - : m_stateZ( stateZ ), m_zUT( zUT ) { + : m_stateZ( stateZ ), m_zUT( zUT ) { m_stateVector[0] = (float)stateVector[0]; m_stateVector[1] = (float)stateVector[1]; @@ -135,9 +135,9 @@ public: // -- See PrFitKsParams to see how these coefficients are derived. const float zMagnet = magnetParams[0] + magnetParams[1] * ty2 + magnetParams[2] * tx2 + - magnetParams[3] * std::abs( stateQoP() ) + /// this is where the old one stopped. - magnetParams[4] * std::abs( stateX() ) + magnetParams[5] * std::abs( stateY() ) + - magnetParams[6] * std::abs( stateTy() ); + magnetParams[3] * std::abs( stateQoP() ) + /// this is where the old one stopped. + magnetParams[4] * std::abs( stateX() ) + magnetParams[5] * std::abs( stateY() ) + + magnetParams[6] * std::abs( stateTy() ); const float dz = zMagnet - stateZ; const float xMagnet = stateX() + dz * stateTx(); @@ -211,7 +211,6 @@ public: return xStart + ( z - zStart ) * m_slopeX; } - inline float yAtZ( float z ) const noexcept { return yMagnet() + m_displY + ( z - zMagnet() ) * slopeY(); } void updateX( float dx, float dsl ) noexcept { @@ -234,9 +233,13 @@ public: return sinTrack * std::abs( momentum() ); } - inline float distance( const Downstream::Hit& hit ) const noexcept { return hit.xAt( yAtZ( hit.z ) ) - xAtZ( hit.z ); } + inline float distance( const Downstream::Hit& hit ) const noexcept { + return hit.xAt( yAtZ( hit.z ) ) - xAtZ( hit.z ); + } - inline float distanceFast( const Downstream::Hit& hit, float zStart, float xStart ) const noexcept { return hit.xAt( yAtZ( hit.z ) ) - xAtZFast( hit.z, zStart, xStart ); } + inline float distanceFast( const Downstream::Hit& hit, float zStart, float xStart ) const noexcept { + return hit.xAt( yAtZ( hit.z ) ) - xAtZFast( hit.z, zStart, xStart ); + } bool isYCompatible( const float tol ) const noexcept { return std ::all_of( m_hits.begin(), m_hits.end(), @@ -250,10 +253,9 @@ public: } private: - // -- using the float versions everywhere Gaudi::TrackVectorF m_stateVector; - float m_stateZ{ 0 }; + float m_stateZ{ 0 }; Gaudi::XYZPointF m_magnet; float m_momentumParam{ 0 }; diff --git a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp index 3ecd841d446..3e256ce6a98 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -197,8 +197,6 @@ public: for ( PrDownTrack& track : trackCandidates ) { if ( has4LayerTrack && track.hits().size() < 4 ) continue; - - addOverlapRegions( track, preSelHits ); if ( track.chi2() > m_maxChi2 || insideBeampipe( track ) || track.hits().size() < m_MinNumUTHits || @@ -382,8 +380,8 @@ private: addVHits( xuTrack, preSelHits[plane_indx[3]] ); if ( xuTrack.hits().size() < 3 ) continue; simplyFit( xuTrack ); - //fitAndRemove( xuTrack ); - if ( xuTrack.chi2() > m_maxChi2 || !xuTrack.isYCompatible( m_yTol ) || insideBeampipe( xuTrack ) ) continue; + // fitAndRemove( xuTrack ); + if ( xuTrack.chi2() > m_maxChi2 || !xuTrack.isYCompatible( m_yTol ) || insideBeampipe( xuTrack ) ) continue; trackCandidates.push_back( std::move( xuTrack ) ); } // Loop over good xu tracks } // Loop over good x tracks @@ -464,8 +462,8 @@ private: if ( xAtYEq0 < lowerBoundX || !( yMin - yTol <= yPredLay && yPredLay <= yMax + yTol ) ) continue; auto xx = mH.get().cast() + y * dxDy; - if ( xPredTol < std::abs(pos - xx) ) continue; // go from -x to +x - if ( xPredTol < xx - pos ) break; // can break if we go out of the right bound + if ( xPredTol < std::abs( pos - xx ) ) continue; // go from -x to +x + if ( xPredTol < xx - pos ) break; // can break if we go out of the right bound preSelHits[layer].emplace_back( &hitHandler, itHit, xx, zLayer, fabs( xx - pos ) ); } @@ -558,7 +556,6 @@ private: if ( 2 < track.hits().size() ) again = true; break; } - const float dist = std::abs( track.distance( hit ) ); hit.projection = dist; @@ -585,9 +582,9 @@ private: // Fit and remove the worst hit, as long as over tolerance // Perform a chi2 fit to the track and remove outliers //========================================================================= - template + template void simplyFit( PrDownTrack& track ) const { - + //== Fit, using the magnet point as constraint. float mat[6], rhs[3]; mat[0] = track.weightXMag(); @@ -600,11 +597,11 @@ private: rhs[1] = 0.; rhs[2] = 0.; - //for ( std::int8_t iHit = 0; iHit < nHits; ++iHit ) { - for ( const auto& hit: track.hits()) { - const float dz = 0.001 * ( hit.z - track.zMagnet() ); - const float w = hit.weight(); - const float t = hit.sin(); + // for ( std::int8_t iHit = 0; iHit < nHits; ++iHit ) { + for ( const auto& hit : track.hits() ) { + const float dz = 0.001 * ( hit.z - track.zMagnet() ); + const float w = hit.weight(); + const float t = hit.sin(); mat[0] += w; mat[1] += w * dz; @@ -613,18 +610,17 @@ private: mat[4] += w * dz * t; mat[5] += w * t * t; - if constexpr(useProjections){ + if constexpr ( useProjections ) { const float dist = hit.projection; rhs[0] += w * dist; rhs[1] += w * dist * dz; rhs[2] += w * dist * t; - }else{ + } else { const float dist = track.distance( hit ); rhs[0] += w * dist; rhs[1] += w * dist * dz; rhs[2] += w * dist * t; } - } // -- solve the equation and update the parameters of the track @@ -645,19 +641,17 @@ private: float chi2 = track.initialChi2(); - for ( auto& hit: track.hits()) { + for ( auto& hit : track.hits() ) { const float dist = track.distance( hit ); - hit.projection = dist;// this is needed later on for fitting with projections + hit.projection = dist; // this is needed later on for fitting with projections chi2 += dist * dist * hit.weight(); // -- only flag this hit as removable if it is not alone in a plane or there are 4 planes that fired } - if( track.hits().size() > 2 ) chi2 /= ( track.hits().size() - 2 ); + if ( track.hits().size() > 2 ) chi2 /= ( track.hits().size() - 2 ); track.setChi2( chi2 ); - } - //========================================================================= // Collect the hits in the other x layer //========================================================================= @@ -735,35 +729,29 @@ private: // -- first select all hits, and then // -- accept until over a tolerance - const float xStart = track.xAtZ(preSelHits[0].z); + const float xStart = track.xAtZ( preSelHits[0].z ); const float zStart = preSelHits[0].z; for ( const auto& hit : preSelHits ) { - //const float dist = track.distance( hit ); + // const float dist = track.distance( hit ); const float dist = track.distanceFast( hit, xStart, zStart ); - if ( std::abs(dist) > tol ) continue; - //hit.projection = dist; + if ( std::abs( dist ) > tol ) continue; + // hit.projection = dist; uHitsTemp.push_back( hit ); uHitsTemp.back().projection = dist; } std::sort( uHitsTemp.begin(), uHitsTemp.end(), Downstream::IncreaseByAbsProj ); - for(auto& hit : track.hits()){ - hit.projection = track.distance( hit ); - } + for ( auto& hit : track.hits() ) { hit.projection = track.distance( hit ); } for ( const auto& hit : uHitsTemp ) { auto& greatTrack = goodXUTracks.emplace_back( track ); - - greatTrack.hits().push_back( hit ); - //fitAndRemove( greatTrack ); + // fitAndRemove( greatTrack ); simplyFit( greatTrack ); - - // -- it's sorted if ( greatTrack.chi2() > minChi2 ) { goodXUTracks.pop_back(); @@ -775,7 +763,6 @@ private: continue; } - if ( goodXUTracks.size() >= m_maxXUTracks ) { break; } } @@ -796,18 +783,18 @@ private: // It's enough to look at the closest hit only, since absolute majority of the tracks // with not a closest hit get killed later anyways // This also doesn't affect efficiency, as the addOverlapRegions function probes all close enough hits as well - float bestDist = 1e9f; + float bestDist = 1e9f; const Downstream::Hit* bestHit = nullptr; - const float xStart = track.xAtZ(preSelHits[0].z); + const float xStart = track.xAtZ( preSelHits[0].z ); const float zStart = preSelHits[0].z; - + for ( auto& hit : preSelHits ) { - //const float dist = track.distance( hit ); + // const float dist = track.distance( hit ); const float dist = track.distanceFast( hit, xStart, zStart ); - if ( std::abs(dist) > tol ) continue; - //hit.projection = dist; - if ( std::abs(dist) < std::abs(bestDist) ) { + if ( std::abs( dist ) > tol ) continue; + // hit.projection = dist; + if ( std::abs( dist ) < std::abs( bestDist ) ) { bestDist = dist; bestHit = &hit; } @@ -818,7 +805,7 @@ private: track.hits().push_back( *bestHit ); track.hits().back().projection = bestDist; simplyFit( track ); - //fitAndRemove( track ); + // fitAndRemove( track ); } track.sortFinalHits(); } @@ -844,13 +831,9 @@ private: continue; } - - tr.hits().push_back( firstHit ); tr.hits().push_back( hit ); - - if ( goodXTracks.size() >= 3 ) { break; } } @@ -913,7 +896,7 @@ private: float yTrack = track.yAtZ( hit.z ); if ( !hit.isYCompatible( yTrack, m_yTol ) ) continue; - + track.hits().push_back( hit ); hitAdded = true; } -- GitLab From e322c645e2d04327d4f36d29d99758f846018524 Mon Sep 17 00:00:00 2001 From: decianm Date: Sat, 17 May 2025 20:53:24 +0200 Subject: [PATCH 05/13] a few more optimization attempts --- Pr/PrAlgorithms/src/PrDownTrack.h | 10 ++++---- Pr/PrAlgorithms/src/PrLongLivedTracking.cpp | 27 ++++++++++++--------- 2 files changed, 21 insertions(+), 16 deletions(-) diff --git a/Pr/PrAlgorithms/src/PrDownTrack.h b/Pr/PrAlgorithms/src/PrDownTrack.h index 55fb8524d23..3b5a8d8bbe0 100644 --- a/Pr/PrAlgorithms/src/PrDownTrack.h +++ b/Pr/PrAlgorithms/src/PrDownTrack.h @@ -201,17 +201,17 @@ public: void setChi2( float chi2 ) noexcept { m_chi2 = chi2; } // functions - inline float xAtZ( float z ) const noexcept { + [[gnu::always_inline]] inline float xAtZ( float z ) const noexcept { const float curvature = 1.6e-5 * ( stateTx() - m_slopeX ); return xMagnet() + ( z - zMagnet() ) * m_slopeX + curvature * ( z - m_zUT ) * ( z - m_zUT ); } // -- this ignores the change in curvature, to be used for small distances - inline float xAtZFast( float z, float xStart, float zStart ) const noexcept { + [[gnu::always_inline]] inline float xAtZFast( float z, float xStart, float zStart ) const noexcept { return xStart + ( z - zStart ) * m_slopeX; } - inline float yAtZ( float z ) const noexcept { return yMagnet() + m_displY + ( z - zMagnet() ) * slopeY(); } + [[gnu::always_inline]] inline float yAtZ( float z ) const noexcept { return yMagnet() + m_displY + ( z - zMagnet() ) * slopeY(); } void updateX( float dx, float dsl ) noexcept { m_displX += dx; @@ -233,11 +233,11 @@ public: return sinTrack * std::abs( momentum() ); } - inline float distance( const Downstream::Hit& hit ) const noexcept { + [[gnu::always_inline]] inline float distance( const Downstream::Hit& hit ) const noexcept { return hit.xAt( yAtZ( hit.z ) ) - xAtZ( hit.z ); } - inline float distanceFast( const Downstream::Hit& hit, float zStart, float xStart ) const noexcept { + [[gnu::always_inline]] inline float distanceFast( const Downstream::Hit& hit, float zStart, float xStart ) const noexcept { return hit.xAt( yAtZ( hit.z ) ) - xAtZFast( hit.z, zStart, xStart ); } diff --git a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp index 3e256ce6a98..19cb4a0c734 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -729,8 +729,9 @@ private: // -- first select all hits, and then // -- accept until over a tolerance - const float xStart = track.xAtZ( preSelHits[0].z ); const float zStart = preSelHits[0].z; + const float xStart = track.xAtZ( zStart ); + for ( const auto& hit : preSelHits ) { // const float dist = track.distance( hit ); @@ -786,9 +787,9 @@ private: float bestDist = 1e9f; const Downstream::Hit* bestHit = nullptr; - const float xStart = track.xAtZ( preSelHits[0].z ); const float zStart = preSelHits[0].z; - + const float xStart = track.xAtZ( zStart ); + for ( auto& hit : preSelHits ) { // const float dist = track.distance( hit ); const float dist = track.distanceFast( hit, xStart, zStart ); @@ -889,15 +890,19 @@ private: for ( const auto& trackHit : track.hits() ) { if ( j++ >= trackNbHits ) break; - for ( auto& hit : preSelHits[trackHit.planeCode()] ) { + const auto planeCode = trackHit.planeCode(); + if( preSelHits[planeCode].empty() ) continue; + + const float zStart = preSelHits[planeCode].front().z; + const float xStart = track.xAtZ( zStart ); + + for ( const auto& hit : preSelHits[planeCode] ) { // -- the displacement in z between overlap modules is larger than 1mm - if ( std::abs( hit.z - trackHit.z ) < 1.0 ) continue; - if ( m_overlapTol > std::abs( track.distance( hit ) ) ) { - float yTrack = track.yAtZ( hit.z ); - - if ( !hit.isYCompatible( yTrack, m_yTol ) ) continue; - - track.hits().push_back( hit ); + // -- given that the overlap will be in the same layer as the hit on the track, we can do some heuristics for the distance before calculating the real one + if ( std::abs( hit.z - trackHit.z ) < 1.0 || std::abs( hit.x - trackHit.x) > 3*m_overlapTol) continue; + if ( m_overlapTol > std::abs( track.distanceFast( hit, xStart, zStart ) ) && hit.isYCompatible( track.yAtZ( hit.z ), m_yTol ) ) { + // -- preSelHits is not used anymore + track.hits().push_back( std::move(hit) ); hitAdded = true; } } -- GitLab From 958461b6d087072e3baa06092b459f43fc76989a Mon Sep 17 00:00:00 2001 From: decianm Date: Mon, 19 May 2025 13:15:08 +0200 Subject: [PATCH 06/13] replace another double -> float instance --- Pr/PrAlgorithms/src/PrDownTrack.h | 10 ++-------- Pr/PrAlgorithms/src/PrLongLivedTracking.cpp | 2 +- 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/Pr/PrAlgorithms/src/PrDownTrack.h b/Pr/PrAlgorithms/src/PrDownTrack.h index 3b5a8d8bbe0..55df01efd83 100644 --- a/Pr/PrAlgorithms/src/PrDownTrack.h +++ b/Pr/PrAlgorithms/src/PrDownTrack.h @@ -119,15 +119,9 @@ public: // using Hits = boost::container::static_vector; // Until we can put a bound on the number of hits, use a small_vector - PrDownTrack( Gaudi::TrackVector stateVector, float stateZ, float zUT, LHCb::span magnetParams, + PrDownTrack( Gaudi::TrackVectorF stateVector, float stateZ, float zUT, LHCb::span magnetParams, LHCb::span yParams, LHCb::span momPar, float magnetScale ) - : m_stateZ( stateZ ), m_zUT( zUT ) { - - m_stateVector[0] = (float)stateVector[0]; - m_stateVector[1] = (float)stateVector[1]; - m_stateVector[2] = (float)stateVector[2]; - m_stateVector[3] = (float)stateVector[3]; - m_stateVector[4] = (float)stateVector[4]; + : m_stateVector( stateVector ), m_stateZ( stateZ ), m_zUT( zUT ) { const auto tx2 = stateTx() * stateTx(); const auto ty2 = stateTy() * stateTy(); diff --git a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp index 19cb4a0c734..1eeb41cd8c4 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -140,7 +140,7 @@ public: // -- Note: You want the state the furthest away from the magnet, as then // the straight-line approximation is the best const auto state = _tr.StatePosDir( LHCb::Event::Enum::State::Location::EndT ); - Gaudi::TrackVector stateVector{ state.x().cast(), state.y().cast(), state.tx().cast(), state.ty().cast(), + Gaudi::TrackVectorF stateVector{ state.x().cast(), state.y().cast(), state.tx().cast(), state.ty().cast(), _tr.qOverP().cast() }; PrDownTrack refTrack( stateVector, state.z().cast(), m_zUT, m_zMagnetParams.value(), m_yParams.value(), m_momentumParams.value(), magScaleFactor * ( -1 ) ); -- GitLab From 7a993a3cdebc36edfd2cbfe4b481e023e1207ca7 Mon Sep 17 00:00:00 2001 From: decianm Date: Tue, 20 May 2025 19:36:40 +0200 Subject: [PATCH 07/13] revert a beampipe cut --- Pr/PrAlgorithms/src/PrLongLivedTracking.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp index 1eeb41cd8c4..f8a1c1cdada 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -381,7 +381,7 @@ private: if ( xuTrack.hits().size() < 3 ) continue; simplyFit( xuTrack ); // fitAndRemove( xuTrack ); - if ( xuTrack.chi2() > m_maxChi2 || !xuTrack.isYCompatible( m_yTol ) || insideBeampipe( xuTrack ) ) continue; + if ( xuTrack.chi2() > m_maxChi2 || !xuTrack.isYCompatible( m_yTol ) ) continue; trackCandidates.push_back( std::move( xuTrack ) ); } // Loop over good xu tracks } // Loop over good x tracks -- GitLab From ce4fd3da8556948e3d765ac1eb2cdf1fc8d8920b Mon Sep 17 00:00:00 2001 From: Gitlab CI Date: Wed, 21 May 2025 06:34:54 +0000 Subject: [PATCH 08/13] pre-commit fixes patch generated by https://gitlab.cern.ch/lhcb/Rec/-/jobs/56191350 --- Pr/PrAlgorithms/src/PrDownTrack.h | 9 ++++--- Pr/PrAlgorithms/src/PrLongLivedTracking.cpp | 29 +++++++++++---------- 2 files changed, 21 insertions(+), 17 deletions(-) diff --git a/Pr/PrAlgorithms/src/PrDownTrack.h b/Pr/PrAlgorithms/src/PrDownTrack.h index 55df01efd83..ddce93c8291 100644 --- a/Pr/PrAlgorithms/src/PrDownTrack.h +++ b/Pr/PrAlgorithms/src/PrDownTrack.h @@ -121,7 +121,7 @@ public: PrDownTrack( Gaudi::TrackVectorF stateVector, float stateZ, float zUT, LHCb::span magnetParams, LHCb::span yParams, LHCb::span momPar, float magnetScale ) - : m_stateVector( stateVector ), m_stateZ( stateZ ), m_zUT( zUT ) { + : m_stateVector( stateVector ), m_stateZ( stateZ ), m_zUT( zUT ) { const auto tx2 = stateTx() * stateTx(); const auto ty2 = stateTy() * stateTy(); @@ -205,7 +205,9 @@ public: return xStart + ( z - zStart ) * m_slopeX; } - [[gnu::always_inline]] inline float yAtZ( float z ) const noexcept { return yMagnet() + m_displY + ( z - zMagnet() ) * slopeY(); } + [[gnu::always_inline]] inline float yAtZ( float z ) const noexcept { + return yMagnet() + m_displY + ( z - zMagnet() ) * slopeY(); + } void updateX( float dx, float dsl ) noexcept { m_displX += dx; @@ -231,7 +233,8 @@ public: return hit.xAt( yAtZ( hit.z ) ) - xAtZ( hit.z ); } - [[gnu::always_inline]] inline float distanceFast( const Downstream::Hit& hit, float zStart, float xStart ) const noexcept { + [[gnu::always_inline]] inline float distanceFast( const Downstream::Hit& hit, float zStart, + float xStart ) const noexcept { return hit.xAt( yAtZ( hit.z ) ) - xAtZFast( hit.z, zStart, xStart ); } diff --git a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp index f8a1c1cdada..2e38091dfd5 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -139,11 +139,11 @@ public: // -- Note: You want the state the furthest away from the magnet, as then // the straight-line approximation is the best - const auto state = _tr.StatePosDir( LHCb::Event::Enum::State::Location::EndT ); + const auto state = _tr.StatePosDir( LHCb::Event::Enum::State::Location::EndT ); Gaudi::TrackVectorF stateVector{ state.x().cast(), state.y().cast(), state.tx().cast(), state.ty().cast(), - _tr.qOverP().cast() }; - PrDownTrack refTrack( stateVector, state.z().cast(), m_zUT, m_zMagnetParams.value(), m_yParams.value(), - m_momentumParams.value(), magScaleFactor * ( -1 ) ); + _tr.qOverP().cast() }; + PrDownTrack refTrack( stateVector, state.z().cast(), m_zUT, m_zMagnetParams.value(), m_yParams.value(), + m_momentumParams.value(), magScaleFactor * ( -1 ) ); if ( ( std::abs( refTrack.momentum() ) < m_initialMinMomentum ) || ( refTrack.pt() < m_initialMinPt ) ) continue; @@ -731,7 +731,6 @@ private: // -- accept until over a tolerance const float zStart = preSelHits[0].z; const float xStart = track.xAtZ( zStart ); - for ( const auto& hit : preSelHits ) { // const float dist = track.distance( hit ); @@ -789,7 +788,7 @@ private: const float zStart = preSelHits[0].z; const float xStart = track.xAtZ( zStart ); - + for ( auto& hit : preSelHits ) { // const float dist = track.distance( hit ); const float dist = track.distanceFast( hit, xStart, zStart ); @@ -891,18 +890,20 @@ private: if ( j++ >= trackNbHits ) break; const auto planeCode = trackHit.planeCode(); - if( preSelHits[planeCode].empty() ) continue; - + if ( preSelHits[planeCode].empty() ) continue; + const float zStart = preSelHits[planeCode].front().z; const float xStart = track.xAtZ( zStart ); - + for ( const auto& hit : preSelHits[planeCode] ) { // -- the displacement in z between overlap modules is larger than 1mm - // -- given that the overlap will be in the same layer as the hit on the track, we can do some heuristics for the distance before calculating the real one - if ( std::abs( hit.z - trackHit.z ) < 1.0 || std::abs( hit.x - trackHit.x) > 3*m_overlapTol) continue; - if ( m_overlapTol > std::abs( track.distanceFast( hit, xStart, zStart ) ) && hit.isYCompatible( track.yAtZ( hit.z ), m_yTol ) ) { - // -- preSelHits is not used anymore - track.hits().push_back( std::move(hit) ); + // -- given that the overlap will be in the same layer as the hit on the track, we can do some heuristics for + // the distance before calculating the real one + if ( std::abs( hit.z - trackHit.z ) < 1.0 || std::abs( hit.x - trackHit.x ) > 3 * m_overlapTol ) continue; + if ( m_overlapTol > std::abs( track.distanceFast( hit, xStart, zStart ) ) && + hit.isYCompatible( track.yAtZ( hit.z ), m_yTol ) ) { + // -- preSelHits is not used anymore + track.hits().push_back( std::move( hit ) ); hitAdded = true; } } -- GitLab From a109b22cad336bac558cff06ea219db2388fd57a Mon Sep 17 00:00:00 2001 From: decianm Date: Wed, 21 May 2025 15:55:11 +0200 Subject: [PATCH 09/13] adapt default number of UT hits and naming of approximate extrapolation --- Pr/PrAlgorithms/src/PrDownTrack.h | 8 ++++---- Pr/PrAlgorithms/src/PrLongLivedTracking.cpp | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/Pr/PrAlgorithms/src/PrDownTrack.h b/Pr/PrAlgorithms/src/PrDownTrack.h index ddce93c8291..07966d9723f 100644 --- a/Pr/PrAlgorithms/src/PrDownTrack.h +++ b/Pr/PrAlgorithms/src/PrDownTrack.h @@ -115,7 +115,7 @@ namespace Downstream { class PrDownTrack final { public: - using Hits = boost::container::small_vector>; + using Hits = boost::container::small_vector>; // using Hits = boost::container::static_vector; // Until we can put a bound on the number of hits, use a small_vector @@ -201,7 +201,7 @@ public: } // -- this ignores the change in curvature, to be used for small distances - [[gnu::always_inline]] inline float xAtZFast( float z, float xStart, float zStart ) const noexcept { + [[gnu::always_inline]] inline float xAtZLinear( float z, float xStart, float zStart ) const noexcept { return xStart + ( z - zStart ) * m_slopeX; } @@ -233,9 +233,9 @@ public: return hit.xAt( yAtZ( hit.z ) ) - xAtZ( hit.z ); } - [[gnu::always_inline]] inline float distanceFast( const Downstream::Hit& hit, float zStart, + [[gnu::always_inline]] inline float distanceLinear( const Downstream::Hit& hit, float zStart, float xStart ) const noexcept { - return hit.xAt( yAtZ( hit.z ) ) - xAtZFast( hit.z, zStart, xStart ); + return hit.xAt( yAtZ( hit.z ) ) - xAtZLinear( hit.z, zStart, xStart ); } bool isYCompatible( const float tol ) const noexcept { diff --git a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp index 2e38091dfd5..2a873852010 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -734,7 +734,7 @@ private: for ( const auto& hit : preSelHits ) { // const float dist = track.distance( hit ); - const float dist = track.distanceFast( hit, xStart, zStart ); + const float dist = track.distanceLinear( hit, xStart, zStart ); if ( std::abs( dist ) > tol ) continue; // hit.projection = dist; uHitsTemp.push_back( hit ); @@ -791,7 +791,7 @@ private: for ( auto& hit : preSelHits ) { // const float dist = track.distance( hit ); - const float dist = track.distanceFast( hit, xStart, zStart ); + const float dist = track.distanceLinear( hit, xStart, zStart ); if ( std::abs( dist ) > tol ) continue; // hit.projection = dist; if ( std::abs( dist ) < std::abs( bestDist ) ) { @@ -900,7 +900,7 @@ private: // -- given that the overlap will be in the same layer as the hit on the track, we can do some heuristics for // the distance before calculating the real one if ( std::abs( hit.z - trackHit.z ) < 1.0 || std::abs( hit.x - trackHit.x ) > 3 * m_overlapTol ) continue; - if ( m_overlapTol > std::abs( track.distanceFast( hit, xStart, zStart ) ) && + if ( m_overlapTol > std::abs( track.distanceLinear( hit, xStart, zStart ) ) && hit.isYCompatible( track.yAtZ( hit.z ), m_yTol ) ) { // -- preSelHits is not used anymore track.hits().push_back( std::move( hit ) ); -- GitLab From 2987e85c822d561e89836a82b72cb8bc01a94ac2 Mon Sep 17 00:00:00 2001 From: Gitlab CI Date: Wed, 21 May 2025 13:56:01 +0000 Subject: [PATCH 10/13] pre-commit fixes patch generated by https://gitlab.cern.ch/lhcb/Rec/-/jobs/56333737 --- Pr/PrAlgorithms/src/PrDownTrack.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Pr/PrAlgorithms/src/PrDownTrack.h b/Pr/PrAlgorithms/src/PrDownTrack.h index 07966d9723f..4d14f96e9b4 100644 --- a/Pr/PrAlgorithms/src/PrDownTrack.h +++ b/Pr/PrAlgorithms/src/PrDownTrack.h @@ -234,7 +234,7 @@ public: } [[gnu::always_inline]] inline float distanceLinear( const Downstream::Hit& hit, float zStart, - float xStart ) const noexcept { + float xStart ) const noexcept { return hit.xAt( yAtZ( hit.z ) ) - xAtZLinear( hit.z, zStart, xStart ); } -- GitLab From f91b1501accbb994f635cc8516fe1f2a78b1ea7f Mon Sep 17 00:00:00 2001 From: decianm Date: Wed, 21 May 2025 16:40:10 +0200 Subject: [PATCH 11/13] fix a few small cosmetic things --- Pr/PrAlgorithms/src/PrLongLivedTracking.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp index 2a873852010..dc5c535c0f4 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -463,8 +463,7 @@ private: auto xx = mH.get().cast() + y * dxDy; if ( xPredTol < std::abs( pos - xx ) ) continue; // go from -x to +x - if ( xPredTol < xx - pos ) break; // can break if we go out of the right bound - + preSelHits[layer].emplace_back( &hitHandler, itHit, xx, zLayer, fabs( xx - pos ) ); } } @@ -579,8 +578,7 @@ private: } while ( again ); } //========================================================================= - // Fit and remove the worst hit, as long as over tolerance - // Perform a chi2 fit to the track and remove outliers + // Simplified fit function that only fits and does not perform outlier removal //========================================================================= template void simplyFit( PrDownTrack& track ) const { -- GitLab From 7f254c843e5aa8b0bda5243ee6ceb465f7b8cf3a Mon Sep 17 00:00:00 2001 From: decianm Date: Wed, 21 May 2025 21:07:57 +0200 Subject: [PATCH 12/13] move momentum and pT to part where the direction of the UT-track segment is better known, apply sanity pT cut of 30MeV --- Pr/PrAlgorithms/src/PrLongLivedTracking.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp index dc5c535c0f4..330a6cf40d6 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -145,8 +145,6 @@ public: PrDownTrack refTrack( stateVector, state.z().cast(), m_zUT, m_zMagnetParams.value(), m_yParams.value(), m_momentumParams.value(), magScaleFactor * ( -1 ) ); - if ( ( std::abs( refTrack.momentum() ) < m_initialMinMomentum ) || ( refTrack.pt() < m_initialMinPt ) ) continue; - // -- Veto particles coming from the beam pipe. if ( insideBeampipe( refTrack ) ) continue; @@ -280,7 +278,7 @@ private: // - InitialMinMomentum: Minimum momentum of the track from initial estimate Gaudi::Property m_initialMinMomentum{ this, "InitialMinMomentum", 1400. * Gaudi::Units::MeV }; // - MinPt: Minimum pT of the track - Gaudi::Property m_minPt{ this, "MinPt", 0. * Gaudi::Units::MeV }; + Gaudi::Property m_minPt{ this, "MinPt", 30. * Gaudi::Units::MeV }; // - MinMomentum: Minimum momentum of the track Gaudi::Property m_minMomentum{ this, "MinMomentum", 0. * Gaudi::Units::MeV }; @@ -362,6 +360,8 @@ private: const float slopeX = ( track.xMagnet() - posX ) / ( track.zMagnet() - meanZ ); track.setSlopeX( slopeX ); + if ( ( std::abs( track.momentum() ) < m_initialMinMomentum ) || ( track.pt() < m_initialMinPt ) ) continue; + // -- Fit x projection findMatchingHits( track, preSelHits[plane_indx[1]], matchingXHits ); fitXProjection( track, myHit, matchingXHits, goodXTracks, skip_second_layer ); -- GitLab From 8678b71559d3a0e5ffd5f31f34ecb99990fbe3df Mon Sep 17 00:00:00 2001 From: Gitlab CI Date: Wed, 21 May 2025 19:08:45 +0000 Subject: [PATCH 13/13] pre-commit fixes patch generated by https://gitlab.cern.ch/lhcb/Rec/-/jobs/56348885 --- Pr/PrAlgorithms/src/PrLongLivedTracking.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp index 330a6cf40d6..3f1002e50c4 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -361,7 +361,7 @@ private: track.setSlopeX( slopeX ); if ( ( std::abs( track.momentum() ) < m_initialMinMomentum ) || ( track.pt() < m_initialMinPt ) ) continue; - + // -- Fit x projection findMatchingHits( track, preSelHits[plane_indx[1]], matchingXHits ); fitXProjection( track, myHit, matchingXHits, goodXTracks, skip_second_layer ); @@ -463,7 +463,7 @@ private: auto xx = mH.get().cast() + y * dxDy; if ( xPredTol < std::abs( pos - xx ) ) continue; // go from -x to +x - + preSelHits[layer].emplace_back( &hitHandler, itHit, xx, zLayer, fabs( xx - pos ) ); } } -- GitLab