diff --git a/Pr/PrAlgorithms/src/PrDownTrack.h b/Pr/PrAlgorithms/src/PrDownTrack.h index 8cb7dcbe506d8d642b109bf86b6c75cc09d5d7a8..4a4391bd7299d47cfabd01691ae083c0c56066fd 100644 --- a/Pr/PrAlgorithms/src/PrDownTrack.h +++ b/Pr/PrAlgorithms/src/PrDownTrack.h @@ -209,7 +209,8 @@ public: PrDownTrack() = default; PrDownTrack( Gaudi::TrackVectorF stateVector, float stateZ, float zUT, LHCb::span magnetParams, - LHCb::span yParams, LHCb::span momPar, float magnetScale ) + LHCb::span yParams, LHCb::span momPar, float magnetScale, + bool useYSlopefromSeed ) : m_stateVector( stateVector ), m_stateZ( stateZ ), m_zUT( zUT ) { const auto tx2 = stateTx() * stateTx(); @@ -228,8 +229,10 @@ public: const float dSlope = std::abs( m_slopeX - stateTx() ); const float dSlope2 = dSlope * dSlope; - const float by = stateY() / ( stateZ + ( yParams[0] * fabs( stateTy() ) * zMagnet + yParams[1] ) * dSlope2 ); - m_slopeY = by * ( 1. + yParams[0] * fabs( by ) * dSlope2 ); + const float by = useYSlopefromSeed + ? stateTy() + : stateY() / ( stateZ + ( yParams[0] * fabs( stateTy() ) * zMagnet + yParams[1] ) * dSlope2 ); + m_slopeY = useYSlopefromSeed ? stateTy() : by * ( 1. + yParams[0] * fabs( by ) * dSlope2 ); const float yMagnet = stateY() + dz * by - yParams[1] * by * dSlope2; @@ -250,7 +253,7 @@ public: m_weightXMag = 1.0 / ( errXMag * errXMag ); m_weightYMag = 1.0 / ( errYMag * errYMag ); - m_magnet = Gaudi::XYZPoint( xMagnet, yMagnet, zMagnet ); + m_magnet = Gaudi::XYZPointF( xMagnet, yMagnet, zMagnet ); //=== Save for reference m_displX = 0.; diff --git a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp index d3dac3208363b13c6ef7b9978c062f07ea3893db..55198b4da177ba28e8b57d12395e4e0062bb3d38 100644 --- a/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp +++ b/Pr/PrAlgorithms/src/PrLongLivedTracking.cpp @@ -122,12 +122,16 @@ public: //========================================================================== // Get the output container //========================================================================== - DownstreamTracks finalTracks{ &InputTracks, Zipping::generateZipIdentifier(), LHCb::getMemResource( evtCtx ) }; + DownstreamTracks finalTracks{ &InputTracks, + m_highMassTuning.value() ? DownstreamTracks::History::PrDownstreamHighMass + : DownstreamTracks::History::PrDownstream, + Zipping::generateZipIdentifier(), LHCb::getMemResource( evtCtx ) }; + finalTracks.reserve( InputTracks.size() ); const float magScaleFactor = magnet.signedRelativeCurrent(); - bool magnetOff = std::abs( magScaleFactor ) > 1e-6 ? false : true; + const bool magnetOff = std::abs( magScaleFactor ) < 1e-6; m_nSeeds += InputTracks.size(); @@ -147,7 +151,7 @@ public: 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.value(), m_zMagnetParams.value(), m_yParams.value(), - m_momentumParams.value(), magScaleFactor * ( -1 ) ); + m_momentumParams.value(), magScaleFactor * ( -1 ), m_highMassTuning.value() ); // -- Veto particles coming from the beam pipe. if ( insideBeampipe( refTrack ) ) continue; @@ -196,13 +200,17 @@ public: std::vector downstream_tracks_mlp_datas; downstream_tracks_mlp_datas.reserve( trackCandidates.size() ); + std::sort( preSelHits[0].begin(), preSelHits[0].end(), Downstream::IncreaseByProj ); + std::sort( preSelHits[3].begin(), preSelHits[3].end(), Downstream::IncreaseByProj ); + for ( PrDownTrack& track : trackCandidates ) { if ( has4LayerTrack && track.hits().size() < 4 ) continue; addOverlapRegions( track, preSelHits ); - if ( track.chi2() > m_maxChi2.value() || insideBeampipe( track ) || - track.hits().size() < m_MinNumUTHits.value() || track.hits().size() > m_MaxNumUTHits.value() ) + if ( ( track.chi2() > m_maxChi2.value() || insideBeampipe( track ) || + track.hits().size() < m_MinNumUTHits.value() || track.hits().size() > m_MaxNumUTHits.value() ) || + ( m_highMassTuning.value() && track.pt() < m_minPt.value() ) ) continue; downstream_tracks_mlp_datas.emplace_back( track, _tr ); @@ -236,11 +244,17 @@ private: // DataObjectReadHandle m_ForwardTracks { this, "ForwardTracks" LHCb::TrackLocation::Forward }; // DataObjectReadHandle m_MatchTracks { this, "MatchTracks", LHCb::TrackLocation::Match }; // + Gaudi::Property m_highMassTuning{ this, "UseHighMassTunings", false }; // - XPredTolConst: x-window for preselection is XPredTolConst/p + XPredTolOffset Gaudi::Property m_xPredTolConst{ this, "XPredTolConst", 200. * Gaudi::Units::mm* Gaudi::Units::GeV }; // - XPredTolOffset: x-window for preselection is XPredTolConst/p + XPredTolOffset Gaudi::Property m_xPredTolOffset{ this, "XPredTolOffset", 6. * Gaudi::Units::mm }; + // - YPredTolConst: y-window for preselection is XPredTol/YPredTolConst + YPredTolOffset + Gaudi::Property m_yPredTolConst{ this, "YPredTolConst", 2.0 }; + // - YPredTolOffset: y-window for preselection is XPredTol/YPredTolConst + YPredTolOffset + Gaudi::Property m_yPredTolOffset{ this, "YPredTolOffset", 7.5 * Gaudi::Units::mm }; + // - TolMatchConst: x-window for matching x hits is TolMatchConst/p + TolMatchOffset Gaudi::Property m_tolMatchConst{ this, "TolMatchConst", 20000. }; // - TolMatchOffset: x-window for matching x hits is TolMatchConst/p + TolMatchOffset @@ -264,7 +278,7 @@ private: // - MinUTy: half-height of of beampipe rectangle Gaudi::Property m_minUTy{ this, "MinUTy", 25. * Gaudi::Units::mm }; // - MaxGhostProb: Maximum ghost probability for tracks - Gaudi::Property m_maxGhostProb{ this, "MaxGhostProb", 0.75 }; + Gaudi::Property m_maxGhostProb{ this, "GhostProbCut", 0.75 }; // Define parameters for MC09 field, zState = 9410 // - ZMagnetParams: Parameters to determine the z-position of the magnet point. Tune with PrKsParams. @@ -279,7 +293,7 @@ private: // - ZUTa: z-position of first UT station Gaudi::Property m_zUTa{ this, "ZUTa", 2350. * Gaudi::Units::mm }; // - InitialMinPt: Minimum pT of the track from initial estimate - Gaudi::Property m_initialMinPt{ this, "InitialMinPt", 0. * Gaudi::Units::MeV }; + Gaudi::Property m_initialMinPt{ this, "InitialMinPt", 20. * Gaudi::Units::MeV }; // - 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 @@ -407,11 +421,11 @@ private: // P dependance + overal tol. auto p = std::abs( track.momentum() ); 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 + const float yTol = xPredTol / m_yPredTolConst.value() + m_yPredTolOffset.value(); // -- a correction turns out to be beneficial // -- maybe to compensate tracks not coming from 0/0 (?) - const float correction = xPosCorrection( track ); + const float correction = m_highMassTuning.value() ? 0.0f : xPosCorrection( track ); for ( auto& i : preSelHits ) i.clear(); @@ -471,12 +485,18 @@ private: if ( xPredTol < std::abs( pos - xx ) ) continue; // go from -x to +x preSelHits[layer].emplace_back( &hitHandler, itHit, xx, mH.get().cast(), - mH.get().cast(), zLayer, fabs( xx - pos ) ); + mH.get().cast(), zLayer, xx - pos ); } } } } + + std::sort( preSelHits[1].begin(), preSelHits[1].end(), Downstream::IncreaseByProj ); + std::sort( preSelHits[2].begin(), preSelHits[2].end(), Downstream::IncreaseByProj ); } + //======================================================================== + // Fill matrices + //========================================================================= template [[gnu::always_inline]] inline void fillFitMatrices2D( LHCb::LinAlg::MatSym& mat, LHCb::LinAlg::Vec& rhs, const HitProxy& hit, @@ -598,7 +618,6 @@ private: 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() ) { fillFitMatrices2D( mat, rhs, hit, track ); } // -- solve the equation and update the parameters of the track @@ -706,7 +725,7 @@ private: const float xPred = track.xAtZ( preSelHits.front().z ); - // Preselect hits + // -- The x search seems sensitive to the ordering, so we do not use a binary search here to find the starting point for ( auto& hit : preSelHits ) { const float adist = std::abs( hit.x - xPred ); if ( adist <= tol ) { @@ -824,8 +843,17 @@ private: const float zStart = preSelHits[0].z; const float xStart = track.xAtZ( zStart ); - for ( const auto& hit : preSelHits ) { - const float dist = track.distanceLinear( hit, xStart, zStart ); + // -- binary search to find the starting point. The sorting is not perfect, that's why we add a 1mm tolerance + const auto itBeg = + std::lower_bound( preSelHits.begin(), preSelHits.end(), tol, [&]( const Downstream::Hit& hit, float tol ) { + return track.distanceLinear( hit, xStart, zStart ) < -tol - 1.0; + } ); + + for ( auto it = itBeg; it < preSelHits.end(); ++it ) { + const Downstream::Hit& hit = *it; + const float dist = track.distanceLinear( hit, xStart, zStart ); + + if ( dist > tol + 1.0 ) break; if ( std::abs( dist ) > tol ) continue; auto newhit = uHitsTemp.emplace_back(); @@ -870,11 +898,17 @@ private: const float zStart = preSelHits[0].z; const float xStart = track.xAtZ( zStart ); - for ( auto& hit : preSelHits ) { - // const float dist = track.distance( hit ); + // -- binary search to find the starting point. The sorting is not perfect, that's why we add a 1mm tolerance + const auto itBeg = + std::lower_bound( preSelHits.begin(), preSelHits.end(), tol, [&]( const Downstream::Hit& hit, float tol ) { + return track.distanceLinear( hit, xStart, zStart ) < -tol - 1.0; + } ); + + for ( auto it = itBeg; it < preSelHits.end(); ++it ) { + const auto& hit = *it; const float dist = track.distanceLinear( hit, xStart, zStart ); + if ( dist > tol + 1.0 ) break; if ( std::abs( dist ) > tol ) continue; - // hit.projection = dist; if ( std::abs( dist ) < std::abs( bestDist ) ) { bestDist = dist; bestHit = &hit; @@ -886,37 +920,10 @@ private: track.hits().push_back( *bestHit ); track.hits().back().projection = bestDist; simplyFit( track ); - // fitAndRemove( track ); } track.sortFinalHits(); } - //========================================================================= - // Check if the new candidate is better than the old one - //========================================================================= - bool acceptCandidate( const PrDownTrack& track, bool magnetOff ) const { - const int nbMeasureOK = track.hits().size(); - - //== Enough measures to have Chi2/ndof. - if ( nbMeasureOK < 3 ) { return false; } - - // -- use a tighter chi2 for 3 hit tracks - // -- as they are more likely to be ghosts - const float maxChi2 = ( nbMeasureOK == 3 ) ? m_maxChi2ThreeHits.value() : m_maxChi2.value(); - - //== Good enough Chi2/ndof - if ( maxChi2 < track.chi2() ) { return false; } - - //== Compatible momentum - const float p = track.momentum(); - const float deltaP = p * track.stateQoP() - 1.; - - if ( maxDeltaP( track ) < fabs( deltaP ) && !magnetOff ) { return false; } - if ( std::abs( p ) < m_minMomentum.value() || track.pt() < m_minPt.value() ) { return false; } - - return true; - } - //============================================================================= // This is needed for tracks which have more than one x hit in one layer // Maybe we could make this smarter and do it for every track and add the 'second best' @@ -937,10 +944,17 @@ private: const float zStart = preSelHits[planeCode].front().z; const float xStart = track.xAtZ( zStart ); - for ( const auto& hit : preSelHits[planeCode] ) { + // -- binary search to find the starting point. The sorting is not perfect, that's why we add a 1mm tolerance + const auto itBeg = + std::lower_bound( preSelHits[planeCode].begin(), preSelHits[planeCode].end(), 3 * m_overlapTol.value(), + [&]( const Downstream::Hit& hit, float tol ) { return hit.x - trackHit.x < -tol - 3.0; } ); + + for ( auto it = itBeg; it < preSelHits[planeCode].end(); ++it ) { // -- 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 + const auto& hit = *it; + if ( hit.x - trackHit.x > 3 * m_overlapTol.value() + 3.0 ) break; 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 ) ) && @@ -958,35 +972,12 @@ private: } } - //============================================================================= - // Evaluate the Fisher discriminant for a preselection of seed tracks - //============================================================================= - /* float evaluateFisher( const Track* track ) { - const unsigned int nbIT = std::count_if( track->lhcbIDs().begin(), track->lhcbIDs().end(), - [](const LHCb::LHCbID id){ return id.isIT();}); - float nbITD = static_cast(nbIT); - float lhcbIDSizeD = static_cast(track->lhcbIDs().size()); - std::array vals = { track->chi2PerDoF(), track->p(), track->pt(), nbITD, lhcbIDSizeD }; - - return getFisher( vals ); - } */ - /// Does this track point inside the beampipe? bool insideBeampipe( const PrDownTrack& track ) const { 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 - /* float getFisher(const std::array vals) { - float c_fishConst = -1.69860581797; - std::array c_fishCoefficients = {-0.241020410138, 3.03197732663e-07, -1.14400162824e-05, - 0.126857153245, 0.122359738469}; - float fishVal = c_fishConst; - for(int i = 0; i < 5; i++) fishVal += vals[i] * c_fishCoefficients[i]; - return fishVal; - } */ - /// Helper to evaluate the maximum discrepancy between momentum from kink and curvature in T-stations float maxDeltaP( const PrDownTrack& track ) const { return m_maxDeltaPConst.value() / std::abs( track.momentum() ) + m_maxDeltaPOffset.value(); diff --git a/Pr/PrAlgorithms/src/PrResidualPrUTHits.cpp b/Pr/PrAlgorithms/src/PrResidualPrUTHits.cpp index 4e383fb81beaf5195b2a906b86d7c787c86224d4..7e95106cc155255c4d6da0b19961b27f0700ffcd 100644 --- a/Pr/PrAlgorithms/src/PrResidualPrUTHits.cpp +++ b/Pr/PrAlgorithms/src/PrResidualPrUTHits.cpp @@ -13,6 +13,8 @@ #include "Event/PrHits.h" #include "Event/PrLongTracks.h" #include "Event/PrUpstreamTracks.h" +#include "Event/Track.h" +#include "UTDAQ/UTDAQHelper.h" #include "UTDAQ/UTInfo.h" #include "UTDet/DeUTDetector.h" @@ -48,11 +50,16 @@ public: { typename base_class_t::KeyValue{ "TracksLocation", "" }, typename base_class_t::KeyValue{ "PrUTHitsLocation", "" } }, typename base_class_t::KeyValue{ "PrUTHitsOutput", "" } ) {} + + mutable Gaudi::Accumulators::AveragingCounter m_inHits{ this, "# Hits in input container" }; + mutable Gaudi::Accumulators::AveragingCounter m_outHits{ this, "# Hits in output container" }; + mutable Gaudi::Accumulators::AveragingCounter m_vetoedHits{ this, "# Hits that are vetoed" }; }; // Declaration of the Algorithm Factory DECLARE_COMPONENT_WITH_ID( PrResidualPrUTHits, "PrResidualPrUTHits" ) DECLARE_COMPONENT_WITH_ID( PrResidualPrUTHits, "PrResidualPrUTHits_Upstream" ) +DECLARE_COMPONENT_WITH_ID( PrResidualPrUTHits, "PrResidualPrUTHits_v1Tracks" ) //============================================================================= // Main execution @@ -62,20 +69,40 @@ LHCb::Pr::UT::Hits PrResidualPrUTHits::operator()( const EventContext& evtCtx const LHCb::Pr::UT::Hits& uthithandler ) const { LHCb::Pr::UT::Hits tmp{ Zipping::generateZipIdentifier(), LHCb::getMemResource( evtCtx ) }; - // mark used UT hits const unsigned int nhits = uthithandler.nHits(); + m_inHits += nhits; tmp.reserve( nhits ); boost::dynamic_bitset<> used{ nhits, false }; - /// mark used SciFi Hits - for ( const auto& track : tracks.scalar() ) { - const int nuthits = track.nUTHits().cast(); - for ( int idx = 0; idx < nuthits; idx++ ) { - const int index = track.ut_index( idx ).cast(); - if ( index >= 0 ) used[index] = true; + /// mark used UT hits + if constexpr ( std::is_same_v ) { + for ( const auto& track : tracks ) { + for ( const auto id : track->lhcbIDs() ) { + if ( id.isUT() ) { + const auto sectorID = id.utID().sectorFullID(); + const auto indexs = uthithandler.indices( sectorID ); + for ( int idx = indexs.first; idx != indexs.second; idx++ ) { + if ( uthithandler.lhcbid( idx ) == id ) { + used[idx] = true; + break; + } + } + } + } + } + } else { + /// mark used UT Hits + for ( const auto& track : tracks.scalar() ) { + const int nuthits = track.nUTHits().cast(); + for ( int idx = 0; idx < nuthits; idx++ ) { + const int index = track.ut_index( idx ).cast(); + if ( index >= 0 ) used[index] = true; + } } } + m_vetoedHits += used.count(); + for ( auto fullchan = 0; fullchan < static_cast( UTInfo::MaxNumberOfSectors ); fullchan++ ) { const auto indexs = uthithandler.indices( fullchan ); @@ -84,5 +111,9 @@ LHCb::Pr::UT::Hits PrResidualPrUTHits::operator()( const EventContext& evtCtx tmp.copyHit( fullchan, idx, uthithandler ); } } + m_outHits += tmp.size(); + + assert( tmp.size() + used.count() == nhits && "# output hits + # vetoed hits != # input hits" ); + return tmp; } diff --git a/Pr/PrConverters/src/fromPrTracksV1Track.cpp b/Pr/PrConverters/src/fromPrTracksV1Track.cpp index 9dd844a21733c06b3317ee37574584a756de2738..855fc5342355ffbe78cc3c8abda42b16dcc084fa 100644 --- a/Pr/PrConverters/src/fromPrTracksV1Track.cpp +++ b/Pr/PrConverters/src/fromPrTracksV1Track.cpp @@ -289,7 +289,10 @@ namespace LHCb::Converters::Track::v1 { outTrack->setType( ConversionInfo::Type( inTrack.backward() ) ); else outTrack->setType( ConversionInfo::Type ); - outTrack->setHistory( ConversionInfo::PrHistory ); + if constexpr ( std::is_same_v ) + outTrack->setHistory( inTracks.history() ); + else + outTrack->setHistory( ConversionInfo::PrHistory ); outTracks.insert( outTrack, inTrack.indices().cast() ); } return outTracks; diff --git a/Tr/PrKalmanFilter/include/PrKalmanFilter/KF.h b/Tr/PrKalmanFilter/include/PrKalmanFilter/KF.h index 8620c1db26ceb935005e986caa13ba3ddb7a9f73..1bb137e4c53ab9b9933f18862ca849c790f8b7b4 100644 --- a/Tr/PrKalmanFilter/include/PrKalmanFilter/KF.h +++ b/Tr/PrKalmanFilter/include/PrKalmanFilter/KF.h @@ -1223,7 +1223,7 @@ namespace LHCb::Pr::Tracks::Fit { } else if constexpr ( isDownstream ) { - return std::make_unique( History::PrDownstream, Type::Downstream, PatRecStatus::PatRecIDs ); + return std::make_unique( tracks.history(), Type::Downstream, PatRecStatus::PatRecIDs ); } else if constexpr ( isUpstream ) { @@ -1308,7 +1308,7 @@ namespace LHCb::Pr::Tracks::Fit { if constexpr ( isLong ) return tracks.history(); else if constexpr ( isDownstream ) - return LHCb::Event::Enum::Track::History::PrDownstream; + return tracks.history(); else if constexpr ( isVelo ) return LHCb::Event::Enum::Track::History::PrPixel; else if constexpr ( isSeed ) diff --git a/Tr/TrackUtils/src/TracksDownstreamConverter.cpp b/Tr/TrackUtils/src/TracksDownstreamConverter.cpp index 71e69d3425ef41e57930f0a73ac3d7a1a5a8357f..67dc6b19bb8d3b4804973b30c309d9613d6d80a2 100644 --- a/Tr/TrackUtils/src/TracksDownstreamConverter.cpp +++ b/Tr/TrackUtils/src/TracksDownstreamConverter.cpp @@ -71,7 +71,7 @@ public: for ( auto const lhcbid : track.lhcbIDs() ) { newTrack.addToLhcbIDs( lhcbid ); } newTrack.setType( Track::Type::Downstream ); - newTrack.setHistory( Track::History::PrDownstream ); + newTrack.setHistory( tracksUT.history() ); newTrack.setPatRecStatus( Track::PatRecStatus::PatRecIDs ); newTrack.addToAncestors( trackSeed ); }