diff --git a/Phys/RelatedInfoTools/CMakeLists.txt b/Phys/RelatedInfoTools/CMakeLists.txt index f71a3106be74dc3593572ca8e0290488b15e91ed..1084efab3d1f0576530b8adf21d30fe29f8c42b6 100644 --- a/Phys/RelatedInfoTools/CMakeLists.txt +++ b/Phys/RelatedInfoTools/CMakeLists.txt @@ -22,6 +22,7 @@ gaudi_add_module(RelatedInfoTools src/RelInfoCylVariables.cpp src/RelInfoPFVariables.cpp src/RelInfoVertexIsolation.cpp + src/ProbeAndLongMatcher.cpp LINK Boost::headers Gaudi::GaudiAlgLib diff --git a/Phys/RelatedInfoTools/src/ProbeAndLongMatcher.cpp b/Phys/RelatedInfoTools/src/ProbeAndLongMatcher.cpp new file mode 100644 index 0000000000000000000000000000000000000000..247a575615d061adf0a1359741c4965827216b6a --- /dev/null +++ b/Phys/RelatedInfoTools/src/ProbeAndLongMatcher.cpp @@ -0,0 +1,161 @@ +/*****************************************************************************\ + * (c) Copyright 2000-2021 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. * +\*****************************************************************************/ + +#include "Event/Particle.h" +#include "GaudiAlg/Transformer.h" +#include "Kernel/HashIDs.h" +#include "Kernel/HeaderMapping.h" +#include "Kernel/LHCbID.h" + +/**** @ class ProbeAndLongMatcher + * + * Function to match the probe daughter of Jpsi composite to a long track + * by comparing the hits in sub-detector + * Needed in the study of Tracking efficiency + * + * + */ + +class ProbeAndLongMatcher + : public Gaudi::Functional::MultiTransformerFilter( + LHCb::Particle::Range const&, LHCb::Particle::Range const& )> { + +public: + ProbeAndLongMatcher( const std::string& name, ISvcLocator* pSvcLocator ); + + std::tuple + operator()( LHCb::Particle::Range const& composites, LHCb::Particle::Range const& longparts ) const override; + +private: + // properties are members of the class + // need to specify owner, name, default value, doc + Gaudi::Property m_thresMuon{this, "MinMuonFrac", 0.4}; // overlap threshold of Muon stations + Gaudi::Property m_thresVP{this, "MinVPFrac", 0.4}; // overlap threshold of VP + Gaudi::Property m_thresUT{this, "MinUTFrac", 0.4}; // overlap threshold of UT + Gaudi::Property m_thresFT{this, "MinFTFrac", 0.4}; // overlap threshold of FT + + Gaudi::Property p_checkMuon{this, "checkMuon", false}; + Gaudi::Property p_checkVP{this, "checkVP", false}; + Gaudi::Property p_checkUT{this, "checkUT", false}; + Gaudi::Property p_checkFT{this, "checkFT", false}; + + mutable Gaudi::Accumulators::StatCounter m_all{this, "#Input composites"}; // Number of JPsi candidates + mutable Gaudi::Accumulators::StatCounter m_matched{ + this, "#Matched composites"}; // Number of JPsi candidates with a matched probe track + mutable Gaudi::Accumulators::MsgCounter m_no_track{ + this, "Charged ProtoParticle with no Track found. This is likely to be a bug."}; + mutable Gaudi::Accumulators::MsgCounter m_no_state{ + this, "Track has empty states. This is likely to be a bug."}; +}; + +DECLARE_COMPONENT( ProbeAndLongMatcher ) + +// Implementation +ProbeAndLongMatcher::ProbeAndLongMatcher( const std::string& name, ISvcLocator* pSvcLocator ) + : MultiTransformerFilter( name, pSvcLocator, {KeyValue{"TwoBodyComposites", ""}, KeyValue{"LongTracks", ""}}, + {KeyValue{"MatchedComposites", ""}, KeyValue{"MatchedLongTracks", ""}} ) {} + +std::tuple ProbeAndLongMatcher:: + operator() // we use a mulitransformer to define multiple outputs + ( LHCb::Particle::Range const& composites, LHCb::Particle::Range const& longparts ) const { + + // initialize counter + m_all += composites.size(); + + // output containers + auto results = std::tuple{}; + auto& [filter_passed, selectedComp, selectedLong] = results; + + for ( const auto composite : composites ) { + + for ( const auto& daughter : composite->daughters() ) { + + const LHCb::ProtoParticle* proto = daughter->proto(); + const LHCb::Track* track = proto->track(); + const auto lhcbids = track->lhcbIDs(); + + // Sanity check from Particle_to_Particle.cpp + if ( proto == nullptr ) { + ++m_no_track; + continue; + } + if ( track->states().empty() ) { + ++m_no_state; + continue; + } + + bool hasVELO = false; + bool hasUT = false; + bool hasFT = false; + for ( auto const id : lhcbids ) { + if ( id.isVP() ) hasVELO = true; + if ( id.isUT() ) hasUT = true; + if ( id.isFT() ) hasFT = true; + } + + // to remove the tag tracks from the daughters (maybe come up with a better way / flag or label? to identify the + // probe track??) + if ( ( !p_checkUT && hasVELO && hasFT ) || ( p_checkUT && hasVELO && hasUT && hasFT ) ) continue; + + const LHCb::Particle* tmp_long = nullptr; + double tmp_frac = 0; + + for ( const auto longpart : longparts ) { + + if ( daughter->charge() != longpart->charge() ) continue; + + // define LHCbIDs and muonPID for the long track + const LHCb::ProtoParticle* longproto = longpart->proto(); + const LHCb::Track* longtrack = longproto->track(); + const auto longIDs = longtrack->lhcbIDs(); + const auto muonpid = longproto->muonPID(); + std::vector muonids; + muonids.reserve( 10 ); + LHCb::HashIDs::lhcbIDs( muonpid, muonids ); + + // evaluate the overlap + std::pair fracVP( 0., 0. ), fracUT( 0., 0. ), fracFT( 0., 0. ), fracMuon( 0., 0. ); + + if ( p_checkVP ) fracVP = LHCb::HashIDs::overlap( lhcbids, longIDs, LHCb::LHCbID::channelIDtype::VP ); + if ( p_checkUT ) fracUT = LHCb::HashIDs::overlap( lhcbids, longIDs, LHCb::LHCbID::channelIDtype::UT ); + if ( p_checkFT ) fracFT = LHCb::HashIDs::overlap( lhcbids, longIDs, LHCb::LHCbID::channelIDtype::FT ); + if ( p_checkMuon ) fracMuon = LHCb::HashIDs::overlap( lhcbids, muonids, LHCb::LHCbID::channelIDtype::Muon ); + + if ( p_checkVP ) debug() << "match fraction of VP hits " << fracVP.first << " " << fracVP.second << endmsg; + if ( p_checkUT ) debug() << "match fraction of UT hits " << fracUT.first << " " << fracUT.second << endmsg; + if ( p_checkFT ) debug() << "match fraction of FT hits " << fracFT.first << " " << fracFT.second << endmsg; + if ( p_checkMuon ) + debug() << "match fraction of muon hits " << fracMuon.first << " " << fracMuon.second << endmsg; + + if ( p_checkVP && fracVP.first < m_thresVP.value() ) continue; + if ( p_checkUT && fracUT.first < m_thresUT.value() ) continue; + if ( p_checkFT && fracFT.first < m_thresFT.value() ) continue; + if ( p_checkMuon && fracMuon.first < m_thresMuon.value() ) continue; + + double frac = fracVP.first + fracUT.first + fracFT.first + fracMuon.first; + + if ( tmp_frac < frac ) { + tmp_frac = frac; + tmp_long = longpart; + } + } // end loop over long tracks + + if ( tmp_frac > 0 ) { + selectedComp.insert( composite ); + selectedLong.insert( tmp_long ); + } + } // end loop over daughters + } // end loop over composites + + filter_passed = selectedComp.size() > 0; + m_matched += selectedComp.size(); + return results; +}