diff --git a/Rich/RichFutureRecCheckers/CMakeLists.txt b/Rich/RichFutureRecCheckers/CMakeLists.txt index e558c30fd52855c7d888064ca79564a3c26749ff..8d31e8e4e219359696ecfb3c8916c41d56c8d7d9 100644 --- a/Rich/RichFutureRecCheckers/CMakeLists.txt +++ b/Rich/RichFutureRecCheckers/CMakeLists.txt @@ -23,17 +23,21 @@ gaudi_add_module(RichFutureRecCheckers src/RichPIDQC.cpp src/RichMCOpticalPhotons.cpp src/RichMCDetectorHits.cpp + src/RichPIDTupleCreatorAlg.cpp LINK Boost::headers Gaudi::GaudiAlgLib Gaudi::GaudiKernel Gaudi::GaudiUtilsLib + LHCb::LHCbAlgsLib LHCb::LHCbKernel LHCb::MCEvent + LHCb::MCInterfaces LHCb::RecEvent LHCb::RichFutureMCUtilsLib LHCb::RichUtils LHCb::TrackEvent + Rec::FunctorCoreLib Rec::RichFutureRecBase Rec::RichFutureRecEvent Rec::RichFutureRecInterfacesLib diff --git a/Rich/RichFutureRecCheckers/src/RichPIDTupleCreatorAlg.cpp b/Rich/RichFutureRecCheckers/src/RichPIDTupleCreatorAlg.cpp new file mode 100644 index 0000000000000000000000000000000000000000..a690c97ae36eb10b3bb3e0e8ff3775caa9a7cb3d --- /dev/null +++ b/Rich/RichFutureRecCheckers/src/RichPIDTupleCreatorAlg.cpp @@ -0,0 +1,236 @@ +/*****************************************************************************\ +* (c) Copyright 2000-2028 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. * +\*****************************************************************************/ +// base class +#include "GaudiAlg/GaudiTupleAlg.h" +#include "RichFutureRecBase/RichRecHistoAlgBase.h" +// Gaudi Functional +#include "LHCbAlgs/Consumer.h" + +// Event model +#include "Event/MCParticle.h" +#include "Event/RichPID.h" +#include "Event/Track.h" +#include "MCInterfaces/IRichMCTruthTool.h" + +// Relations +#include "RichFutureMCUtils/RichMCRelations.h" +#include "RichFutureMCUtils/TrackToMCParticle.h" + +// tools +#include "RichFutureRecInterfaces/IRichPIDPlots.h" +#include "TrackInterfaces/ITrackSelector.h" + +// Boost +#include "boost/lexical_cast.hpp" + +// STL +#include <algorithm> +#include <array> +#include <mutex> +#include <numeric> + +namespace Rich::Future::Rec::MC::Moni { + + /** @class RichPIDTupleAlgCreatorAlg.h + * + * Create RICH PID Ntuple for performance evaluation + * + * @author Sajan Easo + * @date 2024-11-28 + */ + + class RichPIDTupleCreatorAlg final + : public LHCb::Algorithm::Consumer<void( const LHCb::Track::Range&, // + const LHCb::RichPIDs&, // + const Rich::Future::MC::Relations::TkToMCPRels& ), + Gaudi::Functional::Traits::BaseClass_t<GaudiTupleAlg>> { + + public: + /// standard constructor + + RichPIDTupleCreatorAlg( const std::string& name, ISvcLocator* pSvcLocator ) + : Consumer( name, pSvcLocator, + {KeyValue{"TracksLocation", LHCb::TrackLocation::Default}, + KeyValue{"RichPIDsLocation", LHCb::RichPIDLocation::Default}, + KeyValue{"TrackToMCParticlesRelations", Rich::Future::MC::Relations::TrackToMCParticles}} ) {} + + public: + /// Functional operator + void operator()( const LHCb::Track::Range& tracks, // + const LHCb::RichPIDs& pids, // + const Rich::Future::MC::Relations::TkToMCPRels& rels ) const override; + + private: + /// Tools + /** Track selector. **/ + ToolHandle<const ITrackSelector> m_tkSelForPIDNtuple{this, "TrackSelector", "TrackSelector"}; + + private: + /// Enable PID Ntuple creation and filling + + Gaudi::Property<bool> m_fillRichPIDNtuple{this, "FillRichPIDNtuple", true, + "Enable the creation and filling of RICH PID Ntuple"}; + /// Allow reassign PID to below threshold + Gaudi::Property<bool> m_allowBTreassign{this, "AllowBTReassign", true, + "Allow PIDs to be reassigned as Below Threshold"}; + + private: + /// mutex lock + mutable std::mutex m_updatePIDNtupleLock; + + /// RICH Event ID rate + mutable Gaudi::Accumulators::BinomialCounter<> m_eventsWithID{this, "Event ID rate"}; + + /// RICH Track ID rate + mutable Gaudi::Accumulators::BinomialCounter<> m_tracksWithID{this, "Track ID rate"}; + + /// Reco PID reassigned as BT + mutable Gaudi::Accumulators::BinomialCounter<> m_recoPIDBT{this, "Reassigned Reco PID BT"}; + }; + +} // namespace Rich::Future::Rec::MC::Moni +using namespace Rich::Future::Rec::MC::Moni; + +void RichPIDTupleCreatorAlg::operator()( const LHCb::Track::Range& tracks, // + const LHCb::RichPIDs& pids, // + const Rich::Future::MC::Relations::TkToMCPRels& rels ) const { + + auto tksWithID = m_tracksWithID.buffer(); + unsigned int pidCount = 0; + + // count events with at least 1 PID object + m_eventsWithID += !pids.empty(); + + // track selector shortcut + const auto tkSel = m_tkSelForPIDNtuple.get(); + + // count total selected tracks, PIDs etc. + const auto nTks = std::count_if( tracks.begin(), tracks.end(), [&tkSel, &pids, &tksWithID]( const auto* tk ) { + const bool sel = tk && tkSel->accept( *tk ); + if ( sel ) { + tksWithID += std::any_of( pids.begin(), pids.end(), [&tk]( const auto id ) { return id && tk == id->track(); } ); + } + return sel; + } ); + + // info() << "Counted " << nTks << " tracks" << endmsg; + verbose() << "Counted " << nTks << " tracks" << endmsg; + + // helper for the track -> MCP relations + const Rich::Future::MC::Relations::TrackToMCParticle tkToMPCs( rels ); + // the lock + std::lock_guard lock( m_updatePIDNtupleLock ); + + // Create and fill RICH PID Ntuple + if ( m_fillRichPIDNtuple ) { + + // Loop over all PID results + for ( const auto pid : pids ) { + + // Track for this PID + const auto tk = pid->track(); + + // info() << *pid << endmsg; + debug() << *pid << endmsg; + + // is track selected + if ( !tkSel->accept( *tk ) ) { continue; } + // Is the track in the input list + // (if not skip, means we are running on a reduced track list with selection cuts) + if ( std::none_of( tracks.begin(), tracks.end(), [&tk]( const auto* t ) { return t == tk; } ) ) { continue; } + // Count PIDs and tracks + ++pidCount; + verbose() << " -> PID count = " << pidCount << endmsg; + + // Get best reco PID + auto bpid = pid->bestParticleID(); + // if below threshold, set as such + const bool reassignRecoBT = ( m_allowBTreassign && !pid->isAboveThreshold( bpid ) ); + if ( reassignRecoBT ) { + + // info() << " -> Reassigned RecoPID to BT" << endmsg; + verbose() << " -> Reassigned RecoPID to BT" << endmsg; + bpid = Rich::BelowThreshold; + } + m_recoPIDBT += reassignRecoBT; + + // info() << " -> Best Reco PID = " << bpid << endmsg; + verbose() << " -> Best Reco PID = " << bpid << endmsg; + + // Get the MCParticle range for this track + const auto mcPs = tkToMPCs.mcParticleRange( *tk ); + + // info() << " -> Matched to " << mcPs.size() << " MCParticles" << endmsg; + verbose() << " -> Matched to " << mcPs.size() << " MCParticles" << endmsg; + + // For now, filling once per associated MCParticle + + for ( const auto MC : mcPs ) { + // get the MCP + const auto mcP = MC.to(); + + Tuple tuple = nTuple( "RichPIDTuple", "RICH PID Information Ntuple" ); + + // Status Code for filling the ntuple + StatusCode sc = StatusCode::SUCCESS; + + // fill some track info in the ntuple + if ( sc ) sc = tuple->column( "TrackP", tk->p() ); + if ( sc ) sc = tuple->column( "TrackPt", tk->pt() ); + if ( sc ) sc = tuple->column( "TrackChi2PerDof", tk->chi2PerDoF() ); + if ( sc ) sc = tuple->column( "TrackNumDof", tk->nDoF() ); + if ( sc ) sc = tuple->column( "TrackType", static_cast<int>( tk->type() ) ); + if ( sc ) sc = tuple->column( "TrackHistory", static_cast<int>( tk->history() ) ); + if ( sc ) sc = tuple->column( "TrackGhostProb", tk->ghostProbability() ); + if ( sc ) sc = tuple->column( "TrackLikelihood", tk->likelihood() ); + if ( sc ) sc = tuple->column( "TrackCloneDist", tk->info( LHCb::Track::AdditionalInfo::CloneDist, 9e10 ) ); + // fill some RICH info in the ntuple + if ( sc ) sc = tuple->column( "RichDLLe", ( pid ? pid->particleDeltaLL( Rich::Electron ) : 0 ) ); + if ( sc ) sc = tuple->column( "RichDLLmu", ( pid ? pid->particleDeltaLL( Rich::Muon ) : 0 ) ); + if ( sc ) sc = tuple->column( "RichDLLpi", ( pid ? pid->particleDeltaLL( Rich::Pion ) : 0 ) ); + if ( sc ) sc = tuple->column( "RichDLLk", ( pid ? pid->particleDeltaLL( Rich::Kaon ) : 0 ) ); + if ( sc ) sc = tuple->column( "RichDLLp", ( pid ? pid->particleDeltaLL( Rich::Proton ) : 0 ) ); + if ( sc ) sc = tuple->column( "RichDLLd", ( pid ? pid->particleDeltaLL( Rich::Deuteron ) : 0 ) ); + if ( sc ) sc = tuple->column( "RichDLLbt", ( pid ? pid->particleDeltaLL( Rich::BelowThreshold ) : 0 ) ); + if ( sc ) sc = tuple->column( "RichUsedAero", ( pid ? pid->usedAerogel() : false ) ); + if ( sc ) sc = tuple->column( "RichUsedR1Gas", ( pid ? pid->usedRich1Gas() : false ) ); + if ( sc ) sc = tuple->column( "RichUsedR2Gas", ( pid ? pid->usedRich2Gas() : false ) ); + if ( sc ) sc = tuple->column( "RichAboveElThres", ( pid ? pid->electronHypoAboveThres() : false ) ); + if ( sc ) sc = tuple->column( "RichAboveMuThres", ( pid ? pid->muonHypoAboveThres() : false ) ); + if ( sc ) sc = tuple->column( "RichAbovePiThres", ( pid ? pid->pionHypoAboveThres() : false ) ); + if ( sc ) sc = tuple->column( "RichAboveKaThres", ( pid ? pid->kaonHypoAboveThres() : false ) ); + if ( sc ) sc = tuple->column( "RichAbovePrThres", ( pid ? pid->protonHypoAboveThres() : false ) ); + if ( sc ) sc = tuple->column( "RichAboveDeThres", ( pid ? pid->deuteronHypoAboveThres() : false ) ); + if ( sc ) sc = tuple->column( "RichBestPID", ( pid ? static_cast<int>( pid->bestParticleID() ) : -1 ) ); + + // MCParticle information + + if ( sc ) sc = tuple->column( "MCParticleType", mcP ? mcP->particleID().pid() : 0 ); + if ( sc ) sc = tuple->column( "MCParticleP", mcP ? mcP->p() : -99999 ); + if ( sc ) sc = tuple->column( "MCParticlePt", mcP ? mcP->pt() : -99999 ); + if ( sc ) sc = tuple->column( "MCVirtualMass", mcP ? mcP->virtualMass() : -99999 ); + + // write the tuple for this protoP + sc.andThen( [&] { return tuple->write(); } ) + .orThrow( "Failed to fill RICH PID ntuple", "RichPIDTupleCreatorAlg" ); + + } // end of MCParticle loop + + } // end pid loop + + } // end creation and filling of RICH PID ntuple +} +//----------------------------------------------------------------------------- + +// Declaration of the Algorithm Factory +DECLARE_COMPONENT( RichPIDTupleCreatorAlg ) + +//-----------------------------------------------------------------------------