From 79bde340f94ee86815b03832390a868f6096d063 Mon Sep 17 00:00:00 2001 From: Gerhard Raven Date: Fri, 27 Aug 2021 12:18:00 +0200 Subject: [PATCH 1/3] Streamline FastJetBuiler (follow up to Phys!954) * add support for using a few enums as Gaudi::Property * prefer enums to be Gaudi::Property - this allow to have readable string in the configuration as as opposed to having magic integer values, and avoids the need for updateHandlers * follow introduction of LHCb::Algorithm as frontend to Gaudi::Functional * avoid 'OUT&' as input arguments to functions, prefer to return them by value instead --- Phys/JetAccessories/src/FastJetBuilder.cpp | 344 +++++++++++---------- 1 file changed, 181 insertions(+), 163 deletions(-) diff --git a/Phys/JetAccessories/src/FastJetBuilder.cpp b/Phys/JetAccessories/src/FastJetBuilder.cpp index c7285b2aa..331ca3a56 100644 --- a/Phys/JetAccessories/src/FastJetBuilder.cpp +++ b/Phys/JetAccessories/src/FastJetBuilder.cpp @@ -8,48 +8,26 @@ * granted to it by virtue of its status as an Intergovernmental Organization * * or submit itself to any jurisdiction. * \*****************************************************************************/ - -// STL. -#include -#include -#include - -// Gaudi -#include "GaudiAlg/GaudiTool.h" -#include "GaudiAlg/Transformer.h" - -// LHCb #include "DetDesc/DetectorElement.h" #include "DetDesc/GenericConditionAccessorHolder.h" - -// Jets. -#include "Kernel/IParticleCombiner.h" -#include "Kernel/JetEnums.h" - -// Event. #include "Event/Particle.h" #include "Event/RecVertex.h" #include "Event/Track.h" - -// LoKi. +#include "GaudiAlg/GaudiTool.h" +#include "Kernel/IParticleCombiner.h" +#include "Kernel/IRelatedPVFinder.h" +#include "Kernel/JetEnums.h" +#include "Kernel/Particle2Vertex.h" +#include "LHCbAlgs/Transformer.h" #include "LoKi/ParticleContextCuts.h" #include "LoKi/ParticleCuts.h" #include "LoKi/PhysExtract.h" #include "LoKi/VertexCuts.h" - -// Root. #include "Math/VectorUtil.h" - -// Track to vertex relations -#include "Kernel/IRelatedPVFinder.h" -#include "Kernel/Particle2Vertex.h" #include "Relations/Relation1D.h" - #include "fastjet/JetDefinition.hh" #include "fastjet/PseudoJet.hh" - -using namespace std; -using namespace LHCb; +#include /* Run3 Jet building (JB) algorithm class for use in the HLT and offline. @@ -77,12 +55,29 @@ in root files located in ParamFiles package. TODO: Add a JetID filter , probably as a separate transformer. */ -using P2PVTable = unique_ptr>; +using P2PVTable = std::unique_ptr>; // the actual type of input data container using InputJ = LHCb::Particle::ConstVector; // FastJet Jet type using PJet = fastjet::PseudoJet; +namespace FastJetDetails { + enum struct Ordering { by_rapidity = 3, by_pt = 2, by_E = 1 }; + + template + auto sorted_by( Ordering order, In&& input ) { + switch ( order ) { + case Ordering::by_rapidity: + return fastjet::sorted_by_rapidity( std::forward( input ) ); + case Ordering::by_pt: + return fastjet::sorted_by_pt( std::forward( input ) ); + case Ordering::by_E: + return fastjet::sorted_by_E( std::forward( input ) ); + } + throw std::out_of_range( "unknown enumeration value" ); + } +} // namespace FastJetDetails + namespace { int to_user_index( const int index ) { return index + 10000; } int from_user_index( const int index ) { return index - 10000; } @@ -93,12 +88,101 @@ namespace { jetinput.set_user_index( index ); return jetinput; } + + using namespace std::string_view_literals; + + constexpr auto StrategyMap = std::array{std::pair{fastjet::Strategy::N2MinHeapTiled, "N2MinHeapTiled"sv}, + std::pair{fastjet::Strategy::N2Tiled, "N2Tiled"sv}, + std::pair{fastjet::Strategy::N2PoorTiled, "N2PoorTiled"sv}, + std::pair{fastjet::Strategy::N2Plain, "N2Plain"sv}, + std::pair{fastjet::Strategy::N3Dumb, "N3Dumb"sv}, + std::pair{fastjet::Strategy::Best, "Best"sv}, + std::pair{fastjet::Strategy::NlnN, "NlnN"sv}, + std::pair{fastjet::Strategy::NlnN3pi, "NlnN3pi"sv}, + std::pair{fastjet::Strategy::NlnN4pi, "NlnN4pi"sv}, + std::pair{fastjet::Strategy::NlnNCam4pi, "NlnNCam4pi"sv}, + std::pair{fastjet::Strategy::NlnNCam2pi2R, "NlnNCam2pi2R"sv}, + std::pair{fastjet::Strategy::NlnNCam, "NlnNCam"sv}, + std::pair{fastjet::Strategy::plugin_strategy, "plugin_strategy"sv}}; + + constexpr auto JetAlgorithmMap = + std::array{std::pair{fastjet::JetAlgorithm::kt_algorithm, "kt_algorithm"sv}, + std::pair{fastjet::JetAlgorithm::cambridge_algorithm, "cambridge_algorithm"sv}, + std::pair{fastjet::JetAlgorithm::antikt_algorithm, "antikt_algorithm"sv}, + std::pair{fastjet::JetAlgorithm::genkt_algorithm, "genkt_algorithm"sv}, + std::pair{fastjet::JetAlgorithm::cambridge_for_passive_algorithm, "cambridge_for_passive_algorithm"sv}, + std::pair{fastjet::JetAlgorithm::genkt_for_passive_algorithm, "genkt_for_passive_algorithm"sv}, + std::pair{fastjet::JetAlgorithm::ee_kt_algorithm, "ee_kt_algorithm"sv}, + std::pair{fastjet::JetAlgorithm::ee_genkt_algorithm, "ee_genkt_algorithm"sv}, + std::pair{fastjet::JetAlgorithm::plugin_algorithm, "plugin_algorithm"sv}}; + + /// the various recombination schemes + constexpr auto RecombinationSchemeMap = + std::array{std::pair{fastjet::RecombinationScheme::E_scheme, "E_scheme"sv}, + std::pair{fastjet::RecombinationScheme::pt_scheme, "pt_scheme"sv}, + std::pair{fastjet::RecombinationScheme::pt2_scheme, "pt2_scheme"sv}, + std::pair{fastjet::RecombinationScheme::Et_scheme, "Et_scheme"sv}, + std::pair{fastjet::RecombinationScheme::Et2_scheme, "Et2_scheme"sv}, + std::pair{fastjet::RecombinationScheme::BIpt_scheme, "BIpt_scheme"sv}, + std::pair{fastjet::RecombinationScheme::BIpt2_scheme, "BIpt2_scheme"sv}, + std::pair{fastjet::RecombinationScheme::external_scheme, "external_scheme"sv}}; + + constexpr auto OrderingMap = std::array{ + std::pair{FastJetDetails::Ordering::by_rapidity, "by_rapidity"sv}, + std::pair{FastJetDetails::Ordering::by_pt, "by_pt"sv}, + std::pair{FastJetDetails::Ordering::by_E, "by_E"sv}, + }; + + template + std::string toString_( Enum e, std::array, N> const& tbl ) { + auto i = std::find_if( tbl.begin(), tbl.end(), [e]( const auto& entry ) { return entry.first == e; } ); + if ( i == tbl.end() ) throw std::out_of_range( "unknown enumeration value" ); + return std::string{i->second}; + } + + template + StatusCode parse_( Enum& e, std::string_view s, std::array, N> const& tbl ) { + if ( !s.empty() && s.front() == s.back() && ( s.front() == '\'' || s.front() == '\"' ) ) { + s.remove_prefix( 1 ); + s.remove_suffix( 1 ); + } + auto i = std::find_if( tbl.begin(), tbl.end(), [s]( const auto& entry ) { return entry.second == s; } ); + if ( i == tbl.end() ) return StatusCode::FAILURE; + e = i->first; + return StatusCode::SUCCESS; + } } // namespace +namespace FastJetDetails { + std::string toString( Ordering e ) { return toString_( e, OrderingMap ); } + StatusCode parse( Ordering& bt, const std::string& in ) { return parse_( bt, in, OrderingMap ); } + std::ostream& toStream( Ordering e, std::ostream& os ) { return os << std::quoted( toString( e ), '\'' ); } + std::ostream& operator<<( std::ostream& s, Ordering e ) { return toStream( e, s ); } +} // namespace FastJetDetails + +namespace fastjet { + std::string toString( RecombinationScheme e ) { return toString_( e, RecombinationSchemeMap ); } + StatusCode parse( RecombinationScheme& bt, const std::string& in ) { + return parse_( bt, in, RecombinationSchemeMap ); + } + std::ostream& toStream( RecombinationScheme e, std::ostream& os ) { return os << std::quoted( toString( e ), '\'' ); } + std::ostream& operator<<( std::ostream& s, RecombinationScheme e ) { return toStream( e, s ); } + + std::string toString( JetAlgorithm e ) { return toString_( e, JetAlgorithmMap ); } + StatusCode parse( JetFinder& bt, const std::string& in ) { return parse_( bt, in, JetAlgorithmMap ); } + std::ostream& toStream( JetAlgorithm e, std::ostream& os ) { return os << std::quoted( toString( e ), '\'' ); } + std::ostream& operator<<( std::ostream& s, JetFinder e ) { return toStream( e, s ); } + + std::string toString( Strategy e ) { return toString_( e, StrategyMap ); } + StatusCode parse( Strategy& bt, const std::string& in ) { return parse_( bt, in, StrategyMap ); } + std::ostream& toStream( Strategy e, std::ostream& os ) { return os << std::quoted( toString( e ), '\'' ); } + std::ostream& operator<<( std::ostream& s, Strategy e ) { return toStream( e, s ); } +} // namespace fastjet + class FastJetBuilder - : public Gaudi::Functional::Transformer> { + : public LHCb::Algorithm::Transformer> { public: FastJetBuilder( const std::string& name, ISvcLocator* svc ); // Constructor @@ -106,91 +190,44 @@ public: const DetectorElement& ) const override; // Main method private: - mutable Gaudi::Accumulators::SummingCounter<> m_nbInput{this, "Nb of Inputs"}; - mutable Gaudi::Accumulators::SummingCounter<> m_nbJetCounter{this, "Nb of Reconstructed Jets"}; - mutable Gaudi::Accumulators::SummingCounter<> m_nbJetCounter10{this, "Nb of Reconstructed Jets with pT>=10GeV"}; - mutable Gaudi::Accumulators::SummingCounter<> m_nbJetCounter20{this, "Nb of Reconstructed Jets with pT>=20GeV"}; - - ///< Add information to jet - void addInfo( Particle* prt ) const; - - Gaudi::Property m_jetVtx = {this, "JetsByVtx", true, - "If true, build jets for each primary vertex, otherwise " - "build inclusively."}; - - ///< Relate a particles with bestPVs - P2PVTable relateP2Vtx( std::vector Particles, const LHCb::RecVertices& PrimaryVertices ) const; - /// Tool to RelatedPVFinder ToolHandle m_prtVtxTool = { this, "RelatedPVFinder", "GenericParticle2PVRelator__p2PVWithIPChi2_OfflineDistanceCalculatorName_/P2PVWithIPChi2"}; - // Input/output property members (Inputs and Output handled by base class). - std::vector m_inBans; ///< Input locations of Particles to ban. - std::vector m_inTags; ///< Input locations of Particle to tag. + // Combiner to be used + ToolHandle m_combiner{this, "ParticleCombiner", "MomentumCombiner"}; - // Below: FastJet related + Gaudi::Property m_jetVtx = {this, "JetsByVtx", true, + "If true, build jets for each primary vertex, otherwise " + "build inclusively."}; // FastJet configuration - enum m_sort { m_sorted_by_rapidity = 3, m_sorted_by_pt = 2, m_sorted_by_E = 1 }; - Gaudi::Property m_sort = {this, "Sort", 2, "FastJet: sorting criteria"}; - Gaudi::Property m_jetID = {this, "JetID", 98, "LHCb PID number for jets"}; - Gaudi::Property m_r = {this, "RParameter", 0.5, "FastJet: cone size"}; - Gaudi::Property m_ptmin = {this, "PtMin", 5000, + Gaudi::Property m_sort = {this, "Sort", FastJetDetails::Ordering::by_pt, + "FastJet: sorting criteria"}; + Gaudi::Property m_jetID = {this, "JetID", 98, "LHCb PID number for jets"}; + Gaudi::Property m_r = {this, "RParameter", 0.5, "FastJet: cone size"}; + Gaudi::Property m_ptmin = {this, "PtMin", 5000, [=]( auto& ) { if ( m_ptmin.value() < 0 ) throw std::runtime_error{"min pt for FastJet is negative for FastJetBuilder"}; }, "FastJet: min pT"}; - Gaudi::Property m_recom = {this, "Recombination", 0, - [=]( auto& ) { - constexpr auto valid = - std::array{fastjet::E_scheme, fastjet::pt_scheme, fastjet::pt2_scheme, - fastjet::Et_scheme, fastjet::Et2_scheme, fastjet::BIpt_scheme, - fastjet::BIpt2_scheme}; - if ( std::find( valid.begin(), valid.end(), m_strat.value() ) == valid.end() ) - throw std::runtime_error{"Invalid RecombinationScheme for FastJetBuilder"}; - }, - "FastJet: scheme : 0 E, 1 pT, 2 pT^2, ..."}; - - Gaudi::Property m_strat = {this, "Strategy", 1, - [=]( auto& ) { - constexpr auto valid = - std::array{fastjet::N2MinHeapTiled, fastjet::N2Tiled, fastjet::N2PoorTiled, - fastjet::N2Plain, fastjet::N3Dumb, fastjet::Best, - fastjet::NlnN, fastjet::NlnN3pi, fastjet::NlnN4pi, - fastjet::NlnNCam4pi, fastjet::NlnNCam2pi2R, fastjet::NlnNCam, - fastjet::plugin_strategy}; - if ( std::find( valid.begin(), valid.end(), m_strat.value() ) == valid.end() ) - throw std::runtime_error{"Invalid Strategy for FastJetBuilder"}; - }, - "FastJet: strategy: 0 N3Dumb, 1 Best, 2 NlnN, ...,"}; - Gaudi::Property m_type = {this, "Type", 2, - [=]( auto& ) { - constexpr auto valid = std::array{fastjet::kt_algorithm, - fastjet::cambridge_algorithm, - fastjet::genkt_algorithm, - fastjet::cambridge_for_passive_algorithm, - fastjet::genkt_for_passive_algorithm, - fastjet::plugin_algorithm}; - if ( std::find( valid.begin(), valid.end(), m_strat.value() ) == valid.end() ) - throw std::runtime_error{"Invalid JetFinder-algorithm for FastJetBuilder"}; - }, - "FastJet: anti-kt code: 0 kt, 1 cambridge, 2 anti-kt, ..."}; - - // Combiner to be used - ToolHandle m_combiner{this, "ParticleCombiner", "MomentumCombiner"}; + Gaudi::Property m_recom = {this, "Recombination", + fastjet::RecombinationScheme::E_scheme, "FastJet: scheme "}; + Gaudi::Property m_strat = {this, "Strategy", fastjet::Strategy::Best, "FastJet: strategy"}; + Gaudi::Property m_type = {this, "Type", fastjet::JetAlgorithm::antikt_algorithm, + "FastJet: JetAlgorithm"}; - // run fastjet - StatusCode makeFastJets( const std::vector Particles, LHCb::Particle::Vector& recojets, - const DetectorElement& geometry ) const; + mutable Gaudi::Accumulators::SummingCounter<> m_nbInput{this, "Nb of Inputs"}; + mutable Gaudi::Accumulators::SummingCounter<> m_nbJetCounter{this, "Nb of Reconstructed Jets"}; + mutable Gaudi::Accumulators::SummingCounter<> m_nbJetCounter10{this, "Nb of Reconstructed Jets with pT>=10GeV"}; + mutable Gaudi::Accumulators::SummingCounter<> m_nbJetCounter20{this, "Nb of Reconstructed Jets with pT>=20GeV"}; /// MsgSvc counters mutable Gaudi::Accumulators::MsgCounter m_msg_noBestPV{this, "Failed to relate a particle with its bestPV."}; - mutable Gaudi::Accumulators::MsgCounter m_msg_FJfailed{this, "FastJet failed."}; mutable Gaudi::Accumulators::MsgCounter m_msg_invalid_input{this, "Invalid input particle -> FastJet input"}; mutable Gaudi::Accumulators::MsgCounter m_msg_nojetdaughters{this, "Jets have no daughters"}; @@ -198,9 +235,18 @@ private: mutable Gaudi::Accumulators::MsgCounter m_msg_momcomb{this, "Error in momentum combiner with jet daughters"}; + ///< Relate a particles with bestPVs + P2PVTable relateP2Vtx( std::vector Particles, const LHCb::RecVertices& PrimaryVertices ) const; + + ///< Add information to jet + void addInfo( LHCb::Particle& prt ) const; + + // run fastjet + LHCb::Particle::Vector makeFastJets( std::vector const& Particles, + const DetectorElement& geometry ) const; + // prepare fastjet definition - std::pair> - prepare( const std::vector Particles ) const; + std::vector prepare( const std::vector Particles ) const; }; /*----------------------------------------------------------------------------- @@ -264,9 +310,7 @@ LHCb::Particles FastJetBuilder::operator()( const LHCb::Particles& Particles, co for ( const auto prt : tmpneutrals ) tmppart.push_back( prt ); // Make the jets - std::vector tmprecojets; - tmprecojets.reserve( 40 ); - if ( makeFastJets( tmppart, tmprecojets, geometry ).isFailure() ) ++m_msg_FJfailed; + auto tmprecojets = makeFastJets( tmppart, geometry ); if ( msgLevel( MSG::DEBUG ) ) debug() << "Vtx: " << vtx << " #jets: " << tmprecojets.size() << endmsg; recojets.reserve( tmprecojets.size() ); @@ -285,9 +329,7 @@ LHCb::Particles FastJetBuilder::operator()( const LHCb::Particles& Particles, co std::vector tmppart{Particles.begin(), Particles.end()}; // Make the jets - std::vector tmprecojets; - tmprecojets.reserve( 40 ); - if ( makeFastJets( tmppart, tmprecojets, geometry ).isFailure() ) ++m_msg_FJfailed; + auto tmprecojets = makeFastJets( tmppart, geometry ); recojets.reserve( tmprecojets.size() ); for ( auto* jet : tmprecojets ) { @@ -300,7 +342,7 @@ LHCb::Particles FastJetBuilder::operator()( const LHCb::Particles& Particles, co } for ( auto* jet : recojets ) { - addInfo( jet ); + addInfo( *jet ); const auto pt = jet->momentum().Pt(); if ( pt >= 10000 ) { c10 += 1; } if ( pt >= 20000 ) { c20 += 1; } @@ -325,10 +367,10 @@ LHCb::Particles FastJetBuilder::operator()( const LHCb::Particles& Particles, co } // Add additional information. -void FastJetBuilder::addInfo( Particle* prt ) const { - double cpx( 0 ), cpy( 0 ), cptMax( 0 ), nptMax( 0 ), width( 0 ), norm( 0 ), trks( 0 ), pt( prt->momentum().Pt() ); - Particle::ConstVector dtrs; - LoKi::Extract::getParticles( prt, back_inserter( dtrs ), LoKi::Cuts::HASPROTO && LoKi::Cuts::ISBASIC ); +void FastJetBuilder::addInfo( LHCb::Particle& prt ) const { + double cpx( 0 ), cpy( 0 ), cptMax( 0 ), nptMax( 0 ), width( 0 ), norm( 0 ), trks( 0 ), pt( prt.momentum().Pt() ); + LHCb::Particle::ConstVector dtrs; + LoKi::Extract::getParticles( &prt, back_inserter( dtrs ), LoKi::Cuts::HASPROTO && LoKi::Cuts::ISBASIC ); for ( auto dtr : dtrs ) { const Gaudi::LorentzVector& vec = dtr->momentum(); if ( dtr->charge() != 0 ) { @@ -338,17 +380,17 @@ void FastJetBuilder::addInfo( Particle* prt ) const { ++trks; } else if ( vec.Pt() > nptMax ) nptMax = vec.Pt(); - width += ROOT::Math::VectorUtil::DeltaR( vec, prt->momentum() ) * vec.Pt(); + width += ROOT::Math::VectorUtil::DeltaR( vec, prt.momentum() ) * vec.Pt(); norm += vec.Pt(); } - prt->addInfo( LHCb::JetIDInfo::Ntracks, trks ); - prt->addInfo( LHCb::JetIDInfo::MTF, cptMax / pt ); - prt->addInfo( LHCb::JetIDInfo::MNF, nptMax / pt ); - prt->addInfo( LHCb::JetIDInfo::MPT, cptMax ); - prt->addInfo( LHCb::JetIDInfo::CPF, sqrt( cpx * cpx + cpy * cpy ) / pt ); - prt->addInfo( LHCb::JetIDInfo::JetWidth, width ); - prt->addInfo( LHCb::JetIDInfo::JetWidthNorm, width / norm ); // Added to Phys/Phys/JetAccessories/Kernel/JetEnums.h. - // Remove this or the one above after tests + prt.addInfo( LHCb::JetIDInfo::Ntracks, trks ); + prt.addInfo( LHCb::JetIDInfo::MTF, cptMax / pt ); + prt.addInfo( LHCb::JetIDInfo::MNF, nptMax / pt ); + prt.addInfo( LHCb::JetIDInfo::MPT, cptMax ); + prt.addInfo( LHCb::JetIDInfo::CPF, sqrt( cpx * cpx + cpy * cpy ) / pt ); + prt.addInfo( LHCb::JetIDInfo::JetWidth, width ); + prt.addInfo( LHCb::JetIDInfo::JetWidthNorm, width / norm ); // Added to Phys/Phys/JetAccessories/Kernel/JetEnums.h. + // Remove this or the one above after tests } // Relate particles with bestPVs @@ -378,21 +420,8 @@ P2PVTable FastJetBuilder::relateP2Vtx( std::vector Partic } // prepare FastJet -std::pair> -FastJetBuilder::prepare( const std::vector Particles ) const { - - fastjet::JetFinder finder = (fastjet::JetFinder)m_type.value(); - fastjet::RecombinationScheme scheme = (fastjet::RecombinationScheme)m_recom.value(); - fastjet::Strategy strategy = (fastjet::Strategy)m_strat.value(); - - auto ret = std::pair{fastjet::JetDefinition{finder, m_r, scheme, strategy}, std::vector{}}; - auto& [jet_def, jets] = ret; - - if ( msgLevel( MSG::DEBUG ) ) { - debug() << "fastjet::JetDefinition:" << endmsg; - debug() << jet_def.description() << endmsg; - } - +std::vector FastJetBuilder::prepare( const std::vector Particles ) const { + std::vector jets; jets.reserve( Particles.size() ); for ( const auto& [i, p] : LHCb::range::enumerate( Particles ) ) { if ( !p ) { @@ -401,53 +430,43 @@ FastJetBuilder::prepare( const std::vector Particles ) co } jets.push_back( createJetInput( *p, to_user_index( i ) ) ); } - - return ret; + return jets; } + // Run FastJet -StatusCode FastJetBuilder::makeFastJets( const std::vector Particles, - LHCb::Particle::Vector& recojets, const DetectorElement& geometry ) const { +LHCb::Particle::Vector FastJetBuilder::makeFastJets( std::vector const& Particles, + const DetectorElement& geometry ) const { // prepare the input data and define the jets - auto [jetDef, jetinputs] = prepare( Particles ); + auto jetinputs = prepare( Particles ); if ( jetinputs.empty() ) { - recojets.reserve( 0 ); if ( msgLevel( MSG::DEBUG ) ) { debug() << "No particle inputs to reconstruct jets" << endmsg; } - return StatusCode::SUCCESS; + return {}; } // Reconstructed Jets from FastJet - std::vector jets; - auto clusters = fastjet::ClusterSequence( jetinputs, jetDef ); - - switch ( m_sort ) { - case FastJetBuilder::m_sorted_by_rapidity: - jets = fastjet::sorted_by_rapidity( clusters.inclusive_jets( m_ptmin ) ); - break; - case FastJetBuilder::m_sorted_by_pt: - jets = fastjet::sorted_by_pt( clusters.inclusive_jets( m_ptmin ) ); - break; - case FastJetBuilder::m_sorted_by_E: - jets = fastjet::sorted_by_E( clusters.inclusive_jets( m_ptmin ) ); - break; - default: - jets = fastjet::sorted_by_pt( clusters.inclusive_jets( m_ptmin ) ); - break; + auto jetDef = fastjet::JetDefinition{m_type.value(), m_r, m_recom.value(), m_strat.value()}; + if ( msgLevel( MSG::DEBUG ) ) { + debug() << "fastjet::JetDefinition:" << endmsg; + debug() << jetDef.description() << endmsg; } + auto clusters = fastjet::ClusterSequence( jetinputs, jetDef ); + auto jets = sorted_by( m_sort.value(), clusters.inclusive_jets( m_ptmin ) ); + if ( msgLevel( MSG::DEBUG ) ) { debug() << "No jets from fastjet::ClusterSequence" << endmsg; } + LHCb::Particle::Vector recojets; recojets.reserve( jets.size() ); for ( const PJet& jet : jets ) { - const std::vector& constituents = clusters.constituents( jet ); // jet daughter list to be used in the combiner LHCb::Particle::ConstVector daughters; - for ( const PJet& c : constituents ) { + for ( const PJet& c : clusters.constituents( jet ) ) { // find the appropriate input particle const int index = from_user_index( c.user_index() ); if ( 0 > index || (int)jetinputs.size() <= index ) { @@ -473,6 +492,5 @@ StatusCode FastJetBuilder::makeFastJets( const std::vector Date: Sat, 28 Aug 2021 00:20:37 +0200 Subject: [PATCH 2/3] More streamlining * avoid unique_ptr for P2PVTable * explicitly used unique_ptr when lifetime matters * avoid copying vectors when not strictly needed --- Phys/JetAccessories/src/FastJetBuilder.cpp | 140 ++++++++++----------- 1 file changed, 65 insertions(+), 75 deletions(-) diff --git a/Phys/JetAccessories/src/FastJetBuilder.cpp b/Phys/JetAccessories/src/FastJetBuilder.cpp index 331ca3a56..4389c989b 100644 --- a/Phys/JetAccessories/src/FastJetBuilder.cpp +++ b/Phys/JetAccessories/src/FastJetBuilder.cpp @@ -12,7 +12,6 @@ #include "DetDesc/GenericConditionAccessorHolder.h" #include "Event/Particle.h" #include "Event/RecVertex.h" -#include "Event/Track.h" #include "GaudiAlg/GaudiTool.h" #include "Kernel/IParticleCombiner.h" #include "Kernel/IRelatedPVFinder.h" @@ -27,7 +26,6 @@ #include "Relations/Relation1D.h" #include "fastjet/JetDefinition.hh" #include "fastjet/PseudoJet.hh" -#include /* Run3 Jet building (JB) algorithm class for use in the HLT and offline. @@ -55,7 +53,7 @@ in root files located in ParamFiles package. TODO: Add a JetID filter , probably as a separate transformer. */ -using P2PVTable = std::unique_ptr>; +using P2PVTable = LHCb::Relation1D; // the actual type of input data container using InputJ = LHCb::Particle::ConstVector; // FastJet Jet type @@ -79,12 +77,12 @@ namespace FastJetDetails { } // namespace FastJetDetails namespace { - int to_user_index( const int index ) { return index + 10000; } - int from_user_index( const int index ) { return index - 10000; } - PJet createJetInput( const LHCb::Particle& p, const int index ) { + constexpr int to_user_index( int index ) { return index + 10000; } + constexpr int from_user_index( int index ) { return index - 10000; } + PJet createJetInput( const LHCb::Particle& p, const int index ) { // trivial function which "converts" LHCb particle into the fastjet particle - const Gaudi::LorentzVector& v = p.momentum(); - auto jetinput = PJet{v.Px(), v.Py(), v.Pz(), v.E()}; + const auto& v = p.momentum(); + auto jetinput = PJet{v.Px(), v.Py(), v.Pz(), v.E()}; jetinput.set_user_index( index ); return jetinput; } @@ -236,17 +234,17 @@ private: "Error in momentum combiner with jet daughters"}; ///< Relate a particles with bestPVs - P2PVTable relateP2Vtx( std::vector Particles, const LHCb::RecVertices& PrimaryVertices ) const; + P2PVTable relateP2Vtx( LHCb::span, const LHCb::RecVertices& ) const; ///< Add information to jet void addInfo( LHCb::Particle& prt ) const; // run fastjet - LHCb::Particle::Vector makeFastJets( std::vector const& Particles, - const DetectorElement& geometry ) const; + std::vector> makeFastJets( LHCb::span particles, + const DetectorElement& geometry ) const; // prepare fastjet definition - std::vector prepare( const std::vector Particles ) const; + std::vector prepare( LHCb::span Particles ) const; }; /*----------------------------------------------------------------------------- @@ -272,11 +270,6 @@ LHCb::Particles FastJetBuilder::operator()( const LHCb::Particles& Particles, co const DetectorElement& geometry ) const { m_nbInput += Particles.size(); - auto c10 = m_nbJetCounter10.buffer(); - auto c20 = m_nbJetCounter20.buffer(); - - // Output jets - LHCb::Particles recojets; if ( msgLevel( MSG::DEBUG ) ) { debug() << "New event: #Particles: " << Particles.size() << " PVs: " << PrimaryVertices.size() << " PV list ----> "; @@ -289,58 +282,60 @@ LHCb::Particles FastJetBuilder::operator()( const LHCb::Particles& Particles, co std::vector tmpneutrals; tmpparttracks.reserve( 0.7 * Particles.size() ); tmpneutrals.reserve( 0.5 * Particles.size() ); - int pid; for ( auto const* prt : Particles ) { - pid = prt->particleID().pid(); + int pid = prt->particleID().pid(); ( pid == 22 || pid == 111 ) ? tmpneutrals.push_back( prt ) : tmpparttracks.push_back( prt ); } - // Relate particles with tracks (or their composites) to best PV - auto vtx2prtRel = relateP2Vtx( tmpparttracks, PrimaryVertices ); + // Output jets + LHCb::Particles recojets; if ( m_jetVtx.value() ) { // Make jets by PVs - if ( vtx2prtRel ) { - if ( msgLevel( MSG::DEBUG ) ) debug() << "Making jets by PVs" << endmsg; - for ( auto vtx : PrimaryVertices ) { - auto tmppartRange = vtx2prtRel->relations( vtx ); - std::vector tmppart; - tmppart.reserve( tmppartRange.size() + tmpneutrals.size() ); - for ( const auto prt : tmppartRange ) tmppart.push_back( prt.to() ); - // Extend including neutral (They go to every PV) - for ( const auto prt : tmpneutrals ) tmppart.push_back( prt ); - - // Make the jets - auto tmprecojets = makeFastJets( tmppart, geometry ); - - if ( msgLevel( MSG::DEBUG ) ) debug() << "Vtx: " << vtx << " #jets: " << tmprecojets.size() << endmsg; - recojets.reserve( tmprecojets.size() ); - for ( auto* jet : tmprecojets ) { - recojets.insert( jet ); - jet->addInfo( LHCb::JetIDInfo::vtx_x, vtx->position().X() ); - jet->addInfo( LHCb::JetIDInfo::vtx_y, vtx->position().Y() ); - jet->addInfo( LHCb::JetIDInfo::vtx_z, vtx->position().Z() ); - } - } - } else + // Relate particles with tracks (or their composites) to best PV + if ( PrimaryVertices.empty() ) { + if ( msgLevel( MSG::DEBUG ) ) debug() << "No PrimaryVertices in this event" << endmsg; return recojets; + } + auto vtx2prtRel = relateP2Vtx( tmpparttracks, PrimaryVertices ); + if ( msgLevel( MSG::DEBUG ) ) debug() << "Making jets by PVs" << endmsg; + for ( auto const& vtx : PrimaryVertices ) { + auto tmppartRange = vtx2prtRel.relations( vtx ); + std::vector tmppart; + tmppart.reserve( tmppartRange.size() + tmpneutrals.size() ); + for ( const auto& prt : tmppartRange ) tmppart.push_back( prt.to() ); + // Extend including neutral (They go to every PV) + tmppart.insert( tmppart.end(), tmpneutrals.begin(), tmpneutrals.end() ); + + // Make the jets + auto tmprecojets = makeFastJets( tmppart, geometry ); + + if ( msgLevel( MSG::DEBUG ) ) debug() << "Vtx: " << vtx << " #jets: " << tmprecojets.size() << endmsg; + recojets.reserve( tmprecojets.size() ); + for ( auto& jet : tmprecojets ) { + jet->addInfo( LHCb::JetIDInfo::vtx_x, vtx->position().X() ); + jet->addInfo( LHCb::JetIDInfo::vtx_y, vtx->position().Y() ); + jet->addInfo( LHCb::JetIDInfo::vtx_z, vtx->position().Z() ); + recojets.insert( jet.release() ); + } + } } else { // Make jets inclusively if ( msgLevel( MSG::DEBUG ) ) debug() << "Making jets inclusively" << endmsg; - std::vector tmppart{Particles.begin(), Particles.end()}; - // Make the jets - auto tmprecojets = makeFastJets( tmppart, geometry ); - + auto tmprecojets = makeFastJets( std::vector{Particles.begin(), Particles.end()}, geometry ); recojets.reserve( tmprecojets.size() ); - for ( auto* jet : tmprecojets ) { - auto vtx = m_prtVtxTool->relatedPV( jet, PrimaryVertices ); + for ( auto& jet : tmprecojets ) { + auto vtx = m_prtVtxTool->relatedPV( jet.get(), PrimaryVertices ); jet->addInfo( LHCb::JetIDInfo::vtx_x, vtx->position().X() ); jet->addInfo( LHCb::JetIDInfo::vtx_y, vtx->position().Y() ); jet->addInfo( LHCb::JetIDInfo::vtx_z, vtx->position().Z() ); - recojets.insert( jet ); + recojets.insert( jet.release() ); } } + m_nbJetCounter += recojets.size(); + auto c10 = m_nbJetCounter10.buffer(); + auto c20 = m_nbJetCounter20.buffer(); for ( auto* jet : recojets ) { addInfo( *jet ); const auto pt = jet->momentum().Pt(); @@ -362,7 +357,6 @@ LHCb::Particles FastJetBuilder::operator()( const LHCb::Particles& Particles, co } } - m_nbJetCounter += recojets.size(); return recojets; } @@ -371,7 +365,7 @@ void FastJetBuilder::addInfo( LHCb::Particle& prt ) const { double cpx( 0 ), cpy( 0 ), cptMax( 0 ), nptMax( 0 ), width( 0 ), norm( 0 ), trks( 0 ), pt( prt.momentum().Pt() ); LHCb::Particle::ConstVector dtrs; LoKi::Extract::getParticles( &prt, back_inserter( dtrs ), LoKi::Cuts::HASPROTO && LoKi::Cuts::ISBASIC ); - for ( auto dtr : dtrs ) { + for ( auto* dtr : dtrs ) { const Gaudi::LorentzVector& vec = dtr->momentum(); if ( dtr->charge() != 0 ) { if ( vec.Pt() > cptMax ) cptMax = vec.Pt(); @@ -394,16 +388,11 @@ void FastJetBuilder::addInfo( LHCb::Particle& prt ) const { } // Relate particles with bestPVs -P2PVTable FastJetBuilder::relateP2Vtx( std::vector Particles, - const LHCb::RecVertices& PrimaryVertices ) const { - if ( !m_prtVtxTool ) { return nullptr; } - if ( PrimaryVertices.empty() ) { - if ( msgLevel( MSG::DEBUG ) ) debug() << "No PrimaryVertices in this event" << endmsg; - return nullptr; - } - auto vtx2prtRel = std::make_unique>(); - auto vtx = LHCb::VertexBase::ConstVector( PrimaryVertices.begin(), PrimaryVertices.end() ); - for ( auto prt : Particles ) { +P2PVTable FastJetBuilder::relateP2Vtx( LHCb::span Particles, + const LHCb::RecVertices& PrimaryVertices ) const { + auto vtx2prtRel = P2PVTable{}; + auto vtx = LHCb::VertexBase::ConstVector{PrimaryVertices.begin(), PrimaryVertices.end()}; + for ( auto* prt : Particles ) { auto bestPV = m_prtVtxTool->relatedPV( prt, vtx ); if ( !bestPV ) { warning() << "Could not find bestPV for particle " << prt << " Particle ID: " << prt->particleID().pid() @@ -413,17 +402,17 @@ P2PVTable FastJetBuilder::relateP2Vtx( std::vector Partic if ( msgLevel( MSG::DEBUG ) ) debug() << "====> Particle: " << prt << " Particle ID: " << prt->particleID().pid() << " Best PV: " << bestPV << " Chi2/nDofF: " << bestPV->chi2PerDoF() << " PV position: " << bestPV->position() << endmsg; - auto sc = vtx2prtRel->relate( static_cast( bestPV ), prt ); + auto sc = vtx2prtRel.relate( static_cast( bestPV ), prt ); if ( sc.isFailure() ) ++m_msg_noBestPV; } return vtx2prtRel; } // prepare FastJet -std::vector FastJetBuilder::prepare( const std::vector Particles ) const { +std::vector FastJetBuilder::prepare( LHCb::span particles ) const { std::vector jets; - jets.reserve( Particles.size() ); - for ( const auto& [i, p] : LHCb::range::enumerate( Particles ) ) { + jets.reserve( particles.size() ); + for ( const auto& [i, p] : LHCb::range::enumerate( particles ) ) { if ( !p ) { ++m_msg_invalid_input; continue; @@ -434,11 +423,12 @@ std::vector FastJetBuilder::prepare( const std::vector const& Particles, - const DetectorElement& geometry ) const { +std::vector> +FastJetBuilder::makeFastJets( LHCb::span particles, + const DetectorElement& geometry ) const { // prepare the input data and define the jets - auto jetinputs = prepare( Particles ); + auto jetinputs = prepare( particles ); if ( jetinputs.empty() ) { if ( msgLevel( MSG::DEBUG ) ) { debug() << "No particle inputs to reconstruct jets" << endmsg; } @@ -458,7 +448,7 @@ LHCb::Particle::Vector FastJetBuilder::makeFastJets( std::vector> recojets; recojets.reserve( jets.size() ); for ( const PJet& jet : jets ) { @@ -473,7 +463,7 @@ LHCb::Particle::Vector FastJetBuilder::makeFastJets( std::vectorcombine( daughters, pJet, vJet, geometry ); - pJet.setParticleID( LHCb::ParticleID( m_jetID ) ); if ( sc.isFailure() ) { ++m_msg_momcomb; continue; } + pJet.setParticleID( LHCb::ParticleID( m_jetID ) ); // redefine the momentum with the FastJet output - pJet.setMomentum( Gaudi::LorentzVector( jet.px(), jet.py(), jet.pz(), jet.e() ) ); - recojets.push_back( pJet.clone() ); + pJet.setMomentum( {jet.px(), jet.py(), jet.pz(), jet.e()} ); + recojets.emplace_back( pJet.clone() ); } return recojets; } -- GitLab From 5ea16942227010ce0580ef68dd097bf726af8162 Mon Sep 17 00:00:00 2001 From: Gerhard Raven Date: Sat, 28 Aug 2021 00:37:35 +0200 Subject: [PATCH 3/3] delay clone --- Phys/JetAccessories/src/FastJetBuilder.cpp | 47 +++++++++++----------- 1 file changed, 23 insertions(+), 24 deletions(-) diff --git a/Phys/JetAccessories/src/FastJetBuilder.cpp b/Phys/JetAccessories/src/FastJetBuilder.cpp index 4389c989b..5a521bf7b 100644 --- a/Phys/JetAccessories/src/FastJetBuilder.cpp +++ b/Phys/JetAccessories/src/FastJetBuilder.cpp @@ -220,8 +220,8 @@ private: mutable Gaudi::Accumulators::SummingCounter<> m_nbInput{this, "Nb of Inputs"}; mutable Gaudi::Accumulators::SummingCounter<> m_nbJetCounter{this, "Nb of Reconstructed Jets"}; - mutable Gaudi::Accumulators::SummingCounter<> m_nbJetCounter10{this, "Nb of Reconstructed Jets with pT>=10GeV"}; - mutable Gaudi::Accumulators::SummingCounter<> m_nbJetCounter20{this, "Nb of Reconstructed Jets with pT>=20GeV"}; + mutable Gaudi::Accumulators::SummingCounter<> m_nbJetCounter10{this, "Nb of Reconstructed Jets with pT>10GeV"}; + mutable Gaudi::Accumulators::SummingCounter<> m_nbJetCounter20{this, "Nb of Reconstructed Jets with pT>20GeV"}; /// MsgSvc counters mutable Gaudi::Accumulators::MsgCounter m_msg_noBestPV{this, @@ -240,8 +240,8 @@ private: void addInfo( LHCb::Particle& prt ) const; // run fastjet - std::vector> makeFastJets( LHCb::span particles, - const DetectorElement& geometry ) const; + std::vector makeFastJets( LHCb::span particles, + const DetectorElement& geometry ) const; // prepare fastjet definition std::vector prepare( LHCb::span Particles ) const; @@ -312,24 +312,24 @@ LHCb::Particles FastJetBuilder::operator()( const LHCb::Particles& Particles, co if ( msgLevel( MSG::DEBUG ) ) debug() << "Vtx: " << vtx << " #jets: " << tmprecojets.size() << endmsg; recojets.reserve( tmprecojets.size() ); for ( auto& jet : tmprecojets ) { - jet->addInfo( LHCb::JetIDInfo::vtx_x, vtx->position().X() ); - jet->addInfo( LHCb::JetIDInfo::vtx_y, vtx->position().Y() ); - jet->addInfo( LHCb::JetIDInfo::vtx_z, vtx->position().Z() ); - recojets.insert( jet.release() ); + jet.addInfo( LHCb::JetIDInfo::vtx_x, vtx->position().X() ); + jet.addInfo( LHCb::JetIDInfo::vtx_y, vtx->position().Y() ); + jet.addInfo( LHCb::JetIDInfo::vtx_z, vtx->position().Z() ); + recojets.insert( jet.clone() ); } } } else { // Make jets inclusively if ( msgLevel( MSG::DEBUG ) ) debug() << "Making jets inclusively" << endmsg; // Make the jets - auto tmprecojets = makeFastJets( std::vector{Particles.begin(), Particles.end()}, geometry ); + auto tmprecojets = makeFastJets( std::vector( Particles.begin(), Particles.end() ), geometry ); recojets.reserve( tmprecojets.size() ); for ( auto& jet : tmprecojets ) { - auto vtx = m_prtVtxTool->relatedPV( jet.get(), PrimaryVertices ); - jet->addInfo( LHCb::JetIDInfo::vtx_x, vtx->position().X() ); - jet->addInfo( LHCb::JetIDInfo::vtx_y, vtx->position().Y() ); - jet->addInfo( LHCb::JetIDInfo::vtx_z, vtx->position().Z() ); - recojets.insert( jet.release() ); + auto vtx = m_prtVtxTool->relatedPV( &jet, PrimaryVertices ); + jet.addInfo( LHCb::JetIDInfo::vtx_x, vtx->position().X() ); + jet.addInfo( LHCb::JetIDInfo::vtx_y, vtx->position().Y() ); + jet.addInfo( LHCb::JetIDInfo::vtx_z, vtx->position().Z() ); + recojets.insert( jet.clone() ); } } m_nbJetCounter += recojets.size(); @@ -339,8 +339,8 @@ LHCb::Particles FastJetBuilder::operator()( const LHCb::Particles& Particles, co for ( auto* jet : recojets ) { addInfo( *jet ); const auto pt = jet->momentum().Pt(); - if ( pt >= 10000 ) { c10 += 1; } - if ( pt >= 20000 ) { c20 += 1; } + if ( pt > 10 * Gaudi::Units::GeV ) { c10 += 1; } + if ( pt > 20 * Gaudi::Units::GeV ) { c20 += 1; } } if ( msgLevel( MSG::DEBUG ) ) { @@ -423,9 +423,8 @@ std::vector FastJetBuilder::prepare( LHCb::span> -FastJetBuilder::makeFastJets( LHCb::span particles, - const DetectorElement& geometry ) const { +std::vector FastJetBuilder::makeFastJets( LHCb::span particles, + const DetectorElement& geometry ) const { // prepare the input data and define the jets auto jetinputs = prepare( particles ); @@ -448,7 +447,7 @@ FastJetBuilder::makeFastJets( LHCb::span particles, if ( msgLevel( MSG::DEBUG ) ) { debug() << "No jets from fastjet::ClusterSequence" << endmsg; } - std::vector> recojets; + std::vector recojets; recojets.reserve( jets.size() ); for ( const PJet& jet : jets ) { @@ -470,17 +469,17 @@ FastJetBuilder::makeFastJets( LHCb::span particles, ++m_msg_nojetdaughters; continue; } - LHCb::Particle pJet; - LHCb::Vertex vJet; - auto sc = m_combiner->combine( daughters, pJet, vJet, geometry ); + auto& pJet = recojets.emplace_back(); + LHCb::Vertex vJet; + auto sc = m_combiner->combine( daughters, pJet, vJet, geometry ); if ( sc.isFailure() ) { + recojets.pop_back(); ++m_msg_momcomb; continue; } pJet.setParticleID( LHCb::ParticleID( m_jetID ) ); // redefine the momentum with the FastJet output pJet.setMomentum( {jet.px(), jet.py(), jet.pz(), jet.e()} ); - recojets.emplace_back( pJet.clone() ); } return recojets; } -- GitLab