diff --git a/Pr/PrAlgorithms/src/PrDownTrack.h b/Pr/PrAlgorithms/src/PrDownTrack.h index 4d14f96e9b4ddbbbb73c959216f4a3ce3047fb70..8cb7dcbe506d8d642b109bf86b6c75cc09d5d7a8 100644 --- a/Pr/PrAlgorithms/src/PrDownTrack.h +++ b/Pr/PrAlgorithms/src/PrDownTrack.h @@ -32,41 +32,42 @@ namespace Downstream { public: const LHCb::Pr::UT::Hits* hits = nullptr; int hit{ 0 }; - float x{ 0 }, z{ 0 }; + float x{ 0 }, x0{ 0 }, dxdy{ 0 }, z{ 0 }; float projection{ 0 }; using F = SIMDWrapper::scalar::types::float_v; using I = SIMDWrapper::scalar::types::int_v; - Hit( const LHCb::Pr::UT::Hits* _hits, const int _hit, float _x, float _z, float _proj ) - : hits( _hits ), hit( _hit ), x( _x ), z( _z ), projection( _proj ) {} + Hit( const LHCb::Pr::UT::Hits* _hits, const int _hit, float _x, float _x0, float _dxdy, float _z, float _proj ) + : hits( _hits ), hit( _hit ), x( _x ), x0( _x0 ), dxdy( _dxdy ), z( _z ), projection( _proj ) {} - [[nodiscard]] auto lhcbID() const { + [[nodiscard, gnu::always_inline]] inline auto lhcbID() const { const auto mH = getScalarHit(); const auto chanID = mH.get().cast(); return bit_cast( LHCb::LHCbID( LHCb::Detector::UT::ChannelID( chanID ) ).lhcbID() ); } - [[nodiscard]] int planeCode() const { + [[nodiscard, gnu::always_inline]] inline int planeCode() const { const auto mH = getScalarHit(); auto lhcbid = mH.get().cast(); return ( lhcbid & static_cast( UTInfo::MasksBits::HalfLayerMask ) ) >> static_cast( UTInfo::MasksBits::HalfLayerBits ); } - [[nodiscard]] auto weight() const { + [[nodiscard, gnu::always_inline]] inline auto weight() const { const auto mH = getScalarHit(); return mH.get().cast(); } - [[nodiscard]] auto sin() const { + [[nodiscard, gnu::always_inline]] inline auto sin() const { const auto mH = getScalarHit(); - return -mH.get().cast() * mH.get().cast(); + return -dxdy * mH.get().cast(); } - [[nodiscard]] auto zAtYEq0() const { + [[nodiscard, gnu::always_inline]] inline auto zAtYEq0() const { const auto mH = getScalarHit(); return mH.get().cast(); } - [[nodiscard]] bool isYCompatible( const float y, const float tol ) const { + [[nodiscard, gnu::always_inline]] inline auto z0() const { return z; } + [[nodiscard, gnu::always_inline]] inline bool isYCompatible( const float y, const float tol ) const { const auto mH = getScalarHit(); const auto yMin = std::min( mH.get().cast(), mH.get().cast() ); @@ -74,9 +75,20 @@ namespace Downstream { std::max( mH.get().cast(), mH.get().cast() ); return ( ( ( yMin - tol ) <= y ) && ( y <= ( yMax + tol ) ) ); } - [[nodiscard]] inline auto xAt( const float y ) const { + template + [[nodiscard, gnu::always_inline]] inline float_t xAtY( const float_t y ) const { + return x0 + y * dxdy; + } + [[nodiscard, gnu::always_inline]] inline auto dxDy() const { return dxdy; } + [[nodiscard, gnu::always_inline]] inline auto ymin() const { + const auto mH = getScalarHit(); + return std::min( mH.get().cast(), + mH.get().cast() ); + } + [[nodiscard, gnu::always_inline]] inline auto ymax() const { const auto mH = getScalarHit(); - return mH.get().cast() + y * mH.get().cast(); + return std::max( mH.get().cast(), + mH.get().cast() ); } }; @@ -96,6 +108,81 @@ namespace Downstream { lhs.lhcbID() < rhs.lhcbID() ); }; + /** + * @namespace Downstream::SOA + * Provides structure-of-arrays (SOA) representations and proxies, optimized for SIMD processing. + * + * @struct Downstream::SOA::Hits + * SOA collection of UT hits. Only selected variables are buffered. + * Container has an idx field which may store the index of the same hit in hitHandler. + * + * @author Volodymyr Svintozelskyi + * @date 2025-06-08 + */ + namespace SOA { + struct HitsTag { + // Extrapolation + struct x0 : LHCb::Event::float_field {}; + struct z0 : LHCb::Event::float_field {}; + struct dxdy : LHCb::Event::float_field {}; + struct weight : LHCb::Event::float_field {}; + struct ymin : LHCb::Event::float_field {}; + struct ymax : LHCb::Event::float_field {}; + struct idx : LHCb::Event::int_field {}; + + template + using hits_t = LHCb::Event::SOACollection; + }; + + struct Hits : HitsTag::hits_t { + using base_t = typename HitsTag::hits_t; + using tag_t = HitsTag; + using base_t::allocator_type; + Hits( Zipping::ZipFamilyNumber zipIdentifier = Zipping::generateZipIdentifier(), allocator_type alloc = {} ) + : base_t{ std::move( zipIdentifier ), std::move( alloc ) } {} + + // Constructor used by zipping machinery when making a copy of a zip + Hits( Zipping::ZipFamilyNumber zn, Hits const& old ) : base_t{ std::move( zn ), old } {} + + // Define an minimal custom proxy for this track + template + struct HitsProxy : LHCb::Event::Proxy { + using base_t = typename LHCb::Event::Proxy; + using base_t::base_t; + using base_t::width; + + [[nodiscard, gnu::always_inline]] inline auto x0() const { return this->template get(); } + [[nodiscard, gnu::always_inline]] inline auto z0() const { return this->template get(); } + [[nodiscard, gnu::always_inline]] inline auto dxdy() const { return this->template get(); } + [[nodiscard, gnu::always_inline]] inline auto weight() const { return this->template get(); } + [[nodiscard, gnu::always_inline]] inline auto sin() const { return -dxdy() / sqrt( 1.f + dxdy() * dxdy() ); } + [[nodiscard, gnu::always_inline]] inline auto idx() const { return this->template get(); } + [[nodiscard, gnu::always_inline]] inline auto ymin() const { return this->template get(); } + [[nodiscard, gnu::always_inline]] inline auto ymax() const { return this->template get(); } + + template + [[nodiscard, gnu::always_inline]] inline auto xAtY( const float_t y ) const { + return x0() + y * dxdy(); + } + + [[gnu::always_inline]] inline void copyFromScalarHit( const Downstream::Hit& hit, unsigned idx = 0 ) const + requires( simd == SIMDWrapper::InstructionSet::Scalar ) + { + this->template field().set( hit.x0 ); + this->template field().set( hit.z ); + this->template field().set( hit.dxDy() ); + this->template field().set( hit.weight() ); + this->template field().set( idx ); + this->template field().set( hit.ymin() ); + this->template field().set( hit.ymax() ); + } + }; + + template + using proxy_type = HitsProxy; + }; + } // namespace SOA + } // namespace Downstream /** @class PrDownTrack PrDownTrack.h @@ -119,6 +206,8 @@ public: // using Hits = boost::container::static_vector; // Until we can put a bound on the number of hits, use a small_vector + PrDownTrack() = default; + 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 ) { @@ -172,30 +261,33 @@ public: } /// getters - float stateX() const noexcept { return m_stateVector[0]; } - float stateY() const noexcept { return m_stateVector[1]; } - float stateZ() const noexcept { return m_stateZ; } - float stateTx() const noexcept { return m_stateVector[2]; } - float stateTy() const noexcept { return m_stateVector[3]; } - float stateQoP() const noexcept { return m_stateVector[4]; } - Hits& hits() noexcept { return m_hits; } - const Hits& hits() const noexcept { return m_hits; } - float xMagnet() const noexcept { return m_magnet.x(); } - float yMagnet() const noexcept { return m_magnet.y(); } - float zMagnet() const noexcept { return m_magnet.z(); } - float slopeX() const noexcept { return m_slopeX; } - float slopeY() const noexcept { return m_slopeY; } - float weightXMag() const noexcept { return m_weightXMag; } - float weightYMag() const noexcept { return m_weightYMag; } - float chi2() const noexcept { return m_chi2; } - float zUT() const noexcept { return m_zUT; } + [[gnu::always_inline]] inline float stateX() const noexcept { return m_stateVector[0]; } + [[gnu::always_inline]] inline float stateY() const noexcept { return m_stateVector[1]; } + [[gnu::always_inline]] inline float stateZ() const noexcept { return m_stateZ; } + [[gnu::always_inline]] inline float stateTx() const noexcept { return m_stateVector[2]; } + [[gnu::always_inline]] inline float stateTy() const noexcept { return m_stateVector[3]; } + [[gnu::always_inline]] inline float stateQoP() const noexcept { return m_stateVector[4]; } + [[gnu::always_inline]] inline Hits& hits() noexcept { return m_hits; } + [[gnu::always_inline]] inline const Hits& hits() const noexcept { return m_hits; } + [[gnu::always_inline]] inline float xMagnet() const noexcept { return m_magnet.x(); } + [[gnu::always_inline]] inline float yMagnet() const noexcept { return m_magnet.y(); } + [[gnu::always_inline]] inline float zMagnet() const noexcept { return m_magnet.z(); } + [[gnu::always_inline]] inline float slopeX() const noexcept { return m_slopeX; } + [[gnu::always_inline]] inline float slopeY() const noexcept { return m_slopeY; } + [[gnu::always_inline]] inline float weightXMag() const noexcept { return m_weightXMag; } + [[gnu::always_inline]] inline float weightYMag() const noexcept { return m_weightYMag; } + [[gnu::always_inline]] inline float chi2() const noexcept { return m_chi2; } + [[gnu::always_inline]] inline float zUT() const noexcept { return m_zUT; } + [[gnu::always_inline]] inline float displX() const noexcept { return m_displX; } + [[gnu::always_inline]] inline float displY() const noexcept { return m_displY; } /// setters void setSlopeX( float slopeX ) noexcept { m_slopeX = slopeX; } void setChi2( float chi2 ) noexcept { m_chi2 = chi2; } // functions - [[gnu::always_inline]] inline float xAtZ( float z ) const noexcept { + template + [[gnu::always_inline]] inline float_t xAtZ( float_t 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 ); } @@ -205,7 +297,8 @@ public: return xStart + ( z - zStart ) * m_slopeX; } - [[gnu::always_inline]] inline float yAtZ( float z ) const noexcept { + template + [[gnu::always_inline]] inline float_t yAtZ( float_t z ) const noexcept { return yMagnet() + m_displY + ( z - zMagnet() ) * slopeY(); } @@ -230,15 +323,20 @@ public: } [[gnu::always_inline]] inline float distance( const Downstream::Hit& hit ) const noexcept { - return hit.xAt( yAtZ( hit.z ) ) - xAtZ( hit.z ); + return hit.xAtY( yAtZ( hit.z ) ) - xAtZ( hit.z ); + } + + template + [[gnu::always_inline]] inline auto distance( const HitProxy& hit ) const noexcept { + return hit.x0() + hit.dxdy() * yAtZ( hit.z0() ) - xAtZ( hit.z0() ); } [[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 ); + return hit.xAtY( yAtZ( hit.z ) ) - xAtZLinear( hit.z, zStart, xStart ); } - bool isYCompatible( const float tol ) const noexcept { + [[gnu::always_inline]] inline 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 ); } ); } @@ -270,3 +368,171 @@ private: // collection of downstream tracks... From PatDownTrack using PrDownTracks = std::vector; + +template +class PrIncompleteDownTrack final { + using int_t = typename simd_t::int_v; + using float_t = typename simd_t::float_v; + using mask_t = typename simd_t::mask_v; + +public: + PrIncompleteDownTrack( const PrDownTrack& _track, const LHCb::LinAlg::Vec& _fit_res ) + : track( _track ), fit_res( _fit_res ) {} + + [[gnu::always_inline]] inline const float_t dx() const { return fit_res( 0 ); } + + [[gnu::always_inline]] inline const float_t dsl() const { return fit_res( 1 ); } + + [[gnu::always_inline]] inline const float_t dy() const + requires( dim == 3 ) + { + return fit_res( 2 ); + } + + [[gnu::always_inline]] inline const float_t mag_x() const { return track.xMagnet() + dx(); } + + [[gnu::always_inline]] inline const float_t slopeX() const { return track.slopeX() + dsl(); } + + [[gnu::always_inline]] inline const float_t displX() const { return track.displX() + dx(); } + + [[gnu::always_inline]] inline const float_t displY() const { + if constexpr ( dim == 2 ) { + return 0.f; + } else { + return track.displY() + dy(); + }; + } + + [[gnu::always_inline]] inline const float_t curvature() const { return 1.6e-5 * ( track.stateTx() - slopeX() ); } + + template + [[gnu::always_inline]] inline const float_t xAtZ( const arg_t z ) const { + return mag_x() + ( z - track.zMagnet() ) * slopeX() + curvature() * ( z - track.zUT() ) * ( z - track.zUT() ); + } + + template + [[gnu::always_inline]] inline const float_t yAtZ( const arg_t z ) const { + return track.yMagnet() + displY() + ( z - track.zMagnet() ) * track.slopeY(); + } + + template + [[gnu::always_inline]] inline const float_t chi2( const HitProxy& newhit, const float y_tol ) const + requires( dim == 3 ) + { + float_t chi2 = displX() * displX() * track.weightXMag() + displY() * displY() * track.weightYMag(); + mask_t ycompatible( true ); + + for ( const auto& h : track.hits() ) { + const auto y = yAtZ( h.z0() ); + const auto dist = h.xAtY( y ) - xAtZ( h.z0() ); + chi2 += dist * dist * h.weight(); + ycompatible = ycompatible & ( y >= ( h.ymin() - y_tol ) ) & ( y <= ( h.ymax() + y_tol ) ); + } + + const auto track_xAtZ = xAtZ( newhit.z0() ); + const auto track_yAtZ = yAtZ( newhit.z0() ); + const auto hit_xAtY = newhit.xAtY( track_yAtZ ); + chi2 += newhit.weight() * ( hit_xAtY - track_xAtZ ) * ( hit_xAtY - track_xAtZ ); + if ( track.hits().size() > 2 ) chi2 /= ( track.hits().size() - 2 ); + + ycompatible = + ycompatible & ( track_yAtZ >= ( newhit.ymin() - y_tol ) ) & ( track_yAtZ <= ( newhit.ymax() + y_tol ) ); + return select( ycompatible, chi2, std::numeric_limits::max() ); + } + + template + [[gnu::always_inline]] inline const float_t chi2( const HitProxyA& firstHit, const HitProxyB& secondHit ) const + requires( dim == 2 ) + { + float_t chi2 = displX() * displX() * track.weightXMag(); + + const auto getDev2 = [&]( const auto& hit ) { + const auto v = hit.xAtY( yAtZ( hit.z0() ) ) - xAtZ( hit.z0() ); + return v * v; + }; + + chi2 += firstHit.weight() * getDev2( firstHit ); + chi2 += secondHit.weight() * getDev2( secondHit ); + + return chi2; + } + + template + [[gnu::always_inline]] void inline compressstore( PrDownTracks& output, const mask_t hit_ok, const float_t chi2, + const int_t idx, const Downstream::Hits& preSelHits ) const + requires( dim == 3 ) + { + + if ( !any( hit_ok ) ) return; + + std::array chi2_buf, dx_buf, dsl_buf, dy_buf; + std::array idx_buf; + + // Move the results from simd vectors to the arrays, so we can iterate over them + chi2.compressstore( hit_ok, chi2_buf.data() ); + dx().compressstore( hit_ok, dx_buf.data() ); + dsl().compressstore( hit_ok, dsl_buf.data() ); + dy().compressstore( hit_ok, dy_buf.data() ); + idx.compressstore( hit_ok, idx_buf.data() ); + + // Finally store the fitted tracks in the output vector + const short ngoodhits = popcount( hit_ok ); + output.reserve( output.size() + ngoodhits ); + + for ( short i = 0; i < ngoodhits; i++ ) { + auto& greatTrack = output.emplace_back( track ); + greatTrack.hits().push_back( preSelHits[idx_buf[i]] ); + greatTrack.updateX( dx_buf[i], dsl_buf[i] ); + greatTrack.updateY( dy_buf[i] ); + greatTrack.setChi2( chi2_buf[i] ); + + if constexpr ( update_proj ) { + for ( auto& hit : greatTrack.hits() ) { + const float dist = greatTrack.distance( hit ); + hit.projection = dist; + } + } + } + } + + template + [[gnu::always_inline]] void inline compressstore( PrDownTracks& output, const mask_t hit_ok, const float_t chi2, + const int_t idx, const HitProxy& firstHit, + const Downstream::Hits& preSelHits ) const + requires( dim == 2 ) + { + if ( !any( hit_ok ) ) return; + + std::array chi2_buf, dx_buf, dsl_buf; + std::array idx_buf; + + // Move the results from simd vectors to the arrays, so we can iterate over them + chi2.compressstore( hit_ok, chi2_buf.data() ); + dx().compressstore( hit_ok, dx_buf.data() ); + dsl().compressstore( hit_ok, dsl_buf.data() ); + idx.compressstore( hit_ok, idx_buf.data() ); + + // Finally store the fitted tracks in the output vector + const short ngoodhits = popcount( hit_ok ); + output.reserve( output.size() + ngoodhits ); + + for ( short i = 0; i < ngoodhits; i++ ) { + auto& greatTrack = output.emplace_back( track ); + greatTrack.hits().push_back( firstHit ); + greatTrack.hits().push_back( preSelHits[idx_buf[i]] ); + greatTrack.updateX( dx_buf[i], dsl_buf[i] ); + greatTrack.setChi2( chi2_buf[i] ); + + if constexpr ( update_proj ) { + for ( auto& hit : greatTrack.hits() ) { + const float dist = greatTrack.distance( hit ); + hit.projection = dist; + } + } + } + } + +private: + const PrDownTrack& track; + const LHCb::LinAlg::Vec& fit_res; +}; diff --git a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp index 3f1002e50c49c1bd8d7848f817b08c0bd2f3996b..d3dac3208363b13c6ef7b9978c062f07ea3893db 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -69,6 +69,10 @@ class PrLongLivedTracking using SeedTracks = LHCb::Pr::Seeding::Tracks; using DownstreamTracks = LHCb::Pr::Downstream::Tracks; + using simd_t = SIMDWrapper::best::types; + using float_v = simd_t::float_v; + using mask_v = simd_t::mask_v; + public: // - InputLocation: Input location of seed tracks // - OutputLocation: Output location of downstream tracks. @@ -102,18 +106,18 @@ public: const LHCb::UTDAQ::GeomCache& geometry, const DeMagnet& magnet ) const override { // create my state holding all needed mutable variables. std::array preSelHits; - Downstream::Hits matchingXHits; - Downstream::Hits uHitsTemp; + Downstream::SOA::Hits matchingXHits; + Downstream::SOA::Hits uHitsTemp; // -- track collections PrDownTracks goodXTracks; PrDownTracks goodXUTracks; PrDownTracks trackCandidates; - matchingXHits.reserve( 64 ); - trackCandidates.reserve( 16 ); - goodXTracks.reserve( 8 ); - goodXUTracks.reserve( 8 ); - uHitsTemp.reserve( 64 ); + matchingXHits.reserve( 128 ); + trackCandidates.reserve( 128 ); + goodXTracks.reserve( 10 ); + goodXUTracks.reserve( 10 ); + uHitsTemp.reserve( 80 ); //========================================================================== // Get the output container @@ -142,8 +146,8 @@ public: 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 ) ); + PrDownTrack refTrack( stateVector, state.z().cast(), m_zUT.value(), m_zMagnetParams.value(), m_yParams.value(), + m_momentumParams.value(), magScaleFactor * ( -1 ) ); // -- Veto particles coming from the beam pipe. if ( insideBeampipe( refTrack ) ) continue; @@ -172,8 +176,8 @@ public: nXTrack += createTrackCandidates( trackCandidates, refTrack, preSelHits, matchingXHits, uHitsTemp, goodXTracks, goodXUTracks, { 0, 3, 1, 2 }, false ); - bool has4LayerTrack = false; - for ( PrDownTrack& track : trackCandidates ) { has4LayerTrack |= track.hits().size() == 4; } + const bool has4LayerTrack = std::any_of( trackCandidates.begin(), trackCandidates.end(), + []( const PrDownTrack& track ) { return track.hits().size() == 4; } ); if ( !has4LayerTrack ) { nXTrack += createTrackCandidates( trackCandidates, refTrack, preSelHits, matchingXHits, uHitsTemp, goodXTracks, @@ -190,15 +194,15 @@ public: auto bestCandidateScore = std::numeric_limits::infinity(); std::vector downstream_tracks_mlp_datas; - downstream_tracks_mlp_datas.reserve( 16 ); + downstream_tracks_mlp_datas.reserve( trackCandidates.size() ); 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 || - track.hits().size() > m_MaxNumUTHits ) + if ( track.chi2() > m_maxChi2.value() || insideBeampipe( track ) || + track.hits().size() < m_MinNumUTHits.value() || track.hits().size() > m_MaxNumUTHits.value() ) continue; downstream_tracks_mlp_datas.emplace_back( track, _tr ); @@ -211,7 +215,8 @@ public: } ); for ( MLP::PrLongLivedTracking::DataType_t& track_data : downstream_tracks_mlp_datas ) { - if ( bestCandidateScore > track_data.ghost_probability && track_data.ghost_probability < m_maxGhostProb ) { + if ( bestCandidateScore > track_data.ghost_probability && + track_data.ghost_probability < m_maxGhostProb.value() ) { bestCandidate = track_data.m_downstream_track; bestCandidateScore = track_data.ghost_probability; } @@ -289,7 +294,7 @@ private: // -- Parameter to reject seed track which are likely ghosts // - FisherCut: Cut on Fisher-discriminant to reject bad seed tracks. - Gaudi::Property m_seedCut{ this, "FisherCut", -1.0 }; + // Gaudi::Property m_seedCut{ this, "FisherCut", -1.0 }; // -- Parameters for the cut on deltaP (momentum estimate from Seeding and Downstream kink) // - MaxDeltaPConst: Window for deltaP is: MaxDeltaPConst/p + MaxDeltaPOffset @@ -305,9 +310,9 @@ private: // XCorrestionOffset Gaudi::Property m_xCorrectionOffset{ this, "XCorrestionOffset", 0.4 }; // - MaxXTracks: Maximum number of x-tracklets to process further - Gaudi::Property m_maxXTracks{ this, "MaxXTracks", 2 }; + Gaudi::Property m_maxXTracks{ this, "MaxXTracks", 3 }; // - MaxChi2DistXTracks: Maximum chi2 difference to x-tracklet with best chi2 - Gaudi::Property m_maxChi2DistXTracks{ this, "MaxChi2DistXTracks", 0.2 }; + // Gaudi::Property m_maxChi2DistXTracks{ this, "MaxChi2DistXTracks", 0.2 }; // - MaxXUTracks: Maximum number of xu-tracklets to process further Gaudi::Property m_maxXUTracks{ this, "MaxXUTracks", 3 }; Gaudi::Property m_fitXProjChi2Offset{ this, "FitXProjChi2Offset", 4.5 }; @@ -323,11 +328,11 @@ private: // If false only hits and T-tracks from good longtracks are removed. // The criterion for this is the Chi2 of the longtracks from the fit. // - RemoveUsed: Remove seed tracks and used UT hits (with chi2-cut on long track)? - Gaudi::Property m_removeUsed{ this, "RemoveUsed", false }; + // Gaudi::Property m_removeUsed{ this, "RemoveUsed", false }; // - RemoveAll: Remove seed tracks and used UT hits (withoug chi2-cut on long track)? - Gaudi::Property m_removeAll{ this, "RemoveAll", false }; + // Gaudi::Property m_removeAll{ this, "RemoveAll", false }; // - LongChi2: Chi2-cut for the removal - Gaudi::Property m_longChi2{ this, "LongChi2", 1.5 }; + // Gaudi::Property m_longChi2{ this, "LongChi2", 1.5 }; // properties Gaudi::Property m_weightsfilename{ this, "WeightsFileName", @@ -345,9 +350,10 @@ private: // void prepareSeeds(const Tracks& inTracks, std::vector& myInTracks)const; ///< Tag already used T-Seeds int createTrackCandidates( PrDownTracks& trackCandidates, const PrDownTrack& refTrack, - std::array& preSelHits, Downstream::Hits matchingXHits, - Downstream::Hits uHitsTemp, PrDownTracks goodXTracks, PrDownTracks goodXUTracks, + std::array& preSelHits, Downstream::SOA::Hits& matchingXHits, + Downstream::SOA::Hits& uHitsTemp, PrDownTracks& goodXTracks, PrDownTracks& goodXUTracks, const std::array plane_indx, bool skip_second_layer = false ) const { + int nXTrack = 0; for ( auto& myHit : preSelHits[plane_indx[0]] ) { @@ -360,11 +366,11 @@ 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; + if ( ( std::abs( track.momentum() ) < m_initialMinMomentum.value() ) || ( track.pt() < m_initialMinPt.value() ) ) + continue; // -- Fit x projection - findMatchingHits( track, preSelHits[plane_indx[1]], matchingXHits ); - fitXProjection( track, myHit, matchingXHits, goodXTracks, skip_second_layer ); + createXTracks( track, myHit, preSelHits[plane_indx[1]], matchingXHits, goodXTracks, skip_second_layer ); nXTrack += goodXTracks.size(); @@ -379,9 +385,9 @@ private: for ( PrDownTrack& xuTrack : goodXUTracks ) { addVHits( xuTrack, preSelHits[plane_indx[3]] ); if ( xuTrack.hits().size() < 3 ) continue; - simplyFit( xuTrack ); + simplyFit( xuTrack ); // fitAndRemove( xuTrack ); - if ( xuTrack.chi2() > m_maxChi2 || !xuTrack.isYCompatible( m_yTol ) ) continue; + if ( xuTrack.chi2() > m_maxChi2.value() || !xuTrack.isYCompatible( m_yTol.value() ) ) continue; trackCandidates.push_back( std::move( xuTrack ) ); } // Loop over good xu tracks } // Loop over good x tracks @@ -396,11 +402,11 @@ private: const LHCb::Pr::Hits& hitHandler, const LHCb::UTDAQ::GeomCache& geom ) const { // - Max Pt around 100 MeV for strange particle decay -> maximum displacement is in 1/p. - float xPredTol = m_xPredTolOffset; + float xPredTol = m_xPredTolOffset.value(); // P dependance + overal tol. auto p = std::abs( track.momentum() ); - if ( p > 1e-6 ) xPredTol = m_xPredTolConst / p + m_xPredTolOffset; + if ( p > 1e-6 ) xPredTol = m_xPredTolConst.value() / p + m_xPredTolOffset.value(); const float yTol = xPredTol / 2.0 + 7.5; // this is a little vague and not the final word // -- a correction turns out to be beneficial @@ -464,14 +470,39 @@ 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 ) ); + preSelHits[layer].emplace_back( &hitHandler, itHit, xx, mH.get().cast(), + mH.get().cast(), zLayer, fabs( xx - pos ) ); } } } } - - std::sort( preSelHits[1].begin(), preSelHits[1].end(), Downstream::IncreaseByProj ); - std::sort( preSelHits[2].begin(), preSelHits[2].end(), Downstream::IncreaseByProj ); + } + template + [[gnu::always_inline]] inline void fillFitMatrices2D( LHCb::LinAlg::MatSym& mat, + LHCb::LinAlg::Vec& rhs, const HitProxy& hit, + const TrackProxy& track ) const { + const auto dz = 0.001 * ( hit.z0() - track.zMagnet() ); + const auto w = hit.weight(); + const auto t = hit.sin(); + + mat( 0, 0 ) += w; + mat( 1, 0 ) += w * dz; + mat( 1, 1 ) += w * dz * dz; + mat( 2, 0 ) += w * t; + mat( 2, 1 ) += w * dz * t; + mat( 2, 2 ) += w * t * t; + + if constexpr ( useProj ) { + const auto dist = hit.projection; + rhs( 0 ) += w * dist; + rhs( 1 ) += w * dist * dz; + rhs( 2 ) += w * dist * t; + } else { + const auto dist = track.distance( hit ); + rhs( 0 ) += w * dist; + rhs( 1 ) += w * dist * dz; + rhs( 2 ) += w * dist * t; + } } //========================================================================= @@ -487,35 +518,15 @@ private: again = false; //== 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.; + auto mat = LHCb::LinAlg::initialize_with_zeros>(); + auto rhs = LHCb::LinAlg::initialize_with_zeros>(); + mat( 0, 0 ) = track.weightXMag(); + rhs( 0 ) = track.dxMagnet() * track.weightXMag(); std::array differentPlanes = { 0, 0, 0, 0 }; for ( auto& hit : track.hits() ) { - const float dz = 0.001 * ( hit.z - track.zMagnet() ); - const float dist = track.distance( hit ); - 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; - rhs[0] += w * dist; - rhs[1] += w * dist * dz; - rhs[2] += w * dist * t; - + fillFitMatrices2D( mat, rhs, hit, track ); // -- check how many different layers have fired differentPlanes[hit.planeCode()]++; } @@ -525,17 +536,16 @@ private: const std::uint8_t nbUV = differentPlanes[1] + differentPlanes[2]; // -- solve the equation and update the parameters of the track - CholeskyDecomp decomp( mat ); - if ( !decomp ) { + auto [success, inv_M] = mat.invChol(); + if ( !success ) { track.setChi2( 1e42 ); return; - } else { - decomp.Solve( rhs ); } + const auto res = inv_M * rhs; - const float dx = rhs[0]; - const float dsl = 0.001 * rhs[1]; - const float dy = rhs[2]; + const float dx = res( 0 ); + const float dsl = 0.001 * res( 1 ); + const float dy = res( 2 ); if ( nbUV < 4 ) track.updateX( dx, dsl ); track.updateY( dy ); @@ -550,7 +560,7 @@ private: Downstream::Hit& hit = *itH; const float yTrackAtZ = track.yAtZ( hit.z ); - if ( !hit.isYCompatible( yTrackAtZ, m_yTol ) ) { + if ( !hit.isYCompatible( yTrackAtZ, m_yTol.value() ) ) { track.hits().erase( itH ); if ( 2 < track.hits().size() ) again = true; break; @@ -571,7 +581,7 @@ private: if ( track.hits().size() > 2 ) chi2 /= ( track.hits().size() - 2 ); track.setChi2( chi2 ); - if ( m_maxChi2 < chi2 && track.hits().size() > 3 && maxDist > 0 ) { + if ( m_maxChi2.value() < chi2 && track.hits().size() > 3 && maxDist > 0 ) { track.hits().erase( worst ); again = true; } @@ -580,59 +590,28 @@ private: //========================================================================= // 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.; + auto mat = LHCb::LinAlg::initialize_with_zeros>(); + auto rhs = LHCb::LinAlg::initialize_with_zeros>(); + mat( 0, 0 ) = track.weightXMag(); + rhs( 0 ) = track.dxMagnet() * track.weightXMag(); // 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; - } - } + for ( const auto& hit : track.hits() ) { fillFitMatrices2D( mat, rhs, hit, track ); } // -- solve the equation and update the parameters of the track - CholeskyDecomp decomp( mat ); - if ( !decomp ) { + auto [success, inv_M] = mat.invChol(); + if ( !success ) { track.setChi2( 1e42 ); return; - } else { - decomp.Solve( rhs ); } + const auto res = inv_M * rhs; - const float dx = rhs[0]; - const float dsl = 0.001 * rhs[1]; - const float dy = rhs[2]; + const float dx = res( 0 ); + const float dsl = 0.001 * res( 1 ); + const float dy = res( 2 ); track.updateX( dx, dsl ); track.updateY( dy ); @@ -651,26 +630,141 @@ private: } //========================================================================= - // Collect the hits in the other x layer + // Simplified fit function that only fits and does not perform outlier removal //========================================================================= - void findMatchingHits( const PrDownTrack& track, const Downstream::Hits& preSelHits, - Downstream::Hits& matchingXHits ) const { + void simplyFitSIMD( PrDownTrack& track, const Downstream::SOA::Hits& hits_soa, const Downstream::Hits& preSelHits, + PrDownTracks& output, float maxChi2 ) const { + output.clear(); + + auto mat = LHCb::LinAlg::initialize_with_zeros>(); + auto rhs = LHCb::LinAlg::initialize_with_zeros>(); + mat( 0, 0 ) = track.weightXMag(); + rhs( 0 ) = track.dxMagnet() * track.weightXMag(); + + // Fill the matrices using the already stored hits in a track object + // This part is the same for all the U-hits + for ( const auto& hit : track.hits() ) { fillFitMatrices2D( mat, rhs, hit, track ); } + + // Add the U hit and do magic + for ( const auto& hit : hits_soa.simd() ) { + + auto mat_int = mat; + auto rhs_int = rhs; + + fillFitMatrices2D( mat_int, rhs_int, hit, track ); + + auto [success, inv_M] = mat_int.invChol(); + if ( !any( success ) ) { continue; }; + + auto res = inv_M * rhs_int; + res( 1 ) *= 0.001; + + PrIncompleteDownTrack test_track( track, res ); + const auto chi2 = test_track.chi2( hit, m_yTol.value() ); + + // Check the constraints + const auto hit_ok = ( success & ( chi2 < maxChi2 ) & hit.loop_mask() ); + if ( !any( hit_ok ) ) { continue; } + + test_track.compressstore( output, hit_ok, chi2, hit.idx(), preSelHits ); + } + } + + /** + * @brief Creates X tracks by combining two X-hits and performing a fit. + * + * This function attempts to extend a given track with a given firstHit by matching it with hits in the second X + * layer. It first selects hits within a matching window, then performs a fit using SIMD instructions. Tracks passing + * the fit quality criteria are added to the output collection. + * + * @param track The downstream track to be extended. + * @param firstHit The first hit associated with the track. + * @param preSelHits Preselected hits in the X layers to be considered for matching. + * @param matchingXHits Buffer container for hits that match the seed track within the matching window. + * @param goodXTracks Output container for successfully fitted X tracks. + * @param skip_second_layer If true, forces the tracks to have only one X hit, effectively skipping the second layer. + * + */ + void createXTracks( const PrDownTrack& track, const Downstream::Hit& firstHit, const Downstream::Hits& preSelHits, + Downstream::SOA::Hits& matchingXHits, PrDownTracks& goodXTracks, + const bool skip_second_layer ) const { + + // Prepare the buffers + goodXTracks.clear(); matchingXHits.clear(); - if ( preSelHits.empty() ) return; - const float tol = - std::min( m_maxWindow.value(), m_tolMatchOffset + m_tolMatchConst * std::abs( track.stateQoP() ) ); + if ( preSelHits.empty() ) { + // 3-hit case handling + auto& tr = goodXTracks.emplace_back( track ); + tr.hits().push_back( firstHit ); + return; + }; + + const float tol = std::min( m_maxWindow.value(), + m_tolMatchOffset.value() + m_tolMatchConst.value() * std::abs( track.stateQoP() ) ); + const float maxChi2 = m_fitXProjChi2Offset.value() + m_fitXProjChi2Const.value() / std::abs( track.momentum() ); + const float xPred = track.xAtZ( preSelHits.front().z ); + // Preselect hits for ( auto& hit : preSelHits ) { const float adist = std::abs( hit.x - xPred ); - if ( adist <= tol ) { matchingXHits.push_back( hit ); } + if ( adist <= tol ) { + auto newhit = matchingXHits.emplace_back(); + newhit.copyFromScalarHit( hit, std::distance( preSelHits.data(), &hit ) ); + } + } + + // Perform the fit + const auto w1 = firstHit.weight(); + const auto d1 = track.distance( firstHit ); + const auto dz1 = firstHit.z - track.zMagnet(); + + for ( const auto& secondHit : matchingXHits.simd() ) { + + const auto w2 = secondHit.weight(); + const auto d2 = track.distance( secondHit ); + const auto dz2 = secondHit.z0() - track.zMagnet(); + + constexpr float scale = 0.0009765625f; // = 1/1024; just to make it a power of 2 + const auto mat = + std::array{ scale * ( track.weightXMag() + w1 + w2 ), scale * scale * ( w1 * dz1 + w2 * dz2 ), + scale * scale * scale * ( w1 * dz1 * dz1 + w2 * dz2 * dz2 ) }; + const auto rhs = + std::array{ scale * ( w1 * d1 + w2 * d2 ), scale * scale * ( w1 * d1 * dz1 + w2 * d2 * dz2 ) }; + + // Solve linear system + const auto det = mat[0] * mat[2] - mat[1] * mat[1]; + LHCb::LinAlg::Vec fit_res; + fit_res( 0 ) = ( mat[2] * rhs[0] - mat[1] * rhs[1] ) / det; + fit_res( 1 ) = scale * ( mat[0] * rhs[1] - mat[1] * rhs[0] ) / det; + + PrIncompleteDownTrack test_track( track, fit_res ); + const auto chi2 = test_track.chi2( firstHit, secondHit ); + + const auto hit_ok = ( chi2 < maxChi2 ) & secondHit.loop_mask(); + if ( none( hit_ok ) ) continue; + // If we found smth during the second pass we should reject everything, + // as those tracks are supposed to be reconstructed during the first pass + if ( any( hit_ok ) && skip_second_layer ) return; + + test_track.compressstore( goodXTracks, hit_ok, chi2, secondHit.idx(), firstHit, preSelHits ); } - std::sort( matchingXHits.begin(), matchingXHits.end(), - [xPred]( const Downstream::Hit& lhs, const Downstream::Hit& rhs ) { - return std::abs( lhs.x - xPred ) < std::abs( rhs.x - xPred ); - } ); + if ( goodXTracks.size() > m_maxXTracks.value() ) { + // We have more than m_maxXTracks, sort by the distance from the track to a new hit and keep only the best ones + std::partial_sort( goodXTracks.begin(), goodXTracks.begin() + m_maxXTracks.value(), goodXTracks.end(), + [&]( const PrDownTrack& lhs, const PrDownTrack& rhs ) { + return std::abs( lhs.hits().back().x - xPred ) < std::abs( rhs.hits().back().x - xPred ); + } ); + goodXTracks.erase( goodXTracks.begin() + m_maxXTracks.value(), goodXTracks.end() ); + } + + if ( goodXTracks.empty() ) { + // 3-hit case handling + auto& tr = goodXTracks.emplace_back( track ); + tr.hits().push_back( firstHit ); + } } //========================================================================= @@ -701,9 +795,9 @@ private: } // Create a state at zUTa - newTrack.field().setPosition( static_cast( track.xAtZ( m_zUTa ) ), - static_cast( track.yAtZ( m_zUTa ) ), - static_cast( m_zUTa ) ); + newTrack.field().setPosition( static_cast( track.xAtZ( m_zUTa.value() ) ), + static_cast( track.yAtZ( m_zUTa.value() ) ), + static_cast( m_zUTa.value() ) ); newTrack.field().setDirection( static_cast( track.slopeX() ), static_cast( track.slopeY() ) ); newTrack.field().setQOverP( static_cast( 1.0 / track.momentum() ) ); @@ -712,7 +806,7 @@ private: //========================================================================= // Add the U hits. //========================================================================= - void addUHits( PrDownTrack& track, const Downstream::Hits& preSelHits, Downstream::Hits& uHitsTemp, + void addUHits( PrDownTrack& track, const Downstream::Hits& preSelHits, Downstream::SOA::Hits& uHitsTemp, PrDownTracks& goodXUTracks ) const { goodXUTracks.clear(); @@ -720,7 +814,7 @@ private: uHitsTemp.clear(); - const float tol = m_tolUOffset + m_tolUConst / std::abs( track.momentum() ); + const float tol = m_tolUOffset.value() + m_tolUConst.value() / std::abs( track.momentum() ); // -- these numbers are a little arbitrary float minChi2 = ( track.hits().size() == 1 ) ? 800 : 300; @@ -731,41 +825,29 @@ private: 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::IncreaseByAbsProj ); + auto newhit = uHitsTemp.emplace_back(); + newhit.copyFromScalarHit( hit, std::distance( preSelHits.data(), &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 ); - simplyFit( greatTrack ); - - // -- it's sorted - if ( greatTrack.chi2() > minChi2 ) { - goodXUTracks.pop_back(); - break; - } - - if ( !greatTrack.isYCompatible( m_yTol ) ) { - goodXUTracks.pop_back(); - continue; - } + simplyFitSIMD( track, uHitsTemp, preSelHits, goodXUTracks, minChi2 ); - if ( goodXUTracks.size() >= m_maxXUTracks ) { break; } + // If we have more than m_maxXUTracks, sort by chi2 and keep only the best ones + constexpr auto IncreaseByChi2 = []( const PrDownTrack& lhs, const PrDownTrack& rhs ) { + return lhs.chi2() < rhs.chi2(); + }; + if ( goodXUTracks.size() > m_maxXUTracks.value() ) { + std::partial_sort( goodXUTracks.begin(), goodXUTracks.begin() + m_maxXUTracks.value(), goodXUTracks.end(), + IncreaseByChi2 ); + goodXUTracks.erase( goodXUTracks.begin() + m_maxXUTracks.value(), goodXUTracks.end() ); } // 3-hit case handling - if ( ( goodXUTracks.size() == 0 ) && ( track.hits().size() == 2 ) ) { goodXUTracks.emplace_back( track ); } + if ( ( goodXUTracks.empty() ) && ( track.hits().size() == 2 ) ) { goodXUTracks.emplace_back( track ); } } //========================================================================= @@ -775,7 +857,8 @@ private: if ( preSelHits.empty() ) { return; } const auto p = std::abs( track.momentum() ); - const float tol = ( track.hits().size() == 2 ) ? m_tolUOffset + m_tolUConst / p : m_tolVOffset + m_tolVConst / p; + const float tol = ( track.hits().size() == 2 ) ? m_tolUOffset.value() + m_tolUConst.value() / p + : m_tolVOffset.value() + m_tolVConst.value() / p; // Find the closest hit to the track // It's enough to look at the closest hit only, since absolute majority of the tracks @@ -802,51 +885,12 @@ private: if ( bestHit ) { track.hits().push_back( *bestHit ); track.hits().back().projection = bestDist; - simplyFit( track ); + simplyFit( track ); // fitAndRemove( track ); } track.sortFinalHits(); } - // void tagUsedUT( const Track* tr ) const; ///< Tag hits that were already used elsewhere - - //============================================================================= - // Fit the projection in the zx plane, one hit in each x layer - //============================================================================= - void fitXProjection( const PrDownTrack& track, const Downstream::Hit& firstHit, const Downstream::Hits& matchingXHits, - PrDownTracks& goodXTracks, const bool skip_second_layer ) const { - goodXTracks.clear(); - - const float maxChi2 = m_fitXProjChi2Offset + m_fitXProjChi2Const / std::abs( track.momentum() ); - - // Catch if there is no second hit in other station - for ( const auto& hit : matchingXHits ) { - auto& tr = goodXTracks.emplace_back( track ); - xFit( tr, firstHit, hit ); - - if ( tr.chi2() > maxChi2 ) { - goodXTracks.pop_back(); - continue; - } - - tr.hits().push_back( firstHit ); - tr.hits().push_back( hit ); - - if ( goodXTracks.size() >= 3 ) { break; } - } - - if ( skip_second_layer && goodXTracks.size() != 0 ) { - goodXTracks.clear(); - return; - } - - // 3-hit case handling - if ( goodXTracks.size() == 0 ) { - auto& tr = goodXTracks.emplace_back( track ); - tr.hits().push_back( firstHit ); - } - } - //========================================================================= // Check if the new candidate is better than the old one //========================================================================= @@ -858,7 +902,7 @@ private: // -- use a tighter chi2 for 3 hit tracks // -- as they are more likely to be ghosts - const float maxChi2 = ( nbMeasureOK == 3 ) ? m_maxChi2ThreeHits : m_maxChi2; + const float maxChi2 = ( nbMeasureOK == 3 ) ? m_maxChi2ThreeHits.value() : m_maxChi2.value(); //== Good enough Chi2/ndof if ( maxChi2 < track.chi2() ) { return false; } @@ -868,7 +912,7 @@ private: const float deltaP = p * track.stateQoP() - 1.; if ( maxDeltaP( track ) < fabs( deltaP ) && !magnetOff ) { return false; } - if ( std::abs( p ) < m_minMomentum || track.pt() < m_minPt ) { return false; } + if ( std::abs( p ) < m_minMomentum.value() || track.pt() < m_minPt.value() ) { return false; } return true; } @@ -897,9 +941,10 @@ private: // -- 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 ) ) { + if ( std::abs( hit.z - trackHit.z ) < 1.0 || std::abs( hit.x - trackHit.x ) > 3 * m_overlapTol.value() ) + continue; + if ( m_overlapTol.value() > std::abs( track.distanceLinear( hit, xStart, zStart ) ) && + hit.isYCompatible( track.yAtZ( hit.z ), m_yTol.value() ) ) { // -- preSelHits is not used anymore track.hits().push_back( std::move( hit ) ); hitAdded = true; @@ -926,40 +971,10 @@ private: return getFisher( vals ); } */ - //========================================================================= - // Fit hits in x layers, using the magnet point as constraint. - //========================================================================= - void xFit( PrDownTrack& track, const Downstream::Hit& hit1, const Downstream::Hit& hit2 ) const { - - const auto w1 = hit1.weight(); - const auto w2 = hit2.weight(); - - auto d1 = track.distance( hit1 ); - auto d2 = track.distance( hit2 ); - - const auto dz1 = hit1.z - track.zMagnet(); - const auto dz2 = hit2.z - track.zMagnet(); - - auto mat = std::array{ track.weightXMag() + w1 + w2, w1 * dz1 + w2 * dz2, w1 * dz1 * dz1 + w2 * dz2 * dz2 }; - auto rhs = std::array{ w1 * d1 + w2 * d2, w1 * d1 * dz1 + w2 * d2 * dz2 }; - - // Solve linear system - const float det = mat[0] * mat[2] - mat[1] * mat[1]; - const float dx = ( mat[2] * rhs[0] - mat[1] * rhs[1] ) / det; - const float dsl = ( mat[0] * rhs[1] - mat[1] * rhs[0] ) / det; - - track.updateX( dx, dsl ); - - d1 = track.distance( hit1 ); - d2 = track.distance( hit2 ); - - const float chi2 = track.initialChi2() + w1 * d1 * d1 + w2 * d2 * d2; - track.setChi2( chi2 ); - } - /// Does this track point inside the beampipe? bool insideBeampipe( const PrDownTrack& track ) const { - return ( m_minUTx > fabs( track.xAtZ( m_zUTa ) ) ) && ( m_minUTy > fabs( track.yAtZ( m_zUTa ) ) ); + return ( m_minUTx.value() > fabs( track.xAtZ( m_zUTa.value() ) ) ) && + ( m_minUTy.value() > fabs( track.yAtZ( m_zUTa.value() ) ) ); } /// Helper to evaluate the Fisher discriminant @@ -974,13 +989,13 @@ private: /// Helper to evaluate the maximum discrepancy between momentum from kink and curvature in T-stations float maxDeltaP( const PrDownTrack& track ) const { - return m_maxDeltaPConst / std::abs( track.momentum() ) + m_maxDeltaPOffset; + return m_maxDeltaPConst.value() / std::abs( track.momentum() ) + m_maxDeltaPOffset.value(); } /// Helper to evaluate the correction to the x position in the UT float xPosCorrection( const PrDownTrack& track ) const { auto p = track.momentum(); - return std::copysign( m_xCorrectionOffset + m_xCorrectionConst / std::abs( p ), p ); + return std::copysign( m_xCorrectionOffset.value() + m_xCorrectionConst.value() / std::abs( p ), p ); } // -- counters