From 9194d920f0c967c68dba193cf718e3e940189c25 Mon Sep 17 00:00:00 2001 From: Arthur Hennequin Date: Tue, 6 May 2025 16:35:27 +0200 Subject: [PATCH 1/4] Add velo clone killer algorithm --- Pr/PrPixel/CMakeLists.txt | 1 + Pr/PrPixel/src/VeloCloneKiller.cpp | 126 ++++++++++++++++++ .../src/TrackVPOverlapMonitor.cpp | 9 +- 3 files changed, 133 insertions(+), 3 deletions(-) create mode 100644 Pr/PrPixel/src/VeloCloneKiller.cpp diff --git a/Pr/PrPixel/CMakeLists.txt b/Pr/PrPixel/CMakeLists.txt index 1251c911b35..8102e5c70c7 100644 --- a/Pr/PrPixel/CMakeLists.txt +++ b/Pr/PrPixel/CMakeLists.txt @@ -21,6 +21,7 @@ gaudi_add_module(PrPixel src/VeloKalman.cpp src/PrVPHitsToVPLightClusters.cpp src/VeloHeavyFlavourTracking.cpp + src/VeloCloneKiller.cpp LINK Gaudi::GaudiAlgLib LHCb::DAQEventLib diff --git a/Pr/PrPixel/src/VeloCloneKiller.cpp b/Pr/PrPixel/src/VeloCloneKiller.cpp new file mode 100644 index 00000000000..9bd1b2240b8 --- /dev/null +++ b/Pr/PrPixel/src/VeloCloneKiller.cpp @@ -0,0 +1,126 @@ +/*****************************************************************************\ +* (c) Copyright 2000-2025 CERN for the benefit of the LHCb Collaboration * +* * +* This software is distributed under the terms of the GNU General Public * +* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". * +* * +* In applying this licence, CERN does not waive the privileges and immunities * +* granted to it by virtue of its status as an Intergovernmental Organization * +* or submit itself to any jurisdiction. * +\*****************************************************************************/ + +// Gaudi +#include "LHCbAlgs/Transformer.h" + +// LHCb +#include "Event/PrHits.h" +#include "Event/PrVeloTracks.h" +#include "Kernel/LHCbID.h" +#include "Kernel/VPConstants.h" +#include "LHCbMath/MatVec.h" +#include "VPDet/VPDetPaths.h" + +namespace { + template + F segment_distance( const LHCb::LinAlg::Vec& p, const LHCb::LinAlg::Vec& p1, + const LHCb::LinAlg::Vec& p2 ) { + auto u = p2 - p1; + auto AP = p - p1; + return AP.cross( u ).mag() / u.mag(); + } +} // namespace + +class VeloCloneKiller : public LHCb::Algorithm::MultiTransformer( + const EventContext&, const LHCb::Pr::Velo::Tracks&, const LHCb::Pr::VP::Hits& )> { +public: + VeloCloneKiller( const std::string& name, ISvcLocator* pSvcLocator ) + : MultiTransformer( name, pSvcLocator, { KeyValue{ "InputTracksLocation", "" }, KeyValue{ "HitsLocation", "" } }, + { KeyValue{ "OutputTracksLocation", "" } } ) {} + + std::tuple operator()( const EventContext& evtCtx, const LHCb::Pr::Velo::Tracks& input_tracks, + const LHCb::Pr::VP::Hits& velo_hits ) const override { + std::tuple result{ LHCb::Pr::Velo::Tracks{ + input_tracks.backward(), Zipping::generateZipIdentifier(), LHCb::getMemResource( evtCtx ) } }; + auto& [out] = result; + out.reserve( input_tracks.size() ); + + using TracksTag = LHCb::Pr::Velo::Tag; + + std::vector is_clone; + is_clone.resize( input_tracks.size() ); + + auto hits = velo_hits.scalar(); + auto tracks = input_tracks.scalar(); + for ( unsigned i = 0; i < tracks.size(); i++ ) { + if ( is_clone[i] ) continue; + + auto t1 = tracks[i]; + auto t1_hits = t1.template field(); + auto t1p0 = hits[t1_hits[0].get().cast()].template get(); + auto t1p1 = hits[t1_hits[t1_hits.size().cast() - 1].get().cast()] + .template get(); + + bool clone = false; + int clone_idx = -1; + for ( unsigned j = i + 1; j < tracks.size(); j++ ) { + if ( is_clone[j] ) continue; + auto t2 = tracks[j]; + auto t2_hits = t2.template field(); + auto t2p0 = hits[t2_hits[0].get().cast()].template get(); + auto t2p1 = hits[t2_hits[t2_hits.size().cast() - 1].get().cast()] + .template get(); + + auto d0 = segment_distance( t1p0, t2p0, t2p1 ).cast(); + auto d1 = segment_distance( t1p1, t2p0, t2p1 ).cast(); + auto d2 = segment_distance( t2p0, t1p0, t1p1 ).cast(); + auto d3 = segment_distance( t2p1, t1p0, t1p1 ).cast(); + + clone = ( d0 < m_threshold ) && ( d1 < m_threshold ) && ( d2 < m_threshold ) && ( d3 < m_threshold ); + + if ( clone ) { + is_clone[j] = true; + clone_idx = j; + break; + } + } + + int copy_idx = out.size(); + out.copy_back( t1 ); + + if ( clone ) { // merge tracks (keep hits order, remove duplicates) + auto tout = out.scalar()[copy_idx].template field(); + auto t1 = tracks[i].template field(); + auto t2 = tracks[clone_idx].template field(); + tout.resize( 0 ); + int a = 0; + int b = 0; + while ( a < t1.size().cast() || b < t2.size().cast() ) { + auto hout = tout.emplace_back(); + auto ia = ( a < t1.size().cast() ) ? t1[a].get().cast() : 0x7FFFFFFF; + auto ib = ( b < t2.size().cast() ) ? t2[b].get().cast() : 0x7FFFFFFF; + if ( ia > ib ) { + hout.template field().set( ib ); + hout.template field().set( t2[b].get() ); + b++; + } else { + hout.template field().set( ia ); + hout.template field().set( t1[a].get() ); + a++; + b += ( ia == ib ); + } + } + } + } + + m_nbCloneRemoved += input_tracks.size() - out.size(); + + return result; + } + +private: + mutable Gaudi::Accumulators::SummingCounter<> m_nbCloneRemoved{ this, "Nb of removed clones" }; + Gaudi::Property m_threshold{ this, "Threshold", 0.07f }; // Max distance (in mm) of any point on one segment to + // the line defined by the other track +}; + +DECLARE_COMPONENT( VeloCloneKiller ) diff --git a/Tr/TrackMonitors/src/TrackVPOverlapMonitor.cpp b/Tr/TrackMonitors/src/TrackVPOverlapMonitor.cpp index 2becdff241d..dee233a9977 100644 --- a/Tr/TrackMonitors/src/TrackVPOverlapMonitor.cpp +++ b/Tr/TrackMonitors/src/TrackVPOverlapMonitor.cpp @@ -376,9 +376,12 @@ namespace LHCb::Tr::Monitor { ++m_overlap_residual_x_CSO_NSI[{ module, localresidual.x() }]; ++m_overlap_residual_y_CSO_NSI[{ module, localresidual.y() }]; ++module_CSO_NSI_histo[module]; - } else { - warning() << "How can these sensors overlap?: " << sensor1 << " " << sensor2 << endmsg; - } + } /*else { + // This warning can be triggered if we do clone killing on velo tracks + // and merge the clones: + warning() << "How can these sensors overlap?: " << sensor1 << " " << sensor2 + << " rx:" << localresidual.x() << " ry:" << localresidual.y() << endmsg; + }*/ if ( m_expertMode.value() ) { Tuple nodeTuple = nTuple( "TrackVPOverlapMonitor_overlapnodes", "" ); -- GitLab From 69f6ba4e2ce10a0be25d13a049e840dc4bc2180b Mon Sep 17 00:00:00 2001 From: Arthur Hennequin Date: Wed, 7 May 2025 16:27:42 +0200 Subject: [PATCH 2/4] Apply suggestion --- Tr/TrackMonitors/src/TrackVPOverlapMonitor.cpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/Tr/TrackMonitors/src/TrackVPOverlapMonitor.cpp b/Tr/TrackMonitors/src/TrackVPOverlapMonitor.cpp index dee233a9977..63d7f9014f2 100644 --- a/Tr/TrackMonitors/src/TrackVPOverlapMonitor.cpp +++ b/Tr/TrackMonitors/src/TrackVPOverlapMonitor.cpp @@ -302,7 +302,7 @@ namespace LHCb::Tr::Monitor { ++m_Aresidual_y_vs_station[{ nodeA.vpid.station(), nodeA.residualY }]; ++m_Cresidual_x_vs_station[{ nodeA.vpid.station(), nodeC.residualX }]; ++m_Cresidual_y_vs_station[{ nodeA.vpid.station(), nodeC.residualY }]; - } else { + } else if ( inode->vpid.sensor() != jnode->vpid.sensor() ) { // Collect nodes in different sensors ont he same module. These are always consecutive. // Ladder map: { L0: CLI, L1: NLO, L2, NSI, L3: CSO } // Expected overlaps: CLI-NLO, CSO-NSI, CLI-NSI; Not allowed: CSO-CLI, NLO-NSI and CSO-NLO @@ -376,12 +376,10 @@ namespace LHCb::Tr::Monitor { ++m_overlap_residual_x_CSO_NSI[{ module, localresidual.x() }]; ++m_overlap_residual_y_CSO_NSI[{ module, localresidual.y() }]; ++module_CSO_NSI_histo[module]; - } /*else { - // This warning can be triggered if we do clone killing on velo tracks - // and merge the clones: + } else { warning() << "How can these sensors overlap?: " << sensor1 << " " << sensor2 - << " rx:" << localresidual.x() << " ry:" << localresidual.y() << endmsg; - }*/ + << " rx:" << localresidual.x() << " ry:" << localresidual.y() << endmsg; + } if ( m_expertMode.value() ) { Tuple nodeTuple = nTuple( "TrackVPOverlapMonitor_overlapnodes", "" ); -- GitLab From 9f0ffced86716b42809632a9e7547061de3966ae Mon Sep 17 00:00:00 2001 From: Arthur Hennequin Date: Tue, 13 May 2025 08:41:21 +0200 Subject: [PATCH 3/4] One less sqrt in segment distance --- Pr/PrPixel/src/VeloCloneKiller.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Pr/PrPixel/src/VeloCloneKiller.cpp b/Pr/PrPixel/src/VeloCloneKiller.cpp index 9bd1b2240b8..6043fda593c 100644 --- a/Pr/PrPixel/src/VeloCloneKiller.cpp +++ b/Pr/PrPixel/src/VeloCloneKiller.cpp @@ -26,7 +26,7 @@ namespace { const LHCb::LinAlg::Vec& p2 ) { auto u = p2 - p1; auto AP = p - p1; - return AP.cross( u ).mag() / u.mag(); + return sqrt( AP.cross( u ).mag2() / u.mag2() ); } } // namespace -- GitLab From 5a5cfa8702437626c0c65a876296f8460de6b1d6 Mon Sep 17 00:00:00 2001 From: Arthur Hennequin Date: Fri, 11 Jul 2025 11:51:44 +0200 Subject: [PATCH 4/4] Sort track to speed up comparisons --- Pr/PrPixel/src/VeloCloneKiller.cpp | 56 ++++++++++++++++++++++-------- 1 file changed, 41 insertions(+), 15 deletions(-) diff --git a/Pr/PrPixel/src/VeloCloneKiller.cpp b/Pr/PrPixel/src/VeloCloneKiller.cpp index 6043fda593c..b7bf721e730 100644 --- a/Pr/PrPixel/src/VeloCloneKiller.cpp +++ b/Pr/PrPixel/src/VeloCloneKiller.cpp @@ -28,6 +28,16 @@ namespace { auto AP = p - p1; return sqrt( AP.cross( u ).mag2() / u.mag2() ); } + + template + struct TrackRef { + TrackRef( unsigned idx, LHCb::LinAlg::Vec p1, LHCb::LinAlg::Vec p2 ) + : idx( idx ), p1( p1 ), p2( p2 ), tx( ( p2.x() - p1.x() ) / ( p2.z() - p1.z() ) ) {} + unsigned idx; + LHCb::LinAlg::Vec p1; + LHCb::LinAlg::Vec p2; + F tx; + }; } // namespace class VeloCloneKiller : public LHCb::Algorithm::MultiTransformer( @@ -51,24 +61,39 @@ public: auto hits = velo_hits.scalar(); auto tracks = input_tracks.scalar(); + + std::vector> sortedTracks; + sortedTracks.reserve( input_tracks.size() ); + for ( unsigned i = 0; i < tracks.size(); i++ ) { - if ( is_clone[i] ) continue; + auto t = tracks[i]; + auto t_hits = t.template field(); + auto tp0 = hits[t_hits[0].get().cast()].template get(); + auto tp1 = hits[t_hits[t_hits.size().cast() - 1].get().cast()] + .template get(); - auto t1 = tracks[i]; - auto t1_hits = t1.template field(); - auto t1p0 = hits[t1_hits[0].get().cast()].template get(); - auto t1p1 = hits[t1_hits[t1_hits.size().cast() - 1].get().cast()] - .template get(); + sortedTracks.emplace_back( i, tp0, tp1 ); + } + + std::sort( std::begin( sortedTracks ), std::end( sortedTracks ), + []( const auto& a, const auto& b ) { return a.tx.cast() < b.tx.cast(); } ); + + for ( unsigned i = 0; i < sortedTracks.size(); i++ ) { + if ( is_clone[i] ) continue; + auto t1p0 = sortedTracks[i].p1; + auto t1p1 = sortedTracks[i].p2; + auto tx1 = sortedTracks[i].tx; bool clone = false; int clone_idx = -1; - for ( unsigned j = i + 1; j < tracks.size(); j++ ) { + for ( unsigned j = i + 1; j < sortedTracks.size(); j++ ) { if ( is_clone[j] ) continue; - auto t2 = tracks[j]; - auto t2_hits = t2.template field(); - auto t2p0 = hits[t2_hits[0].get().cast()].template get(); - auto t2p1 = hits[t2_hits[t2_hits.size().cast() - 1].get().cast()] - .template get(); + + auto tx2 = sortedTracks[j].tx; + if ( tx2.cast() - tx1.cast() > m_tx_threshold ) break; + + auto t2p0 = sortedTracks[j].p1; + auto t2p1 = sortedTracks[j].p2; auto d0 = segment_distance( t1p0, t2p0, t2p1 ).cast(); auto d1 = segment_distance( t1p1, t2p0, t2p1 ).cast(); @@ -85,12 +110,12 @@ public: } int copy_idx = out.size(); - out.copy_back( t1 ); + out.copy_back( tracks[sortedTracks[i].idx] ); if ( clone ) { // merge tracks (keep hits order, remove duplicates) auto tout = out.scalar()[copy_idx].template field(); - auto t1 = tracks[i].template field(); - auto t2 = tracks[clone_idx].template field(); + auto t1 = tracks[sortedTracks[i].idx].template field(); + auto t2 = tracks[sortedTracks[clone_idx].idx].template field(); tout.resize( 0 ); int a = 0; int b = 0; @@ -121,6 +146,7 @@ private: mutable Gaudi::Accumulators::SummingCounter<> m_nbCloneRemoved{ this, "Nb of removed clones" }; Gaudi::Property m_threshold{ this, "Threshold", 0.07f }; // Max distance (in mm) of any point on one segment to // the line defined by the other track + Gaudi::Property m_tx_threshold{ this, "TxThreshold", 0.001f }; // slope threshold for early exit }; DECLARE_COMPONENT( VeloCloneKiller ) -- GitLab