diff --git a/Tf/TrackSys/python/TrackSys/PrUpgradeChecking.py b/Tf/TrackSys/python/TrackSys/PrUpgradeChecking.py index 73b62951b9b981b712098d4a90a8a78822296822..978d1b7ce229588db4fb88609a96acf6948f81d6 100755 --- a/Tf/TrackSys/python/TrackSys/PrUpgradeChecking.py +++ b/Tf/TrackSys/python/TrackSys/PrUpgradeChecking.py @@ -290,13 +290,6 @@ def PrUpgradeChecking(defTracks={}, seq.Members += [trassociator] GaudiSequencer("MCLinksTrSeq").Members += [seq] - from Configurables import LHCb__Converters__RecVertex__v2__fromVectorLHCbRecVertex as FromVectorLHCbRecVertex - vertexConverter = FromVectorLHCbRecVertex("VertexConverter") - vertexConverter.InputVerticesName = "Rec/Vertex/Vector/Primary" - vertexConverter.InputTracksName = defTracksConverted["Velo"]["Location"] - vertexConverter.OutputVerticesName = "Rec/Vertex/Primary" - GaudiSequencer("MCLinksTrSeq").Members += [vertexConverter] - from Configurables import PrTrackChecker, PrUTHitChecker from Configurables import LoKi__Hybrid__MCTool MCHybridFactory = LoKi__Hybrid__MCTool("MCHybridFactory") @@ -468,9 +461,7 @@ def PrUpgradeChecking(defTracks={}, from Configurables import PrimaryVertexChecker PVCheck = PrimaryVertexChecker("PVChecker") if not PVCheck.isPropertySet('inputVerticesName'): - PVCheck.inputVerticesName = "Rec/Vertex/Primary" - if not PVCheck.isPropertySet('matchByTracks'): - PVCheck.matchByTracks = False + PVCheck.inputVerticesName = "Rec/Vertex/Vector/Primary" if not PVCheck.isPropertySet('nTracksToBeRecble'): PVCheck.nTracksToBeRecble = 4 if not PVCheck.isPropertySet('inputTracksName'): diff --git a/Tr/PatChecker/CMakeLists.txt b/Tr/PatChecker/CMakeLists.txt index 257b46b2ce2cd4ea38d93bf1265de44bafa9d9c0..812a343f44cc7f2654a69319d77864f6c6f984ed 100644 --- a/Tr/PatChecker/CMakeLists.txt +++ b/Tr/PatChecker/CMakeLists.txt @@ -13,7 +13,8 @@ ################################################################################ gaudi_subdir(PatChecker v3r16p1) -gaudi_depends_on_subdirs(Event/DigiEvent +gaudi_depends_on_subdirs(Det/VPDet + Event/DigiEvent Event/LinkerEvent Event/MCEvent Event/PhysEvent @@ -32,5 +33,5 @@ include_directories(SYSTEM ${Boost_INCLUDE_DIRS} ${ROOT_INCLUDE_DIRS}) gaudi_add_module(PatChecker src/*.cpp INCLUDE_DIRS AIDA GSL Event/DigiEvent Kernel/MCInterfaces Tf/TfKernel - LINK_LIBRARIES GSL LinkerEvent MCEvent PhysEvent GaudiAlgLib GaudiKernel) + LINK_LIBRARIES GSL LinkerEvent MCEvent PhysEvent GaudiAlgLib GaudiKernel VPDetLib) diff --git a/Tr/PatChecker/src/PrimaryVertexChecker.cpp b/Tr/PatChecker/src/PrimaryVertexChecker.cpp index f15f5b97b344e0a0aa96120bec09240eddd0bc0d..0e3b53b9a19e6203cf94ef225960c7524821b728 100644 --- a/Tr/PatChecker/src/PrimaryVertexChecker.cpp +++ b/Tr/PatChecker/src/PrimaryVertexChecker.cpp @@ -8,217 +8,250 @@ * granted to it by virtue of its status as an Intergovernmental Organization * * or submit itself to any jurisdiction. * \*****************************************************************************/ -// Include files -// from Gaudi + #include "AIDA/IHistogram1D.h" + #include "DetDesc/Condition.h" +#include "DetDesc/ConditionAccessorHolder.h" + +#include "Event/MCHeader.h" +#include "Event/MCParticle.h" +#include "Event/MCProperty.h" +#include "Event/MCTrackInfo.h" +#include "Event/MCVertex.h" + #include "Event/RecVertex.h" +#include "Event/RecVertex_v2.h" #include "Event/State.h" #include "Event/Track.h" +#include "Event/Track_v2.h" + +#include "GaudiAlg/Consumer.h" +#include "GaudiAlg/GaudiTupleAlg.h" #include "GaudiAlg/Tuples.h" + #include "GaudiKernel/PhysicalConstants.h" #include "GaudiKernel/SystemOfUnits.h" #include "GaudiUtils/HistoStats.h" -#include "Linker/AllLinks.h" -#include "Linker/LinkerWithKey.h" -#include <Event/MCTrackInfo.h> -#include <Linker/LinkedTo.h> -// local -#include "PrimaryVertexChecker.h" +#include "Kernel/STLExtensions.h" -bool sortmlt( MCPVInfo first, MCPVInfo second ) { return first.nRecTracks > second.nRecTracks; } +#include <string> +#include <tuple> +#include <vector> //----------------------------------------------------------------------------- // Implementation file for class : PrimaryVertexChecker //----------------------------------------------------------------------------- -// Declaration of the Algorithm Factory -DECLARE_COMPONENT( PrimaryVertexChecker ) +namespace { + + struct MCPVInfo { + LHCb::MCVertex* pMCPV; // pointer to MC PV + int nRecTracks; // number of reconstructed tracks from this MCPV + int nRecBackTracks; // number of reconstructed backward tracks + int indexRecPVInfo; // index to reconstructed PVInfo (-1 if not reco) + int nCorrectTracks; // correct tracks belonging to reconstructed PV + int multClosestMCPV; // multiplicity of closest reconstructable MCPV + double distToClosestMCPV; // distance to closest reconstructible MCPV + bool decayCharm; // type of mother particle + bool decayBeauty; + bool decayStrange; + std::vector<LHCb::MCParticle*> m_mcPartInMCPV; + std::vector<LHCb::Track*> m_recTracksInMCPV; + }; + + struct RecPVInfo { + public: + int nTracks; // number of tracks + double chi2; + double nDoF; + int mother; + Gaudi::XYZPoint position; // position + Gaudi::XYZPoint positionSigma; // position sigmas + int indexMCPVInfo; // index to MCPVInfo + const LHCb::Event::v2::RecVertex* pRECPV; + }; + + inline const std::string beamSpotCond = "/dd/Conditions/Online/Velo/MotionSystem"; + + struct Beamline_t { + double X = std::numeric_limits<double>::signaling_NaN(); + double Y = std::numeric_limits<double>::signaling_NaN(); + Beamline_t( Condition const& c ) + : X{( c.param<double>( "ResolPosRC" ) + c.param<double>( "ResolPosLA" ) ) / 2} + , Y{c.param<double>( "ResolPosY" )} {} + }; + + bool sortmlt( MCPVInfo const& first, MCPVInfo const& second ) { return first.nRecTracks > second.nRecTracks; } + + enum class recoAs { + all, + isolated, + close, + ntracks_low, + ntracks_high, + z_low, + z_middle, + z_high, + beauty, + charm, + strange, + other, + first, + second, + third, + fourth, + fifth + }; + + constexpr auto All = std::array{ + recoAs::all, recoAs::isolated, recoAs::close, recoAs::ntracks_low, recoAs::ntracks_high, recoAs::z_low, + recoAs::z_middle, recoAs::z_high, recoAs::beauty, recoAs::charm, recoAs::strange, recoAs::other, + recoAs::first, recoAs::second, recoAs::third, recoAs::fourth, recoAs::fifth}; + constexpr auto Part = std::array{recoAs::all, recoAs::beauty, recoAs::charm, recoAs::strange, recoAs::other}; + constexpr auto Basic = std::array{recoAs::all, recoAs::isolated, recoAs::close, recoAs::ntracks_low, + recoAs::ntracks_high, recoAs::z_low, recoAs::z_middle, recoAs::z_high, + recoAs::beauty, recoAs::charm, recoAs::strange, recoAs::other}; + + constexpr int size_basic = static_cast<int>( recoAs::other ) + 1; + constexpr int size_recoAs = static_cast<int>( recoAs::fifth ) + 1; + constexpr int size_multi = size_recoAs - size_basic; + constexpr int begin_multi = static_cast<int>( recoAs::first ); + +} // namespace + +class PrimaryVertexChecker + : public Gaudi::Functional::Consumer<void( const LHCb::Tracks&, LHCb::Event::v2::RecVertices const&, + const LHCb::MCVertices&, const LHCb::MCParticles&, const LHCb::MCHeader&, + const LHCb::MCProperty&, const Beamline_t& ), + LHCb::DetDesc::usesBaseAndConditions<GaudiTupleAlg, Beamline_t>> { +public: + /// Standard constructor + PrimaryVertexChecker( const std::string& name, ISvcLocator* pSvcLocator ) + : Consumer{name, + pSvcLocator, + { + KeyValue{"inputTracksName", LHCb::TrackLocation::Default}, + KeyValue{"inputVerticesName", LHCb::Event::v2::RecVertexLocation::Primary}, + KeyValue{"MCVertexInput", LHCb::MCVertexLocation::Default}, + KeyValue{"MCParticleInput", LHCb::MCParticleLocation::Default}, + KeyValue{"MCHeaderLocation", LHCb::MCHeaderLocation::Default}, + KeyValue{"MCPropertyInput", LHCb::MCPropertyLocation::TrackInfo}, + KeyValue{"BeamSpotLocation", "AlgorithmSpecific-" + name + "-beamspot"}, + }} {} + + StatusCode initialize() override; ///< Algorithm initialization + void operator()( const LHCb::Tracks& tracks, LHCb::Event::v2::RecVertices const& vertices, + const LHCb::MCVertices& mcvtx, const LHCb::MCParticles& mcps, const LHCb::MCHeader& mcheader, + const LHCb::MCProperty& flags, + const Beamline_t& ) const override; ///< Algorithm execution + StatusCode finalize() override; ///< Algorithm finalization + +private: + // bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); } + + Gaudi::Property<int> m_nTracksToBeRecble{this, "nTracksToBeRecble", 4}; // min number of tracks in PV + Gaudi::Property<bool> m_produceHistogram{this, "produceHistogram", true}; // producing histograms (light version) + Gaudi::Property<bool> m_produceNtuple{this, "produceNtuple", false}; // producing NTuples (full version) + Gaudi::Property<bool> m_requireVelo{this, "RequireVelo", true}; // requiring VELO for tracks + Gaudi::Property<double> m_dzIsolated{this, "dzIsolated", 10.0 * Gaudi::Units::mm}; // split close/isolated PVs + Gaudi::Property<int> m_nTracksToPrint{this, "nTracksToPrint", 10}; // split low/high multiplicity PVs + Gaudi::Property<double> m_zToPrint{this, "zToPrint", 50.0 * Gaudi::Units::mm}; // split in z + + mutable Gaudi::Accumulators::Counter<> m_nevt{this, "nEvents"}; + mutable std::map<recoAs, Gaudi::Accumulators::BinomialCounter<>> m_false; // False PVs vs Reco PVs + mutable std::map<recoAs, Gaudi::Accumulators::BinomialCounter<>> m_eff; // MC reconstructible vs Reconstructed + mutable std::map<recoAs, Gaudi::Accumulators::BinomialCounter<>> m_mcpv; // MC vs MC reconstrucible + mutable std::map<recoAs, Gaudi::Accumulators::AveragingCounter<>> m_av_mcp; // average mc particles in MCPV + mutable std::map<recoAs, Gaudi::Accumulators::AveragingCounter<>> m_av_tracks; // average tracks in RecoPV + + std::vector<MCPVInfo>::iterator closestMCPV( std::vector<MCPVInfo>& rblemcpv, const MCPVInfo& mc ) const; + // std::vector<MCPVInfo>::iterator& itmc ) const; + std::vector<MCPVInfo>::iterator closestMCPV( std::vector<MCPVInfo>& rblemcpv, const RecPVInfo& rec ) const; + + void collectProductss( LHCb::MCVertex* mcpv, LHCb::MCVertex* mcvtx, std::vector<LHCb::MCParticle*>& allprods ) const; + void printRat( const std::string name, const Gaudi::Accumulators::BinomialCounter<>& m_mcpv, + const Gaudi::Accumulators::BinomialCounter<>& m_eff, + const Gaudi::Accumulators::BinomialCounter<>& m_false ); + void printAvTracks( const std::string name, const Gaudi::Accumulators::AveragingCounter<>& m_av_tracks, + const Gaudi::Accumulators::AveragingCounter<>& m_av_mcp ); + void printRes( std::string mes, double x, double y, double z ); + std::string toString( const recoAs& n ) const; + bool checkCondition( const MCPVInfo& MCPV, const recoAs& n ) const; + + void match_mc_vertex_by_distance( int ipv, std::vector<RecPVInfo>& rinfo, std::vector<MCPVInfo>& mcpvvec ) const; + int count_velo_tracks( const std::vector<LHCb::Event::v2::WeightedTrack>& tracksPV ); + void count_reconstructible_mc_particles( std::vector<MCPVInfo>& mcpvvec, const LHCb::MCProperty& flags ) const; + double check_histogram( const AIDA::IHistogram1D* h, bool var ); +}; -//============================================================================= -// Standard constructor, initializes variables -//============================================================================= -PrimaryVertexChecker::PrimaryVertexChecker( const std::string& name, ISvcLocator* pSvcLocator ) - : GaudiTupleAlg( name, pSvcLocator ) { - declareProperty( "nTracksToBeRecble", m_nTracksToBeRecble = 5 ); - declareProperty( "produceNtuple", m_produceNtuple = false ); - declareProperty( "RequireVelo", m_requireVelo = true ); - declareProperty( "produceHistogram", m_produceHistogram = true ); - declareProperty( "dzIsolated", m_dzIsolated = 10.0 * Gaudi::Units::mm ); - declareProperty( "matchByTracks", m_matchByTracks = true ); - declareProperty( "inputTracksName", m_inputTracksName = LHCb::TrackLocation::Default ); - declareProperty( "inputVerticesName", m_inputVerticesName = LHCb::RecVertexLocation::Primary ); -} +// Declaration of the Algorithm Factory +DECLARE_COMPONENT_WITH_ID( PrimaryVertexChecker, "PrimaryVertexChecker" ) -//============================================================================= -// Initialization -//============================================================================= StatusCode PrimaryVertexChecker::initialize() { - StatusCode sc = GaudiTupleAlg::initialize(); // Must be executed first - if ( sc.isFailure() ) debug() << "==> Initialize" << endmsg; - - m_nevt = 0; - m_nMCPV = 0; - m_nRecMCPV = 0; - m_nMCPV_isol = 0; - m_nRecMCPV_isol = 0; - m_nMCPV_close = 0; - m_nRecMCPV_close = 0; - m_nFalsePV = 0; - m_nFalsePV_real = 0; - m_nMCPV_1mult = 0; - m_nRecMCPV_1mult = 0; - m_nMCPV_isol_1mult = 0; - m_nRecMCPV_isol_1mult = 0; - m_nMCPV_close_1mult = 0; - m_nRecMCPV_close_1mult = 0; - m_nMCPV_2mult = 0; - m_nRecMCPV_2mult = 0; - m_nMCPV_isol_2mult = 0; - m_nRecMCPV_isol_2mult = 0; - m_nMCPV_close_2mult = 0; - m_nRecMCPV_close_2mult = 0; - m_nMCPV_3mult = 0; - m_nRecMCPV_3mult = 0; - m_nMCPV_isol_3mult = 0; - m_nRecMCPV_isol_3mult = 0; - m_nMCPV_close_3mult = 0; - m_nRecMCPV_close_3mult = 0; - m_nBFalse = 0; - m_nRecBFalse = 0; - - return StatusCode::SUCCESS; + return Consumer::initialize().andThen( + [&] { addConditionDerivation<Beamline_t( Condition const& )>( {beamSpotCond}, inputLocation<Beamline_t>() ); } ); } //============================================================================= // Main execution //============================================================================= -StatusCode PrimaryVertexChecker::execute() { - debug() << "==> Execute" << endmsg; +void PrimaryVertexChecker::operator()( const LHCb::Tracks& tracks, LHCb::Event::v2::RecVertices const& recoVtx, + const LHCb::MCVertices& mcvtx, const LHCb::MCParticles& mcps, + const LHCb::MCHeader& mcheader, const LHCb::MCProperty& flags, + const Beamline_t& beamline ) const { + debug() << "==> Execute" << endmsg; // Event m_nevt++; - // Get input - std::vector<LHCb::Track*> vecOfTracks; - std::vector<LHCb::RecVertex*> vecOfVertices; - std::string trackLoc; - bool tracks_ok = getInputTracks( vecOfTracks, trackLoc ); - bool vertices_ok = getInputVertices( vecOfVertices ); - - if ( !tracks_ok && m_matchByTracks ) { - return StatusCode::SUCCESS; // return SUCCESSS anyway - } - - if ( !vertices_ok ) { return StatusCode::SUCCESS; } - - std::string m_beamSpotCond = "/dd/Conditions/Online/Velo/MotionSystem"; - // Condition *myCond = get<Condition>(detSvc(), m_beamSpotCond ); - const double xRC = 0.0; // myCond -> paramAsDouble ( "ResolPosRC" ) ; - const double xLA = 0.0; // myCond -> paramAsDouble ( "ResolPosLA" ) ; - const double Y = 0.0; // myCond -> paramAsDouble ( "ResolPosY" ) ; - - double m_beamSpotX = ( xRC + xLA ) / 2; - double m_beamSpotY = Y; - - if ( debugLevel() ) - debug() << trackLoc << " # tracks: " << vecOfTracks.size() << " # vertices: " << vecOfVertices.size() << endmsg; + if ( msgLevel( MSG::DEBUG ) ) + debug() << " # tracks: " << tracks.size() << " # vertices: " << recoVtx.size() << endmsg; // Fill reconstucted PV info - std::vector<RecPVInfo> recpvvec; - std::vector<LHCb::RecVertex*>::iterator itRecV; - for ( itRecV = vecOfVertices.begin(); vecOfVertices.end() != itRecV; itRecV++ ) { - LHCb::RecVertex* pv; - pv = *itRecV; + std::vector<RecPVInfo> recpvvec; + recpvvec.reserve( recoVtx.size() ); + + for ( const auto& pv : recoVtx ) { RecPVInfo recinfo; - recinfo.pRECPV = pv; - recinfo.position = pv->position(); - Gaudi::SymMatrix3x3 covPV = pv->covMatrix(); + recinfo.pRECPV = ( &pv ); + recinfo.position = pv.position(); + Gaudi::SymMatrix3x3 covPV = pv.covMatrix(); double sigx = sqrt( covPV( 0, 0 ) ); double sigy = sqrt( covPV( 1, 1 ) ); double sigz = sqrt( covPV( 2, 2 ) ); Gaudi::XYZPoint a3d( sigx, sigy, sigz ); recinfo.positionSigma = a3d; - recinfo.nTracks = pv->tracks().size(); - double minRD = 99999.; - double maxRD = -99999.; - double chi2 = pv->chi2(); - double nDoF = pv->nDoF(); - - std::string container = "container"; - typedef const SmartRefVector<LHCb::Track> PVTRACKS; - PVTRACKS& tracksIn = pv->tracks(); - PVTRACKS::const_iterator itin; - - int mother = 0; - int velo = 0; - int lg = 0; - double d0 = 0; - double mind0 = 99999.0; - double maxd0 = -99999.0; - double trackChi2 = 0.0; - int tr = 0; - - for ( itin = tracksIn.begin(); itin != tracksIn.end(); itin++ ) { - tr++; - const LHCb::Track* ptr = *itin; - - if ( ptr->type() == LHCb::Track::Types::Long ) { lg++; } - if ( ptr->type() == LHCb::Track::Types::Velo ) { velo++; } - - double x = ptr->position().x(); - double y = ptr->position().y(); - double r2 = x * x + y * y; - Gaudi::XYZVector udir = ptr->firstState().slopes().Unit(); - Gaudi::XYZVector distance = ptr->firstState().position() - pv->position(); - Gaudi::XYZVector d0i = udir.Cross( distance.Cross( udir ) ); - - d0 = d0 + d0i.Mag2(); - if ( d0i.Mag2() > maxd0 ) { maxd0 = d0i.Mag2(); } - if ( d0i.Mag2() < mind0 ) { mind0 = d0i.Mag2(); } - - if ( r2 > maxRD ) { maxRD = r2; } - if ( r2 < minRD ) { minRD = r2; } - - double trChi2 = ptr->chi2(); - trackChi2 = trackChi2 + trChi2; - } - - if ( m_matchByTracks ) { mother = check_mother_particle( pv, trackLoc ); } + recinfo.nTracks = pv.tracks().size(); + recinfo.chi2 = pv.chi2(); + recinfo.nDoF = pv.nDoF(); + recinfo.indexMCPVInfo = -1; - recinfo.minTrackRD = minRD; - recinfo.maxTrackRD = maxRD; + int mother = 0; recinfo.mother = mother; - recinfo.chi2 = chi2; - recinfo.nDoF = nDoF; - recinfo.d0 = d0; - recinfo.d0nTr = (double)d0 / (double)tr; - recinfo.chi2nTr = (double)trackChi2 / (double)tr; - recinfo.mind0 = mind0; - recinfo.maxd0 = maxd0; - recinfo.nVeloTracks = velo; - recinfo.nLongTracks = lg; recinfo.indexMCPVInfo = -1; recpvvec.push_back( recinfo ); } // Fill MC PV info std::vector<MCPVInfo> mcpvvec; - LHCb::MCVertices* mcvtx = get<LHCb::MCVertices>( LHCb::MCVertexLocation::Default ); - for ( LHCb::MCVertices::const_iterator itMCV = mcvtx->begin(); mcvtx->end() != itMCV; itMCV++ ) { - const LHCb::MCParticle* motherPart = ( *itMCV )->mother(); + mcpvvec.reserve( mcvtx.size() ); + + for ( const auto& mcpv : mcvtx ) { + const LHCb::MCParticle* motherPart = mcpv->mother(); if ( 0 == motherPart ) { - if ( ( *itMCV )->type() == LHCb::MCVertex::MCVertexType::ppCollision ) { + if ( mcpv->type() == LHCb::MCVertex::MCVertexType::ppCollision ) { MCPVInfo mcprimvert; - mcprimvert.pMCPV = *itMCV; + mcprimvert.pMCPV = mcpv; mcprimvert.nRecTracks = 0; mcprimvert.nRecBackTracks = 0; mcprimvert.indexRecPVInfo = -1; mcprimvert.nCorrectTracks = 0; mcprimvert.multClosestMCPV = 0; mcprimvert.distToClosestMCPV = 999999.; - mcprimvert.decayBeauty = 0; - mcprimvert.decayCharm = 0; + mcprimvert.decayBeauty = false; + mcprimvert.decayCharm = false; + mcprimvert.decayStrange = false; mcprimvert.m_mcPartInMCPV.clear(); mcprimvert.m_recTracksInMCPV.clear(); mcpvvec.push_back( mcprimvert ); @@ -226,356 +259,252 @@ StatusCode PrimaryVertexChecker::execute() { } } - const LHCb::MCParticle::Container* mcps = get<LHCb::MCParticle::Container>( LHCb::MCParticleLocation::Default ); - - int sMCP = 0; - int sTr = vecOfTracks.size(); - for ( LHCb::MCParticle::Container::const_iterator imc = mcps->begin(); mcps->end() != imc; ++imc ) { - const LHCb::MCParticle* dec = *imc; - double z = dec->originVertex()->position().z(); - if ( dec->particleID().threeCharge() != 0 && ( z < 400.0 && z > -600.0 ) ) { sMCP++; } - } - - if ( m_matchByTracks ) { - - count_reconstructed_tracks( mcpvvec, vecOfTracks, trackLoc ); - - } else { - - count_reconstructible_mc_particles( mcpvvec ); - } + count_reconstructible_mc_particles( mcpvvec, flags ); std::sort( mcpvvec.begin(), mcpvvec.end(), sortmlt ); std::vector<MCPVInfo> rblemcpv; + std::vector<MCPVInfo> notrblemcpv; std::vector<MCPVInfo> not_rble_but_visible; std::vector<MCPVInfo> not_rble; int nmrc = 0; std::vector<MCPVInfo>::iterator itmc; for ( itmc = mcpvvec.begin(); mcpvvec.end() != itmc; itmc++ ) { - rblemcpv.push_back( *itmc ); - if ( itmc->nRecTracks < m_nTracksToBeRecble ) { nmrc++; } + if ( itmc->nRecTracks >= m_nTracksToBeRecble ) { + rblemcpv.push_back( *itmc ); + } else { + notrblemcpv.push_back( *itmc ); + } if ( itmc->nRecTracks < m_nTracksToBeRecble && itmc->nRecTracks > 1 ) { not_rble_but_visible.push_back( *itmc ); } if ( itmc->nRecTracks < m_nTracksToBeRecble && itmc->nRecTracks < 2 ) { not_rble.push_back( *itmc ); } } - // match MC and REC PVs - if ( m_matchByTracks ) { - - for ( int ipv = 0; ipv < (int)recpvvec.size(); ipv++ ) { match_mc_vertex_by_tracks( ipv, recpvvec, rblemcpv ); } + for ( int ipv = 0; ipv < static_cast<int>( recpvvec.size() ); ipv++ ) { + match_mc_vertex_by_distance( ipv, recpvvec, rblemcpv ); + } - } else { + rblemcpv.insert( rblemcpv.end(), notrblemcpv.begin(), notrblemcpv.end() ); - for ( int ipv = 0; ipv < (int)recpvvec.size(); ipv++ ) { match_mc_vertex_by_distance( ipv, recpvvec, rblemcpv ); } + debug() << endmsg << " MC vertices " << endmsg; + debug() << " ===================================" << endmsg; + for ( const auto& [ipv, rmcpv] : LHCb::range::enumerate( rblemcpv ) ) { + std::string ff = " "; + LHCb::MCVertex* mcpv = rmcpv.pMCPV; + if ( rmcpv.indexRecPVInfo < 0 ) ff = " NOTRECO"; + debug() << format( " %3d %3d xyz ( %7.4f %7.4f %8.3f ) nrec = %4d", ipv, rmcpv.indexRecPVInfo, + mcpv->position().x(), mcpv->position().y(), mcpv->position().z(), rmcpv.nRecTracks ) + << ff << endmsg; } - if ( debugLevel() ) { - debug() << endmsg << " MC vertices " << endmsg; - debug() << " ===================================" << endmsg; - for ( int imcpv = 0; imcpv < (int)rblemcpv.size(); imcpv++ ) { - std::string ff = " "; - LHCb::MCVertex* mcpv = rblemcpv[imcpv].pMCPV; - if ( rblemcpv[imcpv].indexRecPVInfo < 0 ) ff = " NOTRECO"; - debug() << format( " %3d %3d xyz ( %7.4f %7.4f %8.3f ) nrec = %4d", imcpv, rblemcpv[imcpv].indexRecPVInfo, - mcpv->position().x(), mcpv->position().y(), mcpv->position().z(), rblemcpv[imcpv].nRecTracks ) - << ff << endmsg; - } - debug() << " -----------------------------------" << endmsg << endmsg; - - debug() << endmsg << " REC vertices " << endmsg; - debug() << " ===================================" << endmsg; - for ( int ipv = 0; ipv < (int)recpvvec.size(); ipv++ ) { - std::string ff = " "; - if ( recpvvec[ipv].indexMCPVInfo < 0 ) ff = " FALSE "; - debug() - << format( - " %3d %3d xyz ( %7.4f %7.4f %8.3f ) ntra = %4d sigxyz ( %7.4f %7.4f %8.4f ) chi2/NDF = %7.2f", - ipv, recpvvec[ipv].indexMCPVInfo, recpvvec[ipv].position.x(), recpvvec[ipv].position.y(), - recpvvec[ipv].position.z(), recpvvec[ipv].nTracks, recpvvec[ipv].nVeloTracks, - recpvvec[ipv].positionSigma.x(), recpvvec[ipv].positionSigma.y(), recpvvec[ipv].positionSigma.z(), - recpvvec[ipv].pRECPV->chi2PerDoF() ) - << ff << endmsg; - } - debug() << " -----------------------------------" << endmsg; + debug() << " -----------------------------------" << endmsg << endmsg; + + debug() << endmsg << " REC vertices " << endmsg; + debug() << " ===================================" << endmsg; + for ( const auto& [ipv, recpv] : LHCb::range::enumerate( recpvvec ) ) { + std::string ff = " "; + if ( recpvvec[ipv].indexMCPVInfo < 0 ) ff = " FALSE "; + debug() << format( + " %3d %3d xyz ( %7.4f %7.4f %8.3f ) ntra = %4d sigxyz ( %7.4f %7.4f %8.4f ) chi2/NDF = %7.2f", + ipv, recpv.indexMCPVInfo, recpv.position.x(), recpv.position.y(), recpv.position.z(), recpv.nTracks, + recpv.positionSigma.x(), recpv.positionSigma.y(), recpv.positionSigma.z(), + recpv.pRECPV->chi2PerDoF() ) + << ff << endmsg; } + debug() << " -----------------------------------" << endmsg; // find nr of false PV - - int nFalsePV = 0; int nFalsePV_real = 0; - for ( int ipv = 0; ipv < (int)recpvvec.size(); ipv++ ) { - int fake = 0; - double x = recpvvec[ipv].position.x(); - double y = recpvvec[ipv].position.y(); - double z = recpvvec[ipv].position.z(); - double r = std::sqrt( x * x + y * y ); - double errx = recpvvec[ipv].positionSigma.x(); - double erry = recpvvec[ipv].positionSigma.y(); - double errz = recpvvec[ipv].positionSigma.z(); - double errr = std::sqrt( ( ( x * errx ) * ( x * errx ) + ( y * erry ) * ( y * erry ) ) / ( x * x + y * y ) ); - double minRDTrack = recpvvec[ipv].minTrackRD; - double maxRDTrack = recpvvec[ipv].maxTrackRD; - int mother = recpvvec[ipv].mother; - double velo = recpvvec[ipv].nVeloTracks; - double lg = recpvvec[ipv].nLongTracks; - double d0 = recpvvec[ipv].d0; - double d0nTr = recpvvec[ipv].d0nTr; - double chi2nTr = recpvvec[ipv].chi2nTr; - double mind0 = recpvvec[ipv].mind0; - double maxd0 = recpvvec[ipv].maxd0; - double chi2 = recpvvec[ipv].chi2; - double nDoF = recpvvec[ipv].nDoF; - - if ( recpvvec[ipv].indexMCPVInfo < 0 ) { - nFalsePV++; - fake = 1; + for ( const auto& [ipv, recpv] : LHCb::range::enumerate( recpvvec ) ) { + int fake = 0; + double x = recpv.position.x(); + double y = recpv.position.y(); + double z = recpv.position.z(); + double r = std::sqrt( x * x + y * y ); + double errx = recpv.positionSigma.x(); + double erry = recpv.positionSigma.y(); + double errz = recpv.positionSigma.z(); + double errr = std::sqrt( ( ( x * errx ) * ( x * errx ) + ( y * erry ) * ( y * erry ) ) / ( x * x + y * y ) ); + int mother = recpv.mother; + double chi2 = recpv.chi2; + double nDoF = recpv.nDoF; + + if ( recpv.indexMCPVInfo < 0 ) { + fake = 1; + auto cmc = closestMCPV( rblemcpv, recpv ); + double dist = 999999.; + + if ( cmc != rblemcpv.end() ) { + dist = ( cmc->pMCPV->position() - recpv.pRECPV->position() ).R(); + + for ( const auto& n : Part ) m_false[n] += checkCondition( *cmc, n ); + + if ( dist > m_dzIsolated.value() ) { + m_false[recoAs::isolated] += true; + } else { + m_false[recoAs::close] += true; + } + if ( recpv.nTracks >= m_nTracksToPrint.value() ) { + m_false[recoAs::ntracks_high] += true; + } else { + m_false[recoAs::ntracks_low] += true; + } + if ( recpv.position.z() < -m_zToPrint.value() ) { + m_false[recoAs::z_low] += true; + } else if ( recpv.position.z() < m_zToPrint.value() ) { + m_false[recoAs::z_middle] += true; + } else { + m_false[recoAs::z_high] += true; + } + + int idx = std::distance( rblemcpv.begin(), cmc ); + if ( idx < size_multi ) { m_false[All[begin_multi + idx]] += true; } + } + bool vis_found = false; for ( unsigned int imc = 0; imc < not_rble_but_visible.size(); imc++ ) { if ( not_rble_but_visible[imc].indexRecPVInfo > -1 ) continue; - double dist = fabs( mcpvvec[imc].pMCPV->position().z() - recpvvec[ipv].position.z() ); - if ( dist < 5.0 * recpvvec[ipv].positionSigma.z() ) { + double dist = fabs( mcpvvec[imc].pMCPV->position().z() - recpv.position.z() ); + if ( dist < 5.0 * recpv.positionSigma.z() ) { vis_found = true; not_rble_but_visible[imc].indexRecPVInfo = 10; break; } } // imc + if ( !vis_found ) nFalsePV_real++; } if ( m_produceNtuple ) { Tuple myTuple2 = nTuple( 102, "PV_nTuple2", CLID_ColumnWiseTuple ); - myTuple2->column( "fake", double( fake ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "r", double( r ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "x", double( x ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "y", double( y ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "z", double( z ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "errr", double( errr ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "errz", double( errz ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "errx", double( errx ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "erry", double( erry ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "minRDTrack", double( minRDTrack ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "maxRDTrack", double( maxRDTrack ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "mother", double( mother ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "velo", double( velo ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "long", double( lg ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "d0", double( d0 ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "d0nTr", double( d0nTr ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "chi2nTr", double( chi2nTr ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "mind0", double( mind0 ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "maxd0", double( maxd0 ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "chi2", double( chi2 ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->column( "nDoF", double( nDoF ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple2->write().ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); + myTuple2->column( "fake", double( fake ) ).ignore(); + myTuple2->column( "r", double( r ) ).ignore(); + myTuple2->column( "x", double( x ) ).ignore(); + myTuple2->column( "y", double( y ) ).ignore(); + myTuple2->column( "z", double( z ) ).ignore(); + myTuple2->column( "errr", double( errr ) ).ignore(); + myTuple2->column( "errz", double( errz ) ).ignore(); + myTuple2->column( "errx", double( errx ) ).ignore(); + myTuple2->column( "erry", double( erry ) ).ignore(); + myTuple2->column( "mother", double( mother ) ).ignore(); + myTuple2->column( "chi2", double( chi2 ) ).ignore(); + myTuple2->column( "nDoF", double( nDoF ) ).ignore(); + myTuple2->write().ignore(); } } // Fill distance to closest recble MC PV and its multiplicity - std::vector<MCPVInfo>::iterator itmcl; - for ( itmcl = rblemcpv.begin(); rblemcpv.end() != itmcl; itmcl++ ) { - std::vector<MCPVInfo>::iterator cmc = closestMCPV( rblemcpv, itmcl ); + for ( auto& mcpv : rblemcpv ) { + std::vector<MCPVInfo>::iterator cmc = closestMCPV( rblemcpv, mcpv ); double dist = 999999.; int mult = 0; if ( cmc != rblemcpv.end() ) { - dist = ( cmc->pMCPV->position() - itmcl->pMCPV->position() ).R(); + dist = ( cmc->pMCPV->position() - mcpv.pMCPV->position() ).R(); mult = cmc->nRecTracks; } - itmcl->distToClosestMCPV = dist; - itmcl->multClosestMCPV = mult; - } - - // count non-reconstructible close and isolated PVs - int nmrc_isol = 0; - int nmrc_close = 0; - - // Counters - int nMCPV = rblemcpv.size() - nmrc; - int nRecMCPV = 0; - int nMCPV_isol = 0; - int nRecMCPV_isol = 0; - int nMCPV_close = 0; - int nRecMCPV_close = 0; - int nMCPV_1mult = 0; - int nRecMCPV_1mult = 0; - int nMCPV_isol_1mult = 0; - int nRecMCPV_isol_1mult = 0; - int nMCPV_close_1mult = 0; - int nRecMCPV_close_1mult = 0; - int nRecMCPV_wrong_1mult = 0; - int nMCPV_2mult = 0; - int nRecMCPV_2mult = 0; - int nMCPV_isol_2mult = 0; - int nRecMCPV_isol_2mult = 0; - int nMCPV_close_2mult = 0; - int nRecMCPV_close_2mult = 0; - int nRecMCPV_wrong_2mult = 0; - int nMCPV_3mult = 0; - int nRecMCPV_3mult = 0; - int nMCPV_isol_3mult = 0; - int nRecMCPV_isol_3mult = 0; - int nMCPV_close_3mult = 0; - int nRecMCPV_close_3mult = 0; - int nRecMCPV_wrong_3mult = 0; - - for ( itmc = rblemcpv.begin(); rblemcpv.end() != itmc; itmc++ ) { - if ( itmc->distToClosestMCPV > m_dzIsolated ) nMCPV_isol++; - if ( itmc->distToClosestMCPV > m_dzIsolated && itmc->nRecTracks < m_nTracksToBeRecble ) nmrc_isol++; - if ( itmc->distToClosestMCPV < m_dzIsolated ) nMCPV_close++; - if ( itmc->distToClosestMCPV < m_dzIsolated && itmc->nRecTracks < m_nTracksToBeRecble ) nmrc_close++; - - if ( itmc->indexRecPVInfo > -1 ) { - nRecMCPV++; - if ( itmc->distToClosestMCPV > m_dzIsolated ) nRecMCPV_isol++; - if ( itmc->distToClosestMCPV < m_dzIsolated ) nRecMCPV_close++; - } - } - - // rblemcpv is already sorted - - // highest mult - if ( rblemcpv.size() > 0 ) { - nMCPV_1mult++; - double dist = rblemcpv[0].distToClosestMCPV; - if ( dist > m_dzIsolated ) nMCPV_isol_1mult++; - if ( dist < m_dzIsolated ) nMCPV_close_1mult++; - if ( rblemcpv[0].indexRecPVInfo > -1 ) { - nRecMCPV_1mult++; - if ( dist > m_dzIsolated ) nRecMCPV_isol_1mult++; - if ( dist < m_dzIsolated ) nRecMCPV_close_1mult++; - } else { - nRecMCPV_wrong_1mult++; - } + mcpv.distToClosestMCPV = dist; + mcpv.multClosestMCPV = mult; } - // second highest - if ( rblemcpv.size() > 1 ) { - nMCPV_2mult++; - double dist = rblemcpv[1].distToClosestMCPV; - if ( dist > m_dzIsolated ) nMCPV_isol_2mult++; - if ( dist < m_dzIsolated ) nMCPV_close_2mult++; - if ( rblemcpv[1].indexRecPVInfo > -1 ) { - nRecMCPV_2mult++; - if ( dist > m_dzIsolated ) nRecMCPV_isol_2mult++; - if ( dist < m_dzIsolated ) nRecMCPV_close_2mult++; - } else { - nRecMCPV_wrong_2mult++; - } - } - // third highest - if ( rblemcpv.size() > 2 ) { - nMCPV_3mult++; - double dist = rblemcpv[2].distToClosestMCPV; - if ( dist > m_dzIsolated ) nMCPV_isol_3mult++; - if ( dist < m_dzIsolated ) nMCPV_close_3mult++; - if ( rblemcpv[2].indexRecPVInfo > -1 ) { - nRecMCPV_3mult++; - if ( dist > m_dzIsolated ) nRecMCPV_isol_3mult++; - if ( dist < m_dzIsolated ) nRecMCPV_close_3mult++; - } else { - nRecMCPV_wrong_3mult++; + for ( const auto& itmc : rblemcpv ) { + for ( const auto& n : Basic ) { + bool cut = false; + cut = checkCondition( itmc, n ); + if ( cut ) { + if ( itmc.nRecTracks < m_nTracksToBeRecble ) { + m_mcpv[n] += false; + } else { + m_mcpv[n] += true; + m_av_mcp[n] += itmc.nRecTracks; + if ( itmc.indexRecPVInfo < 0 ) { + m_eff[n] += false; + } else { + m_eff[n] += true; + m_false[n] += false; + } + } + if ( itmc.indexRecPVInfo > -1 ) { m_av_tracks[n] += recpvvec[itmc.indexRecPVInfo].nTracks; } + } } } - nMCPV_isol = nMCPV_isol - nmrc_isol; - nMCPV_close = nMCPV_close - nmrc_close; - - m_nMCPV += nMCPV; - m_nRecMCPV += nRecMCPV; - m_nMCPV_isol += nMCPV_isol; - m_nRecMCPV_isol += nRecMCPV_isol; - m_nMCPV_close += nMCPV_close; - m_nRecMCPV_close += nRecMCPV_close; - m_nFalsePV += nFalsePV; - m_nFalsePV_real += nFalsePV_real; - m_nMCPV_1mult += nMCPV_1mult; - m_nRecMCPV_1mult += nRecMCPV_1mult; - m_nMCPV_isol_1mult += nMCPV_isol_1mult; - m_nRecMCPV_isol_1mult += nRecMCPV_isol_1mult; - m_nMCPV_close_1mult += nMCPV_close_1mult; - m_nRecMCPV_close_1mult += nRecMCPV_close_1mult; - m_nRecMCPV_wrong_1mult += nRecMCPV_wrong_1mult; - m_nMCPV_2mult += nMCPV_2mult; - m_nRecMCPV_2mult += nRecMCPV_2mult; - m_nMCPV_isol_2mult += nMCPV_isol_2mult; - m_nRecMCPV_isol_2mult += nRecMCPV_isol_2mult; - m_nMCPV_close_2mult += nMCPV_close_2mult; - m_nRecMCPV_wrong_2mult += nRecMCPV_wrong_2mult; - m_nRecMCPV_close_2mult += nRecMCPV_close_2mult; - m_nMCPV_3mult += nMCPV_3mult; - m_nRecMCPV_3mult += nRecMCPV_3mult; - m_nMCPV_isol_3mult += nMCPV_isol_3mult; - m_nRecMCPV_isol_3mult += nRecMCPV_isol_3mult; - m_nMCPV_close_3mult += nMCPV_close_3mult; - m_nRecMCPV_close_3mult += nRecMCPV_close_3mult; - m_nRecMCPV_wrong_3mult += nRecMCPV_wrong_3mult; - + int mcpv = 0; int high = 0; - for ( itmc = rblemcpv.begin(); rblemcpv.end() != itmc; itmc++ ) { - double x = -99999.; - double y = -99999.; - double dx = -99999.; - double dy = -99999.; - double dz = -99999.; - double r = -99999.; - double zMC = -99999.; - double yMC = -99999.; - double xMC = -99999.; - double rMC = -99999.; - double z = -99999.; - double errx = -99999.; - double erry = -99999.; - double errz = -99999.; - double errr = -99999.; - double minRDTrack = 99999.; - double maxRDTrack = -99999.; - double chi2 = -999999.; - double nDoF = -999999.; - int indRec = itmc->indexRecPVInfo; - int reconstructed = 0; - int ntracks_pvrec = 0; - int nvelotracks_pvrec = 0; - int ntracks_pvmc = 0; - int dtrcks = 0; - int pvevt = 0; - int mother = 0; - int assoctrks = 0; - int nassoctrks = 0; - - zMC = itmc->pMCPV->position().z(); - yMC = itmc->pMCPV->position().y(); - xMC = itmc->pMCPV->position().x(); - rMC = std::sqrt( ( xMC - m_beamSpotX ) * ( xMC - m_beamSpotX ) + ( yMC - m_beamSpotY ) * ( yMC - m_beamSpotY ) ); + for ( auto& itmc : rblemcpv ) { + double x = -99999.; + double y = -99999.; + double dx = -99999.; + double dy = -99999.; + double dz = -99999.; + double r = -99999.; + double zMC = -99999.; + double yMC = -99999.; + double xMC = -99999.; + double rMC = -99999.; + double z = -99999.; + double errx = -99999.; + double erry = -99999.; + double errz = -99999.; + double errr = -99999.; + double chi2 = -999999.; + double nDoF = -999999.; + int indRec = itmc.indexRecPVInfo; + int reconstructed = 0; + int ntracks_pvrec = 0; + int ntracks_pvmc = 0; + int dtrcks = 0; + int pvevt = 0; + int mother = 0; + int assoctrks = 0; + int nassoctrks = 0; + + zMC = itmc.pMCPV->position().z(); + yMC = itmc.pMCPV->position().y(); + xMC = itmc.pMCPV->position().x(); + rMC = std::sqrt( ( xMC - beamline.X ) * ( xMC - beamline.X ) + ( yMC - beamline.Y ) * ( yMC - beamline.Y ) ); + + if ( mcpv < size_multi ) { + if ( itmc.nRecTracks < m_nTracksToBeRecble ) { + m_mcpv[All[begin_multi + mcpv]] += false; + } else { + m_mcpv[All[begin_multi + mcpv]] += true; + m_av_mcp[All[begin_multi + mcpv]] += itmc.nRecTracks; + if ( itmc.indexRecPVInfo < 0 ) { + m_eff[All[begin_multi + mcpv]] += false; + } else { + m_eff[All[begin_multi + mcpv]] += true; + m_false[All[begin_multi + mcpv]] += false; + } + } + } if ( indRec > -1 ) { high++; pvevt++; reconstructed = 1; - dx = recpvvec[indRec].position.x() - itmc->pMCPV->position().x(); - dy = recpvvec[indRec].position.y() - itmc->pMCPV->position().y(); - dz = recpvvec[indRec].position.z() - itmc->pMCPV->position().z(); + dx = recpvvec[indRec].position.x() - itmc.pMCPV->position().x(); + dy = recpvvec[indRec].position.y() - itmc.pMCPV->position().y(); + dz = recpvvec[indRec].position.z() - itmc.pMCPV->position().z(); x = recpvvec[indRec].position.x(); y = recpvvec[indRec].position.y(); z = recpvvec[indRec].position.z(); - // zMC = itmc->pMCPV->position().z(); - r = std::sqrt( ( x - m_beamSpotX ) * ( x - m_beamSpotX ) + ( y - m_beamSpotY ) * ( y - m_beamSpotY ) ); - errx = recpvvec[indRec].positionSigma.x(); - erry = recpvvec[indRec].positionSigma.y(); - errz = recpvvec[indRec].positionSigma.z(); - errr = std::sqrt( ( ( x * errx ) * ( x * errx ) + ( y * erry ) * ( y * erry ) ) / ( x * x + y * y ) ); - ntracks_pvrec = recpvvec[indRec].nTracks; - nvelotracks_pvrec = recpvvec[indRec].nVeloTracks; - ntracks_pvmc = itmc->pMCPV->products().size(); - dtrcks = ntracks_pvmc - ntracks_pvrec; - minRDTrack = recpvvec[indRec].minTrackRD; - maxRDTrack = recpvvec[indRec].maxTrackRD; - mother = recpvvec[indRec].mother; - chi2 = recpvvec[indRec].chi2; - nDoF = recpvvec[indRec].nDoF; - + r = std::sqrt( ( x - beamline.X ) * ( x - beamline.X ) + ( y - beamline.Y ) * ( y - beamline.Y ) ); + errx = recpvvec[indRec].positionSigma.x(); + erry = recpvvec[indRec].positionSigma.y(); + errz = recpvvec[indRec].positionSigma.z(); + errr = std::sqrt( ( ( x * errx ) * ( x * errx ) + ( y * erry ) * ( y * erry ) ) / ( x * x + y * y ) ); + ntracks_pvrec = recpvvec[indRec].nTracks; + ntracks_pvmc = itmc.pMCPV->products().size(); + dtrcks = ntracks_pvmc - ntracks_pvrec; + mother = recpvvec[indRec].mother; + chi2 = recpvvec[indRec].chi2; + nDoF = recpvvec[indRec].nDoF; + + if ( mcpv < size_multi ) { m_av_tracks[All[begin_multi + mcpv]] += recpvvec[indRec].nTracks; } // Filling histograms if ( m_produceHistogram ) { - plot( itmc->pMCPV->position().x(), 1001, "xmc", -0.25, 0.25, 50 ); - plot( itmc->pMCPV->position().y(), 1002, "ymc", -0.25, 0.25, 50 ); - plot( itmc->pMCPV->position().z(), 1003, "zmc", -20, 20, 50 ); + plot( itmc.pMCPV->position().x(), 1001, "xmc", -0.25, 0.25, 50 ); + plot( itmc.pMCPV->position().y(), 1002, "ymc", -0.25, 0.25, 50 ); + plot( itmc.pMCPV->position().z(), 1003, "zmc", -20, 20, 50 ); plot( recpvvec[indRec].position.x(), 1011, "xrd", -0.25, 0.25, 50 ); plot( recpvvec[indRec].position.y(), 1012, "yrd", -0.25, 0.25, 50 ); plot( recpvvec[indRec].position.z(), 1013, "zrd", -20, 20, 50 ); @@ -585,119 +514,109 @@ StatusCode PrimaryVertexChecker::execute() { plot( dx / errx, 1031, "pullx", -5., 5., 50 ); plot( dy / erry, 1032, "pully", -5., 5., 50 ); plot( dz / errz, 1033, "pullz", -5., 5., 50 ); + if ( itmc.nRecTracks < 10 ) { + plot( dx, 1101, "dx", -0.25, 0.25, 50 ); + plot( dy, 1102, "dy", -0.25, 0.25, 50 ); + plot( dz, 1103, "dz", -0.5, 0.5, 50 ); + } else if ( itmc.nRecTracks >= 10 && itmc.nRecTracks < 30 ) { + plot( dx, 1111, "dx", -0.25, 0.25, 50 ); + plot( dy, 1112, "dy", -0.25, 0.25, 50 ); + plot( dz, 1113, "dz", -0.5, 0.5, 50 ); + } else { + plot( dx, 1121, "dx", -0.25, 0.25, 50 ); + plot( dy, 1122, "dy", -0.25, 0.25, 50 ); + plot( dz, 1123, "dz", -0.5, 0.5, 50 ); + } + + if ( itmc.pMCPV->position().z() < -m_zToPrint ) { + plot( dx, 1201, "dx", -0.25, 0.25, 50 ); + plot( dy, 1202, "dy", -0.25, 0.25, 50 ); + plot( dz, 1203, "dz", -0.5, 0.5, 50 ); + } else if ( itmc.pMCPV->position().z() < m_zToPrint ) { + plot( dx, 1211, "dx", -0.25, 0.25, 50 ); + plot( dy, 1212, "dy", -0.25, 0.25, 50 ); + plot( dz, 1213, "dz", -0.5, 0.5, 50 ); + } else { + plot( dx, 1221, "dx", -0.25, 0.25, 50 ); + plot( dy, 1222, "dy", -0.25, 0.25, 50 ); + plot( dz, 1223, "dz", -0.5, 0.5, 50 ); + } + plot( double( ntracks_pvrec ), 1041, "ntracks_pvrec", 0., 150., 50 ); plot( double( dtrcks ), 1042, "mcrdtracks", 0., 150., 50 ); - plot( double( nvelotracks_pvrec ), 1043, "nvelotracks_pvrec", 0., 150., 50 ); if ( pvevt == 1 ) { plot( double( recpvvec.size() ), 1051, "nPVperEvt", -0.5, 5.5, 6 ); for ( int ipvrec = 0; ipvrec < (int)recpvvec.size(); ipvrec++ ) { assoctrks = assoctrks + recpvvec[ipvrec].nTracks; } - nassoctrks = vecOfTracks.size() - assoctrks; + nassoctrks = tracks.size() - assoctrks; plot( double( nassoctrks ), 1052, "nassoctrks", 0., 150., 50 ); } } } + mcpv++; int isolated = 0; - double dist_closest = itmc->distToClosestMCPV; - if ( dist_closest > m_dzIsolated ) { isolated = 1; } + double dist_closest = itmc.distToClosestMCPV; + if ( dist_closest > m_dzIsolated.value() ) { isolated = 1; } // Filling ntuple if ( m_produceNtuple ) { Tuple myTuple = nTuple( 101, "PV_nTuple", CLID_ColumnWiseTuple ); - myTuple->column( "reco", double( reconstructed ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "isol", double( isolated ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "ntracks", double( ntracks_pvrec ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "nvelotracks", double( nvelotracks_pvrec ) ) - .ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "nrectrmc", double( itmc->nRecTracks ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "dzclose", dist_closest ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "nmcpv", double( rblemcpv.size() ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "nmcallpv", double( mcpvvec.size() ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "nrecpv", double( recpvvec.size() ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "decayCharm", double( itmc->decayCharm ) ) - .ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "decayBeauty", double( itmc->decayBeauty ) ) - .ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "multi", double( high ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "dx", dx ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "dy", dy ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "dz", dz ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "x", x ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "y", y ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "r", r ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "zMC", zMC ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "yMC", yMC ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "xMC", xMC ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "rMC", rMC ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "z", z ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "xBeam", m_beamSpotX ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "yBeam", m_beamSpotY ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "errx", errx ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "erry", erry ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "errz", errz ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "errr", errr ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "minRDTrack", minRDTrack ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "maxRDTrack", maxRDTrack ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "mother", double( mother ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "evtnr", double( m_nevt ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "chi2", double( chi2 ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "nDoF", double( nDoF ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "size_tracks", double( sTr ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "size_mcp", double( sMCP ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->column( "mcpvrec", double( nmrc ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); - myTuple->write().ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); + myTuple->column( "reco", double( reconstructed ) ).ignore(); + myTuple->column( "isol", double( isolated ) ).ignore(); + myTuple->column( "ntracks", double( ntracks_pvrec ) ).ignore(); + myTuple->column( "nrectrmc", double( itmc.nRecTracks ) ).ignore(); + myTuple->column( "dzclose", dist_closest ).ignore(); + myTuple->column( "nmcpv", double( rblemcpv.size() ) ).ignore(); + myTuple->column( "mtruemcpv", double( mcheader.numOfPrimaryVertices() ) ).ignore(); + myTuple->column( "nmcallpv", double( mcpvvec.size() ) ).ignore(); + myTuple->column( "nrecpv", double( recpvvec.size() ) ).ignore(); + myTuple->column( "decayCharm", double( itmc.decayCharm ) ).ignore(); + myTuple->column( "decayBeauty", double( itmc.decayBeauty ) ).ignore(); + myTuple->column( "decayStrange", double( itmc.decayStrange ) ).ignore(); + myTuple->column( "multirec", double( high ) ).ignore(); + myTuple->column( "multimc", double( mcpv ) ).ignore(); + myTuple->column( "dx", dx ).ignore(); + myTuple->column( "dy", dy ).ignore(); + myTuple->column( "dz", dz ).ignore(); + myTuple->column( "x", x ).ignore(); + myTuple->column( "y", y ).ignore(); + myTuple->column( "r", r ).ignore(); + myTuple->column( "zMC", zMC ).ignore(); + myTuple->column( "yMC", yMC ).ignore(); + myTuple->column( "xMC", xMC ).ignore(); + myTuple->column( "rMC", rMC ).ignore(); + myTuple->column( "z", z ).ignore(); + myTuple->column( "xBeam", beamline.X ).ignore(); + myTuple->column( "yBeam", beamline.Y ).ignore(); + myTuple->column( "errx", errx ).ignore(); + myTuple->column( "erry", erry ).ignore(); + myTuple->column( "errz", errz ).ignore(); + myTuple->column( "errr", errr ).ignore(); + myTuple->column( "mother", double( mother ) ).ignore(); + myTuple->column( "evtnr", double( m_nevt.nEntries() ) ).ignore(); + myTuple->column( "chi2", double( chi2 ) ).ignore(); + myTuple->column( "nDoF", double( nDoF ) ).ignore(); + myTuple->column( "size_tracks", double( tracks.size() ) ).ignore(); + myTuple->column( "size_mcp", double( mcps.size() ) ).ignore(); + myTuple->column( "mcpvrec", double( nmrc ) ).ignore(); + myTuple->write().ignore(); } } - return StatusCode::SUCCESS; -} - -void PrimaryVertexChecker::match_mc_vertex_by_tracks( int ipv, std::vector<RecPVInfo>& rinfo, - std::vector<MCPVInfo>& mcpvvec ) { - - LHCb::RecVertex* invtx = rinfo[ipv].pRECPV; - typedef const SmartRefVector<LHCb::Track> PVTRACKS; - PVTRACKS& tracksIn = invtx->tracks(); - PVTRACKS::const_iterator itin; - - int indexmc = -1; - double ratiomax = 0.; - for ( int imc = 0; imc < (int)mcpvvec.size(); imc++ ) { - if ( mcpvvec[imc].indexRecPVInfo > -1 ) continue; - int ntrin = 0; - for ( std::vector<LHCb::Track*>::iterator itr = mcpvvec[imc].m_recTracksInMCPV.begin(); - itr != mcpvvec[imc].m_recTracksInMCPV.end(); itr++ ) { - for ( itin = tracksIn.begin(); itin != tracksIn.end(); itin++ ) { - const LHCb::Track* ptr = *itin; - if ( ptr == *itr ) { - ntrin++; - break; - } - } - } - double ratio = 1. * ntrin / tracksIn.size(); - if ( ratio > ratiomax ) { - ratiomax = ratio; - indexmc = imc; - } - } // imc - if ( ratiomax > 0.05 ) { - rinfo[ipv].indexMCPVInfo = indexmc; - mcpvvec[indexmc].indexRecPVInfo = ipv; - } + return; } void PrimaryVertexChecker::match_mc_vertex_by_distance( int ipv, std::vector<RecPVInfo>& rinfo, - std::vector<MCPVInfo>& mcpvvec ) { + std::vector<MCPVInfo>& mcpvvec ) const { double mindist = 999999.; int indexmc = -1; - for ( int imc = 0; imc < (int)mcpvvec.size(); imc++ ) { - if ( mcpvvec[imc].indexRecPVInfo > -1 ) continue; - double dist = fabs( mcpvvec[imc].pMCPV->position().z() - rinfo[ipv].position.z() ); + for ( const auto& [imc, mcpv] : LHCb::range::enumerate( mcpvvec ) ) { + if ( mcpv.indexRecPVInfo > -1 ) continue; + double dist = fabs( mcpv.pMCPV->position().z() - rinfo[ipv].position.z() ); if ( dist < mindist ) { mindist = dist; indexmc = imc; @@ -711,16 +630,16 @@ void PrimaryVertexChecker::match_mc_vertex_by_distance( int ipv, std::vector<Rec } } -std::vector<MCPVInfo>::iterator PrimaryVertexChecker::closestMCPV( std::vector<MCPVInfo>& rblemcpv, - std::vector<MCPVInfo>::iterator& itmc ) { +std::vector<MCPVInfo>::iterator PrimaryVertexChecker::closestMCPV( std::vector<MCPVInfo>& rblemcpv, + const MCPVInfo& mc ) const { - std::vector<MCPVInfo>::iterator itret = rblemcpv.end(); - double mindist = 999999.; + auto itret = rblemcpv.end(); + double mindist = 999999.; if ( rblemcpv.size() < 2 ) return itret; std::vector<MCPVInfo>::iterator it; for ( it = rblemcpv.begin(); it != rblemcpv.end(); it++ ) { - if ( it->pMCPV != itmc->pMCPV ) { - double dist = ( it->pMCPV->position() - itmc->pMCPV->position() ).R(); + if ( it->pMCPV != mc.pMCPV ) { + double dist = ( it->pMCPV->position() - mc.pMCPV->position() ).R(); if ( dist < mindist ) { mindist = dist; itret = it; @@ -730,20 +649,20 @@ std::vector<MCPVInfo>::iterator PrimaryVertexChecker::closestMCPV( std::vector<M return itret; } -void PrimaryVertexChecker::collectProductss( LHCb::MCVertex* mcpv, LHCb::MCVertex* mcvtx, - std::vector<LHCb::MCParticle*>& allprods ) { - - SmartRefVector<LHCb::MCParticle> daughters = mcvtx->products(); - SmartRefVector<LHCb::MCParticle>::iterator idau; - for ( idau = daughters.begin(); idau != daughters.end(); idau++ ) { - double dv2 = ( mcpv->position() - ( *idau )->originVertex()->position() ).Mag2(); - if ( dv2 > ( 100. * Gaudi::Units::mm ) * ( 100. * Gaudi::Units::mm ) ) continue; - LHCb::MCParticle* pmcp = *idau; - allprods.push_back( pmcp ); - SmartRefVector<LHCb::MCVertex> decays = ( *idau )->endVertices(); - SmartRefVector<LHCb::MCVertex>::iterator ivtx; - for ( ivtx = decays.begin(); ivtx != decays.end(); ivtx++ ) { collectProductss( mcpv, *ivtx, allprods ); } +std::vector<MCPVInfo>::iterator PrimaryVertexChecker::closestMCPV( std::vector<MCPVInfo>& rblemcpv, + const RecPVInfo& rec ) const { + auto itret = rblemcpv.end(); + double mindist = 999999.; + if ( rblemcpv.size() < 2 ) return itret; + std::vector<MCPVInfo>::iterator it; + for ( it = rblemcpv.begin(); it != rblemcpv.end(); it++ ) { + double dist = ( it->pMCPV->position() - rec.pRECPV->position() ).R(); + if ( dist < mindist ) { + mindist = dist; + itret = it; + } } + return itret; } //============================================================================= @@ -752,290 +671,204 @@ void PrimaryVertexChecker::collectProductss( LHCb::MCVertex* mcpv, LHCb::MCVerte StatusCode PrimaryVertexChecker::finalize() { debug() << "==> Finalize" << endmsg; + info() << " ************************************ " << endmsg; + info() << " MC PV is reconstructible if at least " << m_nTracksToBeRecble.value() << " tracks are reconstructed" + << endmsg; + info() << " MC PV is isolated if dz to closest reconstructible MC PV > " << m_dzIsolated.value() << " mm" << endmsg; + info() << " MC efficiency split by tracks with threshold: (" << m_nTracksToBeRecble.value() << "," + << m_nTracksToPrint.value() << "), >= " << m_nTracksToPrint.value() << endmsg; + info() << " MC efficiency split by z position: <-" << m_zToPrint.value() << ", (-" << m_zToPrint.value() << "," + << m_zToPrint.value() << "), >" << m_zToPrint.value() << endmsg; + std::string ff = "by distance"; + info() << " REC and MC vertices matched: " << ff << endmsg; - info() << " ============================================" << endmsg; - info() << " Efficiencies for reconstructible MC vertices: " << endmsg; - info() << " ============================================" << endmsg; + info() << " " << endmsg; + for ( const auto& n : All ) { printRat( toString( n ), m_mcpv[n], m_eff[n], m_false[n] ); } + info() << " " << endmsg; + for ( const auto& n : All ) { printAvTracks( toString( n ), m_av_tracks[n], m_av_mcp[n] ); } info() << " " << endmsg; - info() << " MC PV is reconstructible if at least " << m_nTracksToBeRecble << " tracks are reconstructed" << endmsg; - info() << " MC PV is isolated if dz to closest reconstructible MC PV > " << m_dzIsolated << " mm" << endmsg; - std::string ff = "by counting tracks"; - if ( !m_matchByTracks ) ff = "by dz distance"; - info() << " REC and MC vertices matched: " << ff << endmsg; + printRes( "1_res_all", check_histogram( histo( HistoID( 1021 ) ), true ), + check_histogram( histo( HistoID( 1022 ) ), true ), check_histogram( histo( HistoID( 1023 ) ), true ) ); + + printRes( "2_res_ntracks<10", check_histogram( histo( HistoID( 1101 ) ), true ), + check_histogram( histo( HistoID( 1102 ) ), true ), check_histogram( histo( HistoID( 1103 ) ), true ) ); + printRes( "3_res_ntracks(10,30)", check_histogram( histo( HistoID( 1111 ) ), true ), + check_histogram( histo( HistoID( 1112 ) ), true ), check_histogram( histo( HistoID( 1113 ) ), true ) ); + printRes( "4_res_ntracks>30", check_histogram( histo( HistoID( 1121 ) ), true ), + check_histogram( histo( HistoID( 1122 ) ), true ), check_histogram( histo( HistoID( 1123 ) ), true ) ); + + printRes( "5_res_z<-50", check_histogram( histo( HistoID( 1201 ) ), true ), + check_histogram( histo( HistoID( 1202 ) ), true ), check_histogram( histo( HistoID( 1203 ) ), true ) ); + printRes( "6_res_z(-50,50)", check_histogram( histo( HistoID( 1211 ) ), true ), + check_histogram( histo( HistoID( 1212 ) ), true ), check_histogram( histo( HistoID( 1213 ) ), true ) ); + printRes( "7_res_z>50", check_histogram( histo( HistoID( 1221 ) ), true ), + check_histogram( histo( HistoID( 1222 ) ), true ), check_histogram( histo( HistoID( 1223 ) ), true ) ); info() << " " << endmsg; - printRat( "All", m_nRecMCPV, m_nMCPV ); - printRat( "Isolated", m_nRecMCPV_isol, m_nMCPV_isol ); - printRat( "Close", m_nRecMCPV_close, m_nMCPV_close ); - printRat( "False rate", m_nFalsePV, m_nRecMCPV + m_nFalsePV ); - - if ( debugLevel() ) { printRat( "Real false rate", m_nFalsePV_real, m_nRecMCPV + m_nFalsePV_real ); } - - info() << endmsg; - printRat( "False PV as B", m_nRecBFalse, m_nBFalse ); - info() << endmsg; - - info() << " --------------------------------------------" << endmsg; - info() << " Substatistics: " << endmsg; - info() << " --------------------------------------------" << endmsg; - info() << " 1st PV (highest multiplicity): " << endmsg; - printRat( "All", m_nRecMCPV_1mult, m_nMCPV_1mult ); - printRat( "Isolated", m_nRecMCPV_isol_1mult, m_nMCPV_isol_1mult ); - printRat( "Close", m_nRecMCPV_close_1mult, m_nMCPV_close_1mult ); - - info() << " ---------------------------------------" << endmsg; - info() << " 2nd PV: " << endmsg; - printRat( "All", m_nRecMCPV_2mult, m_nMCPV_2mult ); - printRat( "Isolated", m_nRecMCPV_isol_2mult, m_nMCPV_isol_2mult ); - printRat( "Close", m_nRecMCPV_close_2mult, m_nMCPV_close_2mult ); - - info() << " ---------------------------------------" << endmsg; - info() << " 3rd PV: " << endmsg; - printRat( "All", m_nRecMCPV_3mult, m_nMCPV_3mult ); - printRat( "Isolated", m_nRecMCPV_isol_3mult, m_nMCPV_isol_3mult ); - printRat( "Close", m_nRecMCPV_close_3mult, m_nMCPV_close_3mult ); + printRes( "1_pull_width_all", check_histogram( histo( HistoID( 1031 ) ), true ), + check_histogram( histo( HistoID( 1032 ) ), true ), check_histogram( histo( HistoID( 1033 ) ), true ) ); + printRes( "1_pull_mean_all", check_histogram( histo( HistoID( 1031 ) ), false ), + check_histogram( histo( HistoID( 1032 ) ), false ), check_histogram( histo( HistoID( 1033 ) ), false ) ); info() << " " << endmsg; - if ( debugLevel() ) { - info() << " * Real false rate means: no visible MC PV within 5 sigma of REC PV." - << " Visible MC PV: 2 tracks reconstructed" << endmsg; - info() << " " << endmsg; - } - const AIDA::IHistogram1D* dx = histo( HistoID( 1021 ) ); - const AIDA::IHistogram1D* pullx = histo( HistoID( 1031 ) ); - const AIDA::IHistogram1D* dy = histo( HistoID( 1022 ) ); - const AIDA::IHistogram1D* pully = histo( HistoID( 1032 ) ); - const AIDA::IHistogram1D* dz = histo( HistoID( 1023 ) ); - const AIDA::IHistogram1D* pullz = histo( HistoID( 1033 ) ); - if ( dx ) { - info() << " ---------------------------------------" << endmsg; - info() << "dx: " - << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", dx->mean(), - Gaudi::Utils::HistoStats::meanErr( dx ), dx->rms(), Gaudi::Utils::HistoStats::rmsErr( dx ) ) - << endmsg; - } - if ( dy ) { - info() << "dy: " - << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", dy->mean(), - Gaudi::Utils::HistoStats::meanErr( dy ), dy->rms(), Gaudi::Utils::HistoStats::rmsErr( dy ) ) - << endmsg; - } - if ( dz ) { - info() << "dz: " - << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", dz->mean(), - Gaudi::Utils::HistoStats::meanErr( dz ), dz->rms(), Gaudi::Utils::HistoStats::rmsErr( dz ) ) - << endmsg; - } - info() << " ---------------------------------------" << endmsg; - if ( pullx ) { - info() << "pullx: " - << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", pullx->mean(), - Gaudi::Utils::HistoStats::meanErr( pullx ), pullx->rms(), - Gaudi::Utils::HistoStats::rmsErr( pullx ) ) - << endmsg; - } - if ( pully ) { - info() << "pully: " - << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", pully->mean(), - Gaudi::Utils::HistoStats::meanErr( pully ), pully->rms(), - Gaudi::Utils::HistoStats::rmsErr( pully ) ) - << endmsg; - } - if ( pullz ) { - info() << "pullz: " - << format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f", pullz->mean(), - Gaudi::Utils::HistoStats::meanErr( pullz ), pullz->rms(), - Gaudi::Utils::HistoStats::rmsErr( pullz ) ) - << endmsg; - } - info() << " ============================================" << endmsg; return GaudiTupleAlg::finalize(); // Must be called after all other actions } -void PrimaryVertexChecker::printRat( std::string mes, int a, int b ) { +void PrimaryVertexChecker::printRat( const std::string name, const Gaudi::Accumulators::BinomialCounter<>& m_mcpv, + const Gaudi::Accumulators::BinomialCounter<>& m_eff, + const Gaudi::Accumulators::BinomialCounter<>& m_false ) { + unsigned int len = 25; + std::string pmes = name; + if ( pmes.size() < len ) pmes.resize( len, ' ' ); + pmes += " : "; + info() << pmes + << format( "%8d from %8d (%8d-%-8d) [ %5.2f %], false %4d from reco. %8d (%8d+%-4d) [ %5.2f %] ", + m_eff.nTrueEntries(), m_eff.nEntries(), m_mcpv.nEntries(), m_mcpv.nFalseEntries(), + m_eff.efficiency() * 100, m_false.nTrueEntries(), m_false.nEntries(), m_false.nFalseEntries(), + m_false.nTrueEntries(), m_false.efficiency() * 100 ) + << endmsg; +} - double rat = 0.; - if ( b > 0 ) rat = 1.0 * a / b; +void PrimaryVertexChecker::printAvTracks( const std::string name, + const Gaudi::Accumulators::AveragingCounter<>& m_av_tracks, + const Gaudi::Accumulators::AveragingCounter<>& m_av_mcp ) { + unsigned int len = 25; + std::string pmes = name; + if ( pmes.size() < len ) pmes.resize( len, ' ' ); + pmes += " : "; + info() << pmes << format( "av. PV tracks: %6.2f [MC: %6.2f]", m_av_tracks.mean(), m_av_mcp.mean() ) << endmsg; +} - // reformat message - unsigned int len = 20; +void PrimaryVertexChecker::printRes( std::string mes, double x, double y, double z ) { + unsigned int len = 25; std::string pmes = mes; while ( pmes.length() < len ) { pmes += " "; } pmes += " : "; - - info() << pmes << format( " %6.3f ( %7d / %8d )", rat, a, b ) << endmsg; + info() << pmes << format( " x: %+5.3f, y: %+5.3f, z: %+5.3f", x, y, z ) << endmsg; } -void PrimaryVertexChecker::count_reconstructed_tracks( std::vector<MCPVInfo>& mcpvvec, - std::vector<LHCb::Track*>& vecOfTracks, std::string trackLoc ) { - - LinkedTo<LHCb::MCParticle> trackMClink( eventSvc(), msgSvc(), trackLoc ); - - // Find # of reconstructed tracks of every MC PV - std::vector<MCPVInfo>::iterator itinfomc; - for ( itinfomc = mcpvvec.begin(); mcpvvec.end() != itinfomc; itinfomc++ ) { - LHCb::MCVertex* avtx = itinfomc->pMCPV; - - SmartRefVector<LHCb::MCParticle> parts = avtx->products(); - std::vector<LHCb::MCParticle*> allproducts; - collectProductss( avtx, avtx, allproducts ); +double PrimaryVertexChecker::check_histogram( const AIDA::IHistogram1D* h, bool rms ) { + return h ? ( rms ? h->rms() : h->mean() ) : 0.; +} - LHCb::Track* recTrack = 0; - std::vector<LHCb::MCParticle*>::iterator imcp; - for ( imcp = allproducts.begin(); allproducts.end() != imcp; imcp++ ) { - LHCb::MCParticle* pmcp = *imcp; - int isReco = 0; - for ( std::vector<LHCb::Track*>::iterator itrec = vecOfTracks.begin(); itrec != vecOfTracks.end(); itrec++ ) { - LHCb::MCParticle* partmc = trackMClink.first( ( *itrec )->key() ); - if ( partmc && partmc == pmcp ) { - recTrack = ( *itrec ); - if ( !m_requireVelo || recTrack->hasVelo() ) { - isReco = 1; - break; - } - } - } - if ( pmcp->particleID().threeCharge() != 0 && isReco ) { - double dv2 = ( avtx->position() - pmcp->originVertex()->position() ).Mag2(); - if ( dv2 < 0.0001 ) { - itinfomc->m_mcPartInMCPV.push_back( pmcp ); - itinfomc->m_recTracksInMCPV.push_back( recTrack ); - } - } - } - itinfomc->nRecTracks = itinfomc->m_mcPartInMCPV.size(); - } +int PrimaryVertexChecker::count_velo_tracks( const std::vector<LHCb::Event::v2::WeightedTrack>& tracksPV ) { + return std::count_if( tracksPV.begin(), tracksPV.end(), []( const auto& t ) { return t.track->hasVelo(); } ); } -int PrimaryVertexChecker::count_velo_tracks( LHCb::RecVertex* RecVtx ) { - SmartRefVector<LHCb::Track> vtx_tracks = RecVtx->tracks(); - int nVeloTracks = 0; - for ( unsigned int it = 0; it < vtx_tracks.size(); it++ ) { - const LHCb::Track* ptr = vtx_tracks[it]; - if ( ptr->hasVelo() ) nVeloTracks++; +void PrimaryVertexChecker::collectProductss( LHCb::MCVertex* mcpv, LHCb::MCVertex* mcvtx, + std::vector<LHCb::MCParticle*>& allprods ) const { + + SmartRefVector<LHCb::MCParticle> daughters = mcvtx->products(); + SmartRefVector<LHCb::MCParticle>::iterator idau; + for ( idau = daughters.begin(); idau != daughters.end(); idau++ ) { + double dv2 = ( mcpv->position() - ( *idau )->originVertex()->position() ).Mag2(); + if ( dv2 > ( 100. * Gaudi::Units::mm ) * ( 100. * Gaudi::Units::mm ) ) continue; + LHCb::MCParticle* pmcp = *idau; + allprods.push_back( pmcp ); + SmartRefVector<LHCb::MCVertex> decays = ( *idau )->endVertices(); + SmartRefVector<LHCb::MCVertex>::iterator ivtx; + for ( ivtx = decays.begin(); ivtx != decays.end(); ivtx++ ) { collectProductss( mcpv, *ivtx, allprods ); } } - return nVeloTracks; } -void PrimaryVertexChecker::count_reconstructible_mc_particles( std::vector<MCPVInfo>& mcpvvec ) { +void PrimaryVertexChecker::count_reconstructible_mc_particles( std::vector<MCPVInfo>& mcpvvec, + const LHCb::MCProperty& flags ) const { - const auto trInfo = MCTrackInfo{*get<LHCb::MCProperty>( LHCb::MCPropertyLocation::TrackInfo )}; - for ( auto itinfomc = mcpvvec.begin(); mcpvvec.end() != itinfomc; itinfomc++ ) { - LHCb::MCVertex* avtx = itinfomc->pMCPV; + const MCTrackInfo trInfo = {flags}; + for ( auto& itinfomc : mcpvvec ) { + LHCb::MCVertex* avtx = itinfomc.pMCPV; std::vector<LHCb::MCParticle*> mcPartInMCPV; SmartRefVector<LHCb::MCParticle> parts = avtx->products(); std::vector<LHCb::MCParticle*> allproducts; collectProductss( avtx, avtx, allproducts ); - for ( LHCb::MCParticle* pmcp : allproducts ) { - if ( pmcp->particleID().threeCharge() != 0 && ( !m_requireVelo || trInfo.hasVelo( pmcp ) ) ) { - - if ( pmcp->particleID().hasBottom() ) { itinfomc->decayBeauty = 1.0; } - if ( pmcp->particleID().hasCharm() ) { itinfomc->decayCharm = 1.0; } - if ( pmcp->particleID().threeCharge() != 0 && ( !m_requireVelo || trInfo.hasVelo( pmcp ) ) ) { - - double dv2 = ( avtx->position() - pmcp->originVertex()->position() ).Mag2(); - if ( dv2 < 0.0000001 && pmcp->p() > 100. * Gaudi::Units::MeV ) { mcPartInMCPV.push_back( pmcp ); } - } + for ( const auto pmcp : allproducts ) { + if ( pmcp->particleID().isMeson() || pmcp->particleID().isBaryon() ) { + if ( pmcp->particleID().hasBottom() ) { itinfomc.decayBeauty = true; } + if ( pmcp->particleID().hasCharm() ) { itinfomc.decayCharm = true; } + if ( pmcp->particleID().hasStrange() ) { itinfomc.decayStrange = true; } } - itinfomc->nRecTracks = mcPartInMCPV.size(); - } - } -} -int PrimaryVertexChecker::check_mother_particle( LHCb::RecVertex* pv, std::string trackLoc ) { - - int motherPart = 0; - const auto& tracksIn = pv->tracks(); - - LinkedTo<LHCb::MCParticle, LHCb::Track> link( eventSvc(), msgSvc(), trackLoc ); - for ( const LHCb::Track* ptr : tracksIn ) { - - LHCb::MCParticle* part = nullptr; - part = link.first( ptr->key() ); - if ( part != nullptr ) { - const LHCb::MCParticle* part1 = part; - int i = 0; - while ( part1 != nullptr ) { - // std::cout<<"X"<<std::endl; - const LHCb::MCParticle* mother = part1->mother(); - if ( mother != nullptr ) { - i = i + 1; - part1 = mother; - /// std::cout<<i<<": "<<part1->particleID().pid()<<std::endl; - } else { - break; - } + if ( ( pmcp->particleID().isMeson() || pmcp->particleID().isBaryon() ) && + ( !m_requireVelo || trInfo.hasVelo( pmcp ) ) ) { + double dv2 = ( avtx->position() - pmcp->originVertex()->position() ).Mag2(); + if ( dv2 < 0.0000001 && pmcp->p() > 100. * Gaudi::Units::MeV ) { mcPartInMCPV.push_back( pmcp ); } } - - if ( part1->particleID().hasBottom() == true ) { motherPart = 2.0; } - if ( part1->particleID().hasCharm() == true ) { motherPart = 1.0; } + itinfomc.nRecTracks = mcPartInMCPV.size(); } } - return motherPart; } -bool PrimaryVertexChecker::getInputTracks( std::vector<LHCb::Track*>& vecOfTracks, std::string& trackLoc ) { - - std::string tracksName = m_inputTracksName; - - if ( m_inputTracksName == "Offline" ) { - tracksName = LHCb::TrackLocation::Default; - } else if ( m_inputTracksName == "3D" ) { - tracksName = LHCb::TrackLocation::Velo; - } else if ( m_inputTracksName == "2D" ) { - tracksName = LHCb::TrackLocation::RZVelo; - } - - if ( tracksName == "none" ) { - debug() << " Tracks not specified " << tracksName << endmsg; +bool PrimaryVertexChecker::checkCondition( const MCPVInfo& MCPV, const recoAs& n ) const { + switch ( n ) { + case recoAs::all: + return true; + case recoAs::isolated: + return ( MCPV.distToClosestMCPV > m_dzIsolated ); + case recoAs::close: + return ( MCPV.distToClosestMCPV <= m_dzIsolated ); + case recoAs::ntracks_low: + return ( MCPV.nRecTracks >= m_nTracksToBeRecble && MCPV.nRecTracks < m_nTracksToPrint ); + case recoAs::ntracks_high: + return ( MCPV.nRecTracks >= m_nTracksToPrint ); + case recoAs::z_low: + return ( MCPV.pMCPV->position().z() < -m_zToPrint ); + case recoAs::z_middle: + return ( MCPV.pMCPV->position().z() >= -m_zToPrint && MCPV.pMCPV->position().z() < m_zToPrint ); + case recoAs::z_high: + return ( MCPV.pMCPV->position().z() >= m_zToPrint ); + case recoAs::beauty: + return ( MCPV.decayBeauty ); + case recoAs::charm: + return ( MCPV.decayCharm ); + case recoAs::strange: + return ( MCPV.decayStrange ); + case recoAs::other: + return ( !( MCPV.decayBeauty ) && !( MCPV.decayCharm ) && !( MCPV.decayStrange ) ); + default: return false; } - - trackLoc = tracksName; - - LHCb::Tracks* usedTracks; - - usedTracks = getOrCreate<LHCb::Tracks, LHCb::Tracks>( tracksName ); - if ( usedTracks->size() == 0 ) return false; - - std::vector<LHCb::Track*>::const_iterator itT; - for ( itT = usedTracks->begin(); usedTracks->end() != itT; itT++ ) { - LHCb::Track* ptr = ( *itT ); - vecOfTracks.push_back( ptr ); - } - - return true; } -bool PrimaryVertexChecker::getInputVertices( std::vector<LHCb::RecVertex*>& vecOfVertices ) { - - std::string verticesName = m_inputVerticesName; - - if ( m_inputVerticesName == "Offline" ) { - verticesName = LHCb::RecVertexLocation::Primary; - } else if ( m_inputVerticesName == "3D" ) { - verticesName = LHCb::RecVertexLocation::Velo3D; - } else if ( m_inputVerticesName == "2D" ) { - verticesName = LHCb::RecVertexLocation::Velo2D; - } - - if ( verticesName == "none" ) { - debug() << " Vertices not specified " << verticesName << endmsg; - return false; +std::string PrimaryVertexChecker::toString( const recoAs& n ) const { + switch ( n ) { + case recoAs::all: + return format( "0%d all", int( recoAs::all ) ); + case recoAs::isolated: + return format( "0%d isolated", int( recoAs::isolated ) ); + case recoAs::close: + return format( "0%d close", int( recoAs::close ) ); + case recoAs::ntracks_low: + return format( "0%d ntracks<%d", int( recoAs::ntracks_low ), int( m_nTracksToPrint ) ); + case recoAs::ntracks_high: + return format( "0%d ntracks>=%d", int( recoAs::ntracks_high ), int( m_nTracksToPrint ) ); + case recoAs::z_low: + return format( "0%d z<%2.1f", int( recoAs::z_low ), -m_zToPrint ); + case recoAs::z_middle: + return format( "0%d z in (%2.1f, %2.1f)", int( recoAs::z_middle ), -m_zToPrint, +m_zToPrint ); + case recoAs::z_high: + return format( "0%d z >=%2.1f", int( recoAs::z_high ), +m_zToPrint ); + case recoAs::beauty: + return format( "0%d decayBeauty", int( recoAs::beauty ) ); + case recoAs::charm: + return format( "0%d decayCharm", int( recoAs::charm ) ); + case recoAs::strange: + return format( "%d decayStrange", int( recoAs::strange ) ); + case recoAs::other: + return format( "%d other", int( recoAs::other ) ); + case recoAs::first: + return format( "%d 1MCPV", int( recoAs::first ) ); + case recoAs::second: + return format( "%d 2MCPV", int( recoAs::second ) ); + case recoAs::third: + return format( "%d 3MCPV", int( recoAs::third ) ); + case recoAs::fourth: + return format( "%d 4MCPV", int( recoAs::fourth ) ); + case recoAs::fifth: + return format( "%d 5MCPV", int( recoAs::fifth ) ); + default: + return "not defined"; } - - LHCb::RecVertices* recoVertices; - - recoVertices = getOrCreate<LHCb::RecVertices, LHCb::RecVertices>( verticesName ); - - std::vector<LHCb::RecVertex*>::const_iterator itVer; - for ( itVer = recoVertices->begin(); recoVertices->end() != itVer; itVer++ ) { - LHCb::RecVertex* ppv = ( *itVer ); - vecOfVertices.push_back( ppv ); - } - - return true; } diff --git a/Tr/PatChecker/src/PrimaryVertexChecker.h b/Tr/PatChecker/src/PrimaryVertexChecker.h deleted file mode 100644 index b5d62bd502903e1df6e3fd205b15253964dd814b..0000000000000000000000000000000000000000 --- a/Tr/PatChecker/src/PrimaryVertexChecker.h +++ /dev/null @@ -1,111 +0,0 @@ -/*****************************************************************************\ -* (c) Copyright 2000-2018 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. * -\*****************************************************************************/ -#ifndef PRIMARYVERTEXCHECKER_H -#define PRIMARYVERTEXCHECKER_H 1 - -#include "Event/MCVertex.h" -#include "Event/RecVertex.h" -#include "Event/Track.h" -#include "MCInterfaces/IForcedBDecayTool.h" - -#include "GaudiAlg/GaudiTupleAlg.h" - -typedef struct { - LHCb::MCVertex* pMCPV; // pointer to MC PV - int nRecTracks; // number of reconstructed tracks from this MCPV - int nRecBackTracks; // number of reconstructed backward tracks - int indexRecPVInfo; // index to reconstructed PVInfo (-1 if not reco) - int nCorrectTracks; // correct tracks belonging to reconstructed PV - int multClosestMCPV; // multiplicity of closest reconstructable MCPV - double distToClosestMCPV; // distance to closest reconstructible MCPV - int decayCharm; // type of mother particle - int decayBeauty; - std::vector<LHCb::MCParticle*> m_mcPartInMCPV; - std::vector<LHCb::Track*> m_recTracksInMCPV; -} MCPVInfo; - -typedef struct { -public: - int nTracks; // number of tracks - int nVeloTracks; // number of velo tracks in a vertex - int nLongTracks; - double minTrackRD; // - double maxTrackRD; // - double chi2; - double nDoF; - double d0; - double d0nTr; - double chi2nTr; - double mind0; - double maxd0; - int mother; - Gaudi::XYZPoint position; // position - Gaudi::XYZPoint positionSigma; // position sigmas - int indexMCPVInfo; // index to MCPVInfo - LHCb::RecVertex* pRECPV; // pointer to REC PV -} RecPVInfo; - -class PrimaryVertexChecker : public GaudiTupleAlg { -public: - /// Standard constructor - PrimaryVertexChecker( const std::string& name, ISvcLocator* pSvcLocator ); - - StatusCode initialize() override; ///< Algorithm initialization - StatusCode execute() override; ///< Algorithm execution - StatusCode finalize() override; ///< Algorithm finalization - -protected: - bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ); } - -private: - int m_nTracksToBeRecble; - bool m_produceHistogram; - bool m_produceNtuple; - bool m_requireVelo; - double m_dzIsolated; - bool m_matchByTracks; - std::string m_inputTracksName; - std::string m_inputVerticesName; - - int m_nevt; - int m_nMCPV, m_nRecMCPV; - int m_nMCPV_isol, m_nRecMCPV_isol; - int m_nMCPV_close, m_nRecMCPV_close; - int m_nFalsePV, m_nFalsePV_real; - int m_nMCPV_1mult, m_nRecMCPV_1mult; - int m_nMCPV_isol_1mult, m_nRecMCPV_isol_1mult; - int m_nMCPV_close_1mult, m_nRecMCPV_close_1mult; - int m_nRecMCPV_wrong_1mult; - int m_nMCPV_2mult, m_nRecMCPV_2mult; - int m_nMCPV_isol_2mult, m_nRecMCPV_isol_2mult; - int m_nMCPV_close_2mult, m_nRecMCPV_close_2mult; - int m_nRecMCPV_wrong_2mult; - int m_nMCPV_3mult, m_nRecMCPV_3mult; - int m_nMCPV_isol_3mult, m_nRecMCPV_isol_3mult; - int m_nMCPV_close_3mult, m_nRecMCPV_close_3mult; - int m_nRecMCPV_wrong_3mult; - int m_nBFalse, m_nRecBFalse; - - std::vector<MCPVInfo>::iterator closestMCPV( std::vector<MCPVInfo>& rblemcpv, std::vector<MCPVInfo>::iterator& itmc ); - - void collectProductss( LHCb::MCVertex* mcpv, LHCb::MCVertex* mcvtx, std::vector<LHCb::MCParticle*>& allprods ); - void printRat( std::string mes, int a, int b ); - void match_mc_vertex_by_tracks( int ipv, std::vector<RecPVInfo>& rinfo, std::vector<MCPVInfo>& mcpvvec ); - void match_mc_vertex_by_distance( int ipv, std::vector<RecPVInfo>& rinfo, std::vector<MCPVInfo>& mcpvvec ); - void count_reconstructed_tracks( std::vector<MCPVInfo>& mcpvvec, std::vector<LHCb::Track*>& vecOfTracks, - std::string trackLoc ); - int count_velo_tracks( LHCb::RecVertex* RecVtx ); - void count_reconstructible_mc_particles( std::vector<MCPVInfo>& mcpvvec ); - bool getInputTracks( std::vector<LHCb::Track*>& vecOfTracks, std::string& trackLoc ); - bool getInputVertices( std::vector<LHCb::RecVertex*>& vecOfVertices ); - int check_mother_particle( LHCb::RecVertex* RecVtx, std::string trackLoc ); -}; -#endif // PRIMARYVERTEXCHECKER_H