diff --git a/Phys/ParticleCombiners/include/CombKernel/ParticleCombiner.h b/Phys/ParticleCombiners/include/CombKernel/ParticleCombiner.h index 706fa455cfb70647ff03fa28477fab483b7c5f21..6e46c4e729945bc12a5f2dca89f598bb4dddf796 100644 --- a/Phys/ParticleCombiners/include/CombKernel/ParticleCombiner.h +++ b/Phys/ParticleCombiners/include/CombKernel/ParticleCombiner.h @@ -150,6 +150,7 @@ namespace Combiner { using base_t = typename base_helper::type; } // namespace types + template class CheckOverlap { public: /** @@ -169,50 +170,54 @@ namespace Combiner { * @return true if overlapping * @return false if no common ancestors */ - bool findOverlap( const LHCb::Particle* particle1, const LHCb::Particle* particle2 ) { + bool findOverlap( const LHCb::Particle* particle1, const LHCb::Particle* particle2, const unsigned cont1, + const unsigned cont2, const unsigned idx1, const unsigned idx2 ) { if ( !particle1 || !particle2 ) { throw GaudiException( "Particle is nullptr.", __PRETTY_FUNCTION__, StatusCode::FAILURE ); } if ( m_tool ) { return m_tool->foundOverlap( particle1, particle2 ); } - return findOverlapImpl( particle1, particle2 ); + return findOverlapImpl( cont1, cont2, idx1, idx2 ); + } + + template + void prepare( const std::tuple& inputs ) { + if ( m_tool ) return; + std::apply( + [&]( auto&&... input ) { + unsigned i = 0; + ( prepare( i++, input ), ... ); + }, + inputs ); + } + + template + void prepare( unsigned i, const InputT& input ) { + m_protos_range[i].reserve( input.size() ); + for ( const auto& p : input ) { + unsigned start = m_protos.size(); + collectAllProtos( p ); + unsigned end = m_protos.size(); + m_protos_range[i].emplace_back( start, end ); + } } private: /** * @brief Search for common ancestors among all the two particles' proto particles. */ - bool findOverlapImpl( const LHCb::Particle* particle1, const LHCb::Particle* particle2 ) { - - if ( particle1->isBasicParticle() && particle2->isBasicParticle() ) { - return findOverlapBasics( particle1->proto(), particle2->proto() ); - } - - m_protos.clear(); - collectAllProtos( particle1 ); - collectAllProtos( particle2 ); - - for ( auto proto1 = m_protos.begin(); proto1 != m_protos.end(); ++proto1 ) { - for ( auto proto2 = std::next( proto1 ); proto2 != m_protos.end(); ++proto2 ) { - if ( findOverlapBasics( *proto1, *proto2 ) ) { return true; } + bool findOverlapImpl( const unsigned container1, const unsigned container2, const unsigned particle1, + const unsigned particle2 ) { + const auto& range_p1 = m_protos_range[container1][particle1]; + const auto& range_p2 = m_protos_range[container2][particle2]; + for ( unsigned i = range_p1.first; i < range_p1.second; i++ ) { + const auto& proto1 = m_protos[i]; + for ( unsigned j = range_p2.first; j < range_p2.second; j++ ) { + if ( proto1 == m_protos[j] ) { return true; } } } return false; } - /** - * @brief Check if basic proto particles overlap. - * - * @return true if pointer to protoparticle is same or underlying track is same - * @return false else - */ - bool findOverlapBasics( const LHCb::ProtoParticle* proto1, const LHCb::ProtoParticle* proto2 ) const { - // We should only perform a track comparison with non-null track pointers, - // else we will determine neutrals to overlap with each other (they have nullptr tracks) - // to avoid this it's enough to check one proto for nullptr track - return ( proto1 == proto2 ) || ( proto1->track() && ( proto1->track() == proto2->track() ) ) || - ( proto1->neutralPID() && ( proto1->neutralPID() == proto2->neutralPID() ) ); - } - /** * @brief Recursively get all leaves (proto particles) from particle's ancestors. * @@ -220,7 +225,16 @@ namespace Combiner { void collectAllProtos( const LHCb::Particle* particle ) { if ( particle->isBasicParticle() ) { assert( particle->proto() ); - m_protos.push_back( particle->proto() ); + const auto* proto = particle->proto(); + // Select correct pointer depending on the type, and encode the type + // in the unused lower bits of the pointer + uintptr_t ptr = reinterpret_cast( proto ); + if ( proto->track() ) { + ptr = ( reinterpret_cast( proto->track() ) & ~3 ) | 1; + } else if ( proto->neutralPID() ) { + ptr = ( reinterpret_cast( proto->neutralPID() ) & ~3 ) | 2; + } + m_protos.push_back( ptr ); } else { for ( const auto& p : particle->daughters() ) { collectAllProtos( p ); } } @@ -231,9 +245,9 @@ namespace Combiner { // storage for all proto particles related to the particles // one goal of this implementation is to reuse this memory instead of // allocating at every overlap check like the CheckOverlap tool - std::vector m_protos{}; + std::vector m_protos{}; + std::array>, NBody> m_protos_range{}; }; - } // namespace Combiner template @@ -377,10 +391,10 @@ private: template void process_decays( InputContainer const&... particles, Combiner::types::Output_t& storage, - CombinationCuts combination_cuts, CompositeCut composite_cut, IGeometryInfo const& geometry, - LHCb::PrimaryVertices const& pvs, Counters npassed_combination_cuts, Counter npassed_vertex_fit, - Counter npassed_composite_cut, std::index_sequence, - std::index_sequence ) const { + CombinationCuts const& combination_cuts, CompositeCut const& composite_cut, + IGeometryInfo const& geometry, LHCb::PrimaryVertices const& pvs, + Counters npassed_combination_cuts, Counter npassed_vertex_fit, Counter npassed_composite_cut, + std::index_sequence, std::index_sequence ) const { // Cannot use a structured binding (auto& [ps, vs] = storage) because Clang // cannot compile the vertex fit lambda when it uses a variable from a // structured binding (see LLVM bug #39963) @@ -498,9 +512,9 @@ private: IGeometryInfo const& geometry, std::index_sequence, std::index_sequence, combinatorics& /* has to be passed by reference! */ per_event ) const { - // Indicies into the first container + // Indices into the first container auto const& indices0 = matching_indices[0].second.get(); - // Indicies into the second container + // Indices into the second container auto const& indices1 = matching_indices[1].second.get(); // Whether we should do a 'triangle loop' to generate all 2-body combinations of the above indicies bool const do_triangle_loop = matching_indices[1].first; @@ -511,7 +525,7 @@ private: this->debug() << "indices0.size() == " << indices0.size() << ", indices1.size() == " << indices1.size() << ", do_triangle_loop == " << do_triangle_loop << endmsg; } - auto check_overlap = Combiner::CheckOverlap{ m_check_overlap.get() }; + auto check_overlap = Combiner::CheckOverlap{ m_check_overlap.get() }; // We will need these objects for all combinations passing the overlap cut, // which is likely to be most combinations, so create them once outside the // loop @@ -519,6 +533,7 @@ private: auto& npassed_comb_cut = std::get<0>( npassed_combination_cuts ); auto const inputs = std::tuple{ particles... }; + check_overlap.prepare( inputs ); for ( auto i0 = 0u; i0 < indices0.size(); ++i0 ) { // Load the particles matching the first child of the current decay // descriptor from the first container @@ -532,7 +547,7 @@ private: this->debug() << "Processing 2-body combination; indices = " << combination_indices << endmsg; } - const auto has_overlap = check_overlap.findOverlap( p0, p1 ); + const auto has_overlap = check_overlap.findOverlap( p0, p1, 0, 1, p0_index, p1_index ); if ( this->msgLevel( MSG::DEBUG ) ) { this->debug() << "Overlap check " << ( has_overlap ? "failed" : "passed" ) << endmsg; @@ -588,7 +603,7 @@ private: std::array const& current_combination_indices, VertexFit const& do_vertex_fit, CombinationCuts const& combination_cuts, CombinationCutCounters& npassed_combination_cuts, IGeometryInfo const& geometry, std::index_sequence, - combinatorics& /*by ref!*/ per_event, Combiner::CheckOverlap& check_overlap ) const { + combinatorics& /*by ref!*/ per_event, Combiner::CheckOverlap& check_overlap ) const { // Combination indices are the indices into input containers of the elements // of the current combination // NBodies is the size of the combination we've been given and that we will @@ -650,7 +665,14 @@ private: // We only need to check the overlap between each child in the current // N-body object and the current particle (i.e. we don't need to check // the overlap *within* the N-body object; that has already been done) - const auto has_overlap = ( check_overlap.findOverlap( &combination.template get(), pN1 ) || ... ); + const auto has_overlap = ( [&]() -> bool { + for ( unsigned c = 0; c < extended_combination_indices.size() - 1; c++ ) { + if ( check_overlap.findOverlap( &combination[c], pN1, c, NBodies, extended_combination_indices[c], + pN1_index ) ) + return true; + } + return false; + } )(); if ( this->msgLevel( MSG::DEBUG ) ) { this->debug() << "Overlap check " << ( has_overlap ? "failed" : "passed" ) << endmsg;