diff --git a/Pr/PrAlgorithms/src/PrDownTrack.h b/Pr/PrAlgorithms/src/PrDownTrack.h index 67a5b180bac4441d77abfb6dfdfd8a01d6b7108a..4d14f96e9b4ddbbbb73c959216f4a3ce3047fb70 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(); } @@ -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 @@ -107,13 +115,14 @@ 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 - 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_stateVector( stateVector ), m_stateZ( stateZ ), m_zUT( zUT ) { + const auto tx2 = stateTx() * stateTx(); const auto ty2 = stateTy() * stateTy(); m_momentumParam = ( momPar[0] + momPar[1] * tx2 + momPar[2] * ty2 ) * magnetScale; @@ -186,12 +195,19 @@ public: void setChi2( float chi2 ) noexcept { m_chi2 = chi2; } // functions - 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 ); } - float yAtZ( float z ) const noexcept { return yMagnet() + m_displY + ( z - zMagnet() ) * slopeY(); } + // -- this ignores the change in curvature, to be used for small distances + [[gnu::always_inline]] inline float xAtZLinear( float z, float xStart, float zStart ) const noexcept { + return xStart + ( z - zStart ) * m_slopeX; + } + + [[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; @@ -213,7 +229,14 @@ public: return sinTrack * std::abs( momentum() ); } - float distance( const Downstream::Hit& hit ) const noexcept { return hit.xAt( yAtZ( hit.z ) ) - xAtZ( hit.z ); } + [[gnu::always_inline]] inline float distance( const Downstream::Hit& hit ) const noexcept { + return hit.xAt( yAtZ( hit.z ) ) - xAtZ( hit.z ); + } + + [[gnu::always_inline]] inline float distanceLinear( const Downstream::Hit& hit, float zStart, + float xStart ) const noexcept { + return hit.xAt( yAtZ( hit.z ) ) - xAtZLinear( hit.z, zStart, xStart ); + } bool isYCompatible( const float tol ) const noexcept { return std ::all_of( m_hits.begin(), m_hits.end(), @@ -227,9 +250,10 @@ public: } private: - Gaudi::TrackVector m_stateVector; - float m_stateZ{ 0 }; - Gaudi::XYZPoint m_magnet; + // -- using the float versions everywhere + Gaudi::TrackVectorF m_stateVector; + float m_stateZ{ 0 }; + 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 3361d140b2d198c54c25855871ccf1966e747136..3f1002e50c49c1bd8d7848f817b08c0bd2f3996b 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -139,13 +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 ); - Gaudi::TrackVector 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 ) ); - - if ( ( std::abs( refTrack.momentum() ) < m_initialMinMomentum ) || ( refTrack.pt() < m_initialMinPt ) ) continue; + 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 ) ); // -- 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 ); @@ -379,7 +379,8 @@ private: 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 +462,7 @@ 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 preSelHits[layer].emplace_back( &hitHandler, itHit, xx, zLayer, fabs( xx - pos ) ); } @@ -478,7 +478,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 +498,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 +520,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 +549,18 @@ 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,14 +571,84 @@ 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; } } while ( again ); } + //========================================================================= + // Simplified fit function that only fits and does not perform outlier removal + //========================================================================= + 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 ); + + float chi2 = track.initialChi2(); + + 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 + } + + if ( track.hits().size() > 2 ) chi2 /= ( track.hits().size() - 2 ); + track.setChi2( chi2 ); + } //========================================================================= // Collect the hits in the other x layer @@ -645,7 +712,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,31 +727,40 @@ 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; + const float zStart = preSelHits[0].z; + const float xStart = track.xAtZ( zStart ); + + for ( const auto& hit : preSelHits ) { + // const float dist = track.distance( hit ); + const float dist = track.distanceLinear( hit, xStart, zStart ); + if ( std::abs( dist ) > tol ) continue; + // hit.projection = dist; uHitsTemp.push_back( hit ); + uHitsTemp.back().projection = dist; } - std::sort( uHitsTemp.begin(), uHitsTemp.end(), Downstream::IncreaseByProj ); + std::sort( uHitsTemp.begin(), uHitsTemp.end(), Downstream::IncreaseByAbsProj ); - for ( auto& hit : uHitsTemp ) { + 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 ); - - if ( !greatTrack.isYCompatible( m_yTol ) ) { - goodXUTracks.pop_back(); - continue; - } + // fitAndRemove( greatTrack ); + simplyFit( greatTrack ); // -- 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; } } @@ -695,23 +771,28 @@ 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; + float bestDist = 1e9f; + const Downstream::Hit* bestHit = nullptr; + + const float zStart = preSelHits[0].z; + const float xStart = track.xAtZ( zStart ); + 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 ); + const float dist = track.distanceLinear( hit, xStart, zStart ); + if ( std::abs( dist ) > tol ) continue; + // hit.projection = dist; + if ( std::abs( dist ) < std::abs( bestDist ) ) { bestDist = dist; bestHit = &hit; } @@ -720,7 +801,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(); } @@ -730,7 +813,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(); @@ -745,6 +828,7 @@ private: goodXTracks.pop_back(); continue; } + tr.hits().push_back( firstHit ); tr.hits().push_back( hit ); @@ -794,7 +878,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(); @@ -803,24 +887,29 @@ private: for ( const auto& trackHit : track.hits() ) { if ( j++ >= trackNbHits ) break; - for ( auto& hit : preSelHits[trackHit.planeCode()] ) { - 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; - } + 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 + // -- 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.distanceLinear( 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; } } } if ( hitAdded ) { track.sortFinalHits(); - fitAndRemove( track ); + fitAndRemove( track ); } }