diff --git a/Tr/PatPV/CMakeLists.txt b/Tr/PatPV/CMakeLists.txt index 470918c58b61ef0e34d13732a931e3424cdfa000..c307832d7529daea2622394fa0f9069c6f553a91 100644 --- a/Tr/PatPV/CMakeLists.txt +++ b/Tr/PatPV/CMakeLists.txt @@ -31,7 +31,6 @@ gaudi_add_module(PatPV src/SimplePVSeedTool.cpp src/TrackBeamLineVertexFinder.cpp src/TrackBeamLineVertexFinderSoA.cpp - src/TrackPVFinderUtils.cpp src/TrackUnbiasedPVFinderSoA.cpp src/PatPVLeftRightFinderMerger.cpp src/PatPVConverters.cpp diff --git a/Tr/PatPV/src/AdaptivePV3DFitter.cpp b/Tr/PatPV/src/AdaptivePV3DFitter.cpp index f8354ed1b309b1809b1d2fc1bcf09d54c4941b36..27bc56bcf749c686f8b69b313059e7e4ddb64e8c 100644 --- a/Tr/PatPV/src/AdaptivePV3DFitter.cpp +++ b/Tr/PatPV/src/AdaptivePV3DFitter.cpp @@ -23,8 +23,9 @@ public: using extends::extends; // Fitting StatusCode fitVertex( const Gaudi::XYZPoint& seedPoint, LHCb::span tracks, - LHCb::RecVertex& vtx, std::vector& tracks2remove, - IGeometryInfo const& geometry ) const override; + LHCb::RecVertex& vtx, std::vector& tracks2remove ) const override; + + using IPVFitter::fitVertex; private: Gaudi::Property m_minTr{ this, "MinTracks", 4, "Minimum number of tracks to make a vertex" }; @@ -136,7 +137,7 @@ namespace { //============================================================================= StatusCode AdaptivePV3DFitter::fitVertex( const Gaudi::XYZPoint& seedPoint, LHCb::span rTracks, LHCb::RecVertex& vtx, - std::vector& tracks2remove, IGeometryInfo const& ) const { + std::vector& tracks2remove ) const { tracks2remove.clear(); // position at which derivatives are evaluated diff --git a/Tr/PatPV/src/LSAdaptPV3DFitter.cpp b/Tr/PatPV/src/LSAdaptPV3DFitter.cpp index 747451bc1b7a573500078fa322c39562ac0b40c6..8b78271c19f98c6cba69a193af37980564177f26 100644 --- a/Tr/PatPV/src/LSAdaptPV3DFitter.cpp +++ b/Tr/PatPV/src/LSAdaptPV3DFitter.cpp @@ -24,8 +24,9 @@ public: using extends::extends; // Fitting StatusCode fitVertex( const Gaudi::XYZPoint& seedPoint, LHCb::span tracks, - LHCb::RecVertex& vtx, std::vector& tracks2remove, - IGeometryInfo const& geometry ) const override; + LHCb::RecVertex& vtx, std::vector& tracks2remove ) const override; + + using IPVFitter::fitVertex; private: Gaudi::Property m_maxDeltaZ{ this, "maxDeltaZ", 0.0005 * Gaudi::Units::mm, "Fit convergence condition" }; @@ -96,8 +97,7 @@ DECLARE_COMPONENT( LSAdaptPV3DFitter ) // Least square adaptive fitting method //============================================================================= StatusCode LSAdaptPV3DFitter::fitVertex( const Gaudi::XYZPoint& seedPoint, LHCb::span rTracks, - LHCb::RecVertex& vtx, std::vector& tracks2remove, - IGeometryInfo const& ) const { + LHCb::RecVertex& vtx, std::vector& tracks2remove ) const { if ( msgLevel( MSG::DEBUG ) ) debug() << "================Test==================" << endmsg; tracks2remove.clear(); diff --git a/Tr/PatPV/src/LSAdaptPVFitter.cpp b/Tr/PatPV/src/LSAdaptPVFitter.cpp index fab67f13c2519f0eeda4ab8543b14bd8547469b8..f28eb6e328525c2ec630c4d00ec943e6c8326a4d 100644 --- a/Tr/PatPV/src/LSAdaptPVFitter.cpp +++ b/Tr/PatPV/src/LSAdaptPVFitter.cpp @@ -26,6 +26,12 @@ public: LHCb::RecVertex& vtx, std::vector& tracks2remove, IGeometryInfo const& geometry ) const override; + // Fitting + StatusCode fitVertex( const Gaudi::XYZPoint&, LHCb::span, LHCb::RecVertex&, + std::vector& ) const override { + throw GaudiException( __PRETTY_FUNCTION__, "not implemented", StatusCode::FAILURE ); + } + private: Gaudi::Property m_minTr{ this, "MinTracks", 5, [this]( auto& ) { diff --git a/Tr/PatPV/src/ParticleUnbiasedPVAdder.cpp b/Tr/PatPV/src/ParticleUnbiasedPVAdder.cpp index 2b671364fa7a8b8dfec2db1ac1f8481e91ae1db1..fee34fb7632bc743348d867f25dcd8f47df9ea3f 100644 --- a/Tr/PatPV/src/ParticleUnbiasedPVAdder.cpp +++ b/Tr/PatPV/src/ParticleUnbiasedPVAdder.cpp @@ -12,8 +12,7 @@ #include "GaudiKernel/KeyedContainer.h" #include "LHCbAlgs/Transformer.h" #include "TrackKernel/PrimaryVertexUtils.h" - -#include "TrackPVFinderUtils.h" +#include "TrackKernel/TrackPVFinderUtils.h" using namespace LHCb::TrackPVFinder; using Particles = LHCb::Particles; diff --git a/Tr/PatPV/src/PatPVConverters.cpp b/Tr/PatPV/src/PatPVConverters.cpp index e576e71adf8606ad1482a1f9a9f4262df136bd46..97aa803010b8e272b048b7cbdd444f230d885ec5 100644 --- a/Tr/PatPV/src/PatPVConverters.cpp +++ b/Tr/PatPV/src/PatPVConverters.cpp @@ -14,7 +14,7 @@ #include "Event/RecVertex_v2.h" #include "Event/Track.h" #include "LHCbAlgs/Transformer.h" -#include "TrackPVFinderUtils.h" +#include "TrackKernel/TrackPVFinderUtils.h" // boost includes #include @@ -196,7 +196,8 @@ public: // create the output container LHCb::Event::PV::PrimaryVertexContainer pvcontainer; // fill it - LHCb::TrackPVFinder::populateFromRecVertices( pvcontainer, recvertices ); + LHCb::TrackPVFinder::populateFromRecVertices( + pvcontainer, LHCb::make_span( std::to_address( recvertices.begin() ), recvertices.size() ) ); // in debug mode: call the vertex fit if ( msgLevel() == MSG::DEBUG ) { debug() << "Test: refit and compare result for all vertices" << endmsg; diff --git a/Tr/PatPV/src/RecVertexAsExtendedPrimaryVertexFitter.cpp b/Tr/PatPV/src/RecVertexAsExtendedPrimaryVertexFitter.cpp index 9ddbc355f251933b08dc9f54747484a64daa494d..33be6823dc7d202ec282b89274c6fe1fd8b3a76c 100644 --- a/Tr/PatPV/src/RecVertexAsExtendedPrimaryVertexFitter.cpp +++ b/Tr/PatPV/src/RecVertexAsExtendedPrimaryVertexFitter.cpp @@ -12,7 +12,7 @@ #include "Event/RecVertex.h" #include "GaudiAlg/GaudiTool.h" #include "TrackInterfaces/IPVFitter.h" -#include "TrackPVFinderUtils.h" +#include "TrackKernel/TrackPVFinderUtils.h" using namespace LHCb::Event::PV; @@ -33,9 +33,9 @@ public: // Standard constructor using extends::extends; // Fitting + using IPVFitter::fitVertex; StatusCode fitVertex( const Gaudi::XYZPoint& seedPoint, LHCb::span tracks, - LHCb::RecVertex& recvertex, std::vector& tracks2remove, - IGeometryInfo const& ) const override { + LHCb::RecVertex& recvertex, std::vector& tracks2remove ) const override { // convert the RecVertex to en extendedprimaryvertex diff --git a/Tr/PatPV/src/SimplePVFitter.cpp b/Tr/PatPV/src/SimplePVFitter.cpp index acb2b6fb297e66e0d93408873e083fa217a7d55c..c583d3cf7d533498a2082cc7891a9d608b9b77ee 100644 --- a/Tr/PatPV/src/SimplePVFitter.cpp +++ b/Tr/PatPV/src/SimplePVFitter.cpp @@ -26,6 +26,12 @@ public: LHCb::RecVertex& vtx, std::vector& tracks2remove, IGeometryInfo const& geometry ) const override; + // Fitting + StatusCode fitVertex( const Gaudi::XYZPoint&, LHCb::span, LHCb::RecVertex&, + std::vector& ) const override { + throw GaudiException( __PRETTY_FUNCTION__, "not implemented", StatusCode::FAILURE ); + } + private: Gaudi::Property m_minTr{ this, "MinTracks", 5, [this]( auto& ) { diff --git a/Tr/PatPV/src/TrackBeamLineVertexFinderSoA.cpp b/Tr/PatPV/src/TrackBeamLineVertexFinderSoA.cpp index 7fd13710eadf7f085fe3255dbc1cece175ae04ec..f81782d70d7435600150dd72b8475a8d3bf3c675 100644 --- a/Tr/PatPV/src/TrackBeamLineVertexFinderSoA.cpp +++ b/Tr/PatPV/src/TrackBeamLineVertexFinderSoA.cpp @@ -25,7 +25,7 @@ # include "Timer.h" #endif -#include "TrackPVFinderUtils.h" +#include "TrackKernel/TrackPVFinderUtils.h" // boost includes #include diff --git a/Tr/PatPV/src/TrackUnbiasedPVFinderSoA.cpp b/Tr/PatPV/src/TrackUnbiasedPVFinderSoA.cpp index 85dce148401d4a92c9e45fdc85b209eabe58d2d2..58b48998a184d87ce05841f80073c953eaa144c2 100644 --- a/Tr/PatPV/src/TrackUnbiasedPVFinderSoA.cpp +++ b/Tr/PatPV/src/TrackUnbiasedPVFinderSoA.cpp @@ -27,7 +27,7 @@ # include "TrackKernel/Timer.h" #endif -#include "TrackPVFinderUtils.h" +#include "TrackKernel/TrackPVFinderUtils.h" // std includes #include diff --git a/Tr/TrackFitEvent/src/PrKalmanFitResult.cpp b/Tr/TrackFitEvent/src/PrKalmanFitResult.cpp index b2c52a6d984c1ba59666f3c5f2247f7ada41a05b..35e570da4cf0fb4ef73eb13b4642bd1f65720d29 100644 --- a/Tr/TrackFitEvent/src/PrKalmanFitResult.cpp +++ b/Tr/TrackFitEvent/src/PrKalmanFitResult.cpp @@ -40,20 +40,19 @@ namespace LHCb { } ChiSquare PrKalmanFitResult::chi2Upstream() const { - auto upstream_chi2 = ChiSquare{}; + auto upstream_chi2 = ChiSquare{ 0, -number_of_trackparameters }; for ( auto const& node : fitnodes ) { - if ( node.type() == Pr::Tracks::Fit::Node::Type::UTHit ) { + if ( node.type() == Pr::Tracks::Fit::Node::Type::UTHit || node.type() == Pr::Tracks::Fit::Node::Type::VPHit ) { upstream_chi2 += node.delta_chi2[Pr::Tracks::Fit::Node::backward]; } } - upstream_chi2 += chi2Velo(); return upstream_chi2; } ChiSquare PrKalmanFitResult::chi2Downstream() const { auto down_chi2 = ChiSquare{ 0, -number_of_trackparameters }; for ( auto const& node : fitnodes ) { - if ( node.type() == Pr::Tracks::Fit::Node::Type::FTHit ) { + if ( node.type() == Pr::Tracks::Fit::Node::Type::FTHit || node.type() == Pr::Tracks::Fit::Node::Type::MuonHit ) { down_chi2 += node.delta_chi2[Pr::Tracks::Fit::Node::forward]; } } @@ -63,5 +62,13 @@ namespace LHCb { ChiSquare PrKalmanFitResult::chi2Match() const { return chi2() - chi2Upstream() - chi2Downstream(); } /// TODO PrKalmanFilter does not support muon tracks at the moment - ChiSquare PrKalmanFitResult::chi2Muon() const { return ChiSquare{}; } + ChiSquare PrKalmanFitResult::chi2Muon() const { + auto muon_chi2 = ChiSquare{ 0, -number_of_trackparameters }; + for ( auto const& node : fitnodes ) { + if ( node.type() == Pr::Tracks::Fit::Node::Type::MuonHit ) { + muon_chi2 += node.delta_chi2[Pr::Tracks::Fit::Node::forward]; + } + } + return muon_chi2; + } } // namespace LHCb diff --git a/Tr/TrackInterfaces/include/TrackInterfaces/IPVFitter.h b/Tr/TrackInterfaces/include/TrackInterfaces/IPVFitter.h index 5c25b35724ba038cf84b48267abc8aae730bbafd..cca943769b7c71202695d51204c94cc2e7b45ebf 100644 --- a/Tr/TrackInterfaces/include/TrackInterfaces/IPVFitter.h +++ b/Tr/TrackInterfaces/include/TrackInterfaces/IPVFitter.h @@ -25,15 +25,21 @@ namespace LHCb { struct IPVFitter : extend_interfaces { DeclareInterfaceID( IPVFitter, 4, 0 ); + + virtual StatusCode fitVertex( const Gaudi::XYZPoint& seedPoint, LHCb::span tracks, + LHCb::RecVertex& vtx, std::vector& tracks2remove ) const = 0; + virtual StatusCode fitVertex( const Gaudi::XYZPoint& seedPoint, LHCb::span tracks, LHCb::RecVertex& vtx, std::vector& tracks2remove, - IGeometryInfo const& geometry ) const = 0; + IGeometryInfo const& ) const { + return fitVertex( seedPoint, tracks, vtx, tracks2remove ); + } - std::optional fit( const Gaudi::XYZPoint& seedPoint, LHCb::span tracks, - IGeometryInfo const& geometry ) const { + std::optional fit( const Gaudi::XYZPoint& seedPoint, + LHCb::span tracks ) const { std::vector tracks2remove; LHCb::RecVertex vtx; - if ( fitVertex( seedPoint, tracks, vtx, tracks2remove, geometry ).isSuccess() ) return vtx; + if ( fitVertex( seedPoint, tracks, vtx, tracks2remove ).isSuccess() ) return vtx; return {}; } }; diff --git a/Tr/TrackInterfaces/include/TrackInterfaces/ITrackVertexer.h b/Tr/TrackInterfaces/include/TrackInterfaces/ITrackVertexer.h index 10600dbd975103a51d721080110f383dbe36a43b..9559fc299973b6c526fcd0c547eba767611e6397 100644 --- a/Tr/TrackInterfaces/include/TrackInterfaces/ITrackVertexer.h +++ b/Tr/TrackInterfaces/include/TrackInterfaces/ITrackVertexer.h @@ -23,6 +23,7 @@ namespace LHCb { class TwoProngVertex; class State; class RecVertex; + class VertexBase; } // namespace LHCb /** @@ -34,19 +35,18 @@ struct ITrackVertexer : extend_interfaces { DeclareInterfaceID( ITrackVertexer, 4, 0 ); - virtual std::unique_ptr fit( const LHCb::State& stateA, const LHCb::State& stateB, - IGeometryInfo const& geometry ) const = 0; - virtual std::unique_ptr fit( LHCb::span states, - IGeometryInfo const& geometry ) const = 0; + virtual std::unique_ptr fit( const LHCb::State& stateA, const LHCb::State& stateB ) const = 0; + virtual std::unique_ptr fit( LHCb::span states ) const = 0; virtual std::unique_ptr fit( LHCb::span tracks, - IGeometryInfo const& geometry ) const = 0; - virtual bool computeDecayLength( const LHCb::TwoProngVertex& vertex, const LHCb::RecVertex& pv, double& chi2, - double& decaylength, double& decaylengtherr ) const = 0; + IGeometryInfo const& geometry ) const = 0; + virtual bool computeDecayLength( const LHCb::TwoProngVertex& vertex, const LHCb::VertexBase& pv, double& chi2, + double& decaylength, double& decaylengtherr ) const = 0; /// Return the ip chi2 for a track (uses stateprovider, not good for /// HLT: better call routine below with track->firstState()) - virtual double ipchi2( const LHCb::Track& track, const LHCb::RecVertex& pv, IGeometryInfo const& geometry ) const = 0; + virtual double ipchi2( const LHCb::Track& track, const LHCb::VertexBase& pv, + IGeometryInfo const& geometry ) const = 0; /// Return the ip chi2 for a track state - virtual double ipchi2( const LHCb::State& state, const LHCb::RecVertex& pv, IGeometryInfo const& geometry ) const = 0; + virtual double ipchi2( const LHCb::State& state, const LHCb::VertexBase& pv ) const = 0; }; diff --git a/Tr/TrackKernel/CMakeLists.txt b/Tr/TrackKernel/CMakeLists.txt index 8795a61a8d13b84fb7b33f686c5b534307a96810..10ac7f3ef7b69e6e79345210f788350b99747eda 100644 --- a/Tr/TrackKernel/CMakeLists.txt +++ b/Tr/TrackKernel/CMakeLists.txt @@ -24,6 +24,7 @@ gaudi_add_library(TrackKernel src/TrackStateVertex.cpp src/TrackTraj.cpp src/TrajVertex.cpp + src/TrackPVFinderUtils.cpp LINK PUBLIC Gaudi::GaudiKernel diff --git a/Tr/PatPV/src/TrackPVFinderUtils.h b/Tr/TrackKernel/include/TrackKernel/TrackPVFinderUtils.h similarity index 99% rename from Tr/PatPV/src/TrackPVFinderUtils.h rename to Tr/TrackKernel/include/TrackKernel/TrackPVFinderUtils.h index 4fe4074a7b1df3f440b917402f14007ff5c857bf..1ab5b8b7dc824620be041e1fbb89c6abd5699a40 100644 --- a/Tr/PatPV/src/TrackPVFinderUtils.h +++ b/Tr/TrackKernel/include/TrackKernel/TrackPVFinderUtils.h @@ -49,7 +49,7 @@ namespace LHCb::TrackPVFinder { void populateFromVeloTracks( PrimaryVertexContainer& pvdata ); /// Fill PrimaryVertex data object from RecVertices and associated tracks - void populateFromRecVertices( PrimaryVertexContainer& pvcontainer, const RecVertex::Range& recvertices ); + void populateFromRecVertices( PrimaryVertexContainer& pvcontainer, LHCb::span recvertices ); /// Reassigns the labels from the PV tracks to the vertices. This is needed if a selection has been made to the /// vertices. diff --git a/Tr/PatPV/src/TrackPVFinderUtils.cpp b/Tr/TrackKernel/src/TrackPVFinderUtils.cpp similarity index 98% rename from Tr/PatPV/src/TrackPVFinderUtils.cpp rename to Tr/TrackKernel/src/TrackPVFinderUtils.cpp index 220b67eed890c14377a02148a4149068abc646a0..11ee2b9277ea39fb9611842df439b9128546368f 100644 --- a/Tr/PatPV/src/TrackPVFinderUtils.cpp +++ b/Tr/TrackKernel/src/TrackPVFinderUtils.cpp @@ -8,7 +8,7 @@ * granted to it by virtue of its status as an Intergovernmental Organization * * or submit itself to any jurisdiction. * \*****************************************************************************/ -#include "TrackPVFinderUtils.h" +#include "TrackKernel/TrackPVFinderUtils.h" namespace LHCb::TrackPVFinder { @@ -50,7 +50,7 @@ namespace LHCb::TrackPVFinder { populateFromVeloTracks( pvdata ); } - void populateFromRecVertices( PrimaryVertexContainer& pvcontainer, const RecVertex::Range& recvertices ) { + void populateFromRecVertices( PrimaryVertexContainer& pvcontainer, LHCb::span recvertices ) { auto& vertices = pvcontainer.vertices; vertices.clear(); vertices.resize( recvertices.size() ); diff --git a/Tr/TrackMonitors/CMakeLists.txt b/Tr/TrackMonitors/CMakeLists.txt index adc5f4e7a5c1fbc39e48d021021c45ecd6e5e655..abb9c9a236aff2c8d0c79149513d4f456a27ad05 100644 --- a/Tr/TrackMonitors/CMakeLists.txt +++ b/Tr/TrackMonitors/CMakeLists.txt @@ -40,6 +40,7 @@ gaudi_add_module(TrackMonitors src/VPHitEfficiencyMonitor.cpp src/VertexCompare.cpp src/BeamSpotMonitor.cpp + src/TrackMonitorNT.cpp LINK AIDA::aida Gaudi::GaudiAlgLib diff --git a/Tr/TrackMonitors/src/TrackFitMatchMonitor.cpp b/Tr/TrackMonitors/src/TrackFitMatchMonitor.cpp index ae8c8b8a3bf450bb03ceed22a10ded931f4504e0..86e7c260e5da14804fe5a69eed8af38dfa5125d9 100644 --- a/Tr/TrackMonitors/src/TrackFitMatchMonitor.cpp +++ b/Tr/TrackMonitors/src/TrackFitMatchMonitor.cpp @@ -43,7 +43,7 @@ public: private: // Trackers that have plots produced (may need to change in order to get the relevant Trackers from some centralized // variable rather than hard-coding them) - enum struct Trackers { VeloUT = 0, TUT, VeloT }; + enum struct Trackers { VeloUT = 0, TUT, VeloT, MuonT }; template void plotDelta( Trackers tracker, const TNode& node, bool upstream ) const; @@ -206,9 +206,10 @@ private: , dty_pull_ty_profile{ owner, Tracker + "/dty pull vs ty", Tracker + " dty pull vs ty", { 20, -0.25, 0.25 } } {} }; // map to associate tracking sub-detector to the relevant struct of histograms - std::array, 3> m_histograms = [&] { + std::array, 4> m_histograms = [&] { auto create = [&]( Trackers t ) { return std::make_unique( this, toString( t ) ); }; - return std::array{ create( Trackers::VeloUT ), create( Trackers::TUT ), create( Trackers::VeloT ) }; + return std::array{ create( Trackers::VeloUT ), create( Trackers::TUT ), create( Trackers::VeloT ), + create( Trackers::MuonT ) }; }(); /** Friend function for converting enum values to strings **/ @@ -220,6 +221,8 @@ private: return "T-UT"; case Trackers::VeloT: return "Velo-T"; + case Trackers::MuonT: + return "Muon-T"; } throw GaudiException( "Unknown Tracker enum value", __func__, StatusCode::FAILURE ); } @@ -381,7 +384,8 @@ template void TrackFitMatchMonitor::fill( const TFitResult& fitResult, const LHCb::Track& track ) const { if ( nodes( fitResult ).size() <= 0 ) return; - const typename TFitResult::NodeType *lastVelo( nullptr ), *firstUT( nullptr ), *lastUT( nullptr ), *firstT( nullptr ); + const typename TFitResult::NodeType *lastVelo( nullptr ), *firstUT( nullptr ), *lastUT( nullptr ), *firstT( nullptr ), + *firstMuon{ nullptr }; // The appropriate function node( TFitResult* ) will be found using ADL. // For TrackFitResult or KalmanFitResult it returns Range of LHCb::FitNodes, see TrackFitResult.h and FitNode.h // For PrKalmanFitResult it returns std::span of PrFitNodes, see PrKalmanFitResult.h and PrFitNode.h @@ -396,6 +400,8 @@ void TrackFitMatchMonitor::fill( const TFitResult& fitResult, const LHCb::Track& if ( !lastUT || lastUT->z() < node.z() ) lastUT = &node; } else if ( node.isFT() ) { if ( !firstT || firstT->z() > node.z() ) firstT = &node; + } else if ( node.isMuon() ) { + if ( !firstMuon || firstMuon->z() > node.z() ) firstMuon = &node; } } if ( lastVelo ) { @@ -406,6 +412,7 @@ void TrackFitMatchMonitor::fill( const TFitResult& fitResult, const LHCb::Track& plotDelta( Trackers::VeloT, *firstT, true ); } } + if ( firstMuon && firstT ) plotDelta( Trackers::MuonT, *firstMuon, true ); // inspired by the problems we see in the field. see also UT field study LHCb::HitPattern hitpattern{ track.lhcbIDs() }; diff --git a/Tr/TrackMonitors/src/TrackMonitor.cpp b/Tr/TrackMonitors/src/TrackMonitor.cpp index d603348b861b0d3cdfb31cf26dfc3afaa0d4a2bc..52020fb9a586d9f00957f2f103579805454ab5cc 100644 --- a/Tr/TrackMonitors/src/TrackMonitor.cpp +++ b/Tr/TrackMonitors/src/TrackMonitor.cpp @@ -60,7 +60,7 @@ namespace { const std::vector HitTypeName{ "VPX", "VPY", "VP2D", "UT", "FT", "Muon" }; // TODO: not sure if these values still make sense for Run3? - constexpr auto HitTypeMaxRes = std::array{ 0.1, 0.1, 0.1, 0.5, 1.0, 10.0 }; + constexpr auto HitTypeMaxRes = std::array{ 0.1, 0.1, 0.1, 0.5, 1.0, 200.0 }; template inline HitType hittypemap( const TNode& node ) { if ( node.isMuon() ) @@ -346,10 +346,6 @@ StatusCode TrackMonitor::initialize() { } void TrackMonitor::initializeHistogramMap() { - // range for residuals for different hittypes - const std::vector> range{ - { { "8", 0., 100. }, { "9", 0., 50. }, { "10", 0., 100. }, { "11", 0., 50. }, { "12", 0., 100. } } }; - if ( m_splitByType ) { // make seperate histogrammer for each requested type for ( const auto& t : m_typesToMonitor ) { diff --git a/Tr/TrackMonitors/src/TrackMonitorNT.cpp b/Tr/TrackMonitors/src/TrackMonitorNT.cpp new file mode 100644 index 0000000000000000000000000000000000000000..54d0f4d0421bd039a47c498cbfac9e981be93193 --- /dev/null +++ b/Tr/TrackMonitors/src/TrackMonitorNT.cpp @@ -0,0 +1,221 @@ +/*****************************************************************************\ +* (c) Copyright 2000-2019 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/KalmanFitResult.h" +#include "Event/PrKalmanFitResult.h" +#include "Event/PrimaryVertices.h" +#include "Event/State.h" +#include "Event/Track.h" +#include "GaudiAlg/GaudiTupleAlg.h" +#include "Kernel/HitPattern.h" +#include "Kernel/LHCbID.h" +#include "LHCbAlgs/Consumer.h" +#include "TrackKernel/PrimaryVertexUtils.h" + +/* + * Class for track monitoring + * @date 7-4-2025 + */ + +namespace { + + template + void fillT( Tuples::Tuple& theTuple, const char* column, const S& value ) { + theTuple->column( column, T( value ) ).ignore( /* AUTOMATICALLY ADDED FOR gaudi/Gaudi!763 */ ); + } + + // FIXME: we could template this on the PV container, if needed. + using Vertices = LHCb::Event::PV::PrimaryVertexContainer; + + void fillChi2PerDoF( Tuples::Tuple& tuple, const char* column, const LHCb::ChiSquare& chi2 ) { + fillT( tuple, column, ( chi2.nDoF() > 0 ? chi2.chi2PerDoF() : 0.0 ) ); + } + + template + void fillFitResult( Tuples::Tuple& tuple, const T& fitResult ) { + fillT( tuple, "nActiveMeasurements", fitResult.nActiveMeasurements() ); + fillT( tuple, "nFitIter", fitResult.nIter() ); + fillT( tuple, "pScatter", fitResult.pScatter() ); + const auto fitnodes = nodes( fitResult ); + fillT( + tuple, "noutlier", + std::count_if( std::begin( fitnodes ), std::end( fitnodes ), []( const auto& n ) { return n.isOutlier(); } ) ); + fillChi2PerDoF( tuple, "chi2dofVelo", fitResult.chi2Velo() ); + fillChi2PerDoF( tuple, "chi2dofDownstream", fitResult.chi2Downstream() ); + fillChi2PerDoF( tuple, "chi2dofMatch", fitResult.chi2Match() ); + fillChi2PerDoF( tuple, "chi2dofMuon", fitResult.chi2Muon() ); + } + + void fillTuple( Tuples::Tuple tuple, const LHCb::Track::Range& tracks, + const LHCb::Event::PV::PrimaryVertexContainer& pvs ) { + + int itrack = 0; + + // compute some extra info per pv. rather expensive so do only once + struct PVExtraInfo { + unsigned nTracks{ 0 }; + float totalWeight{ 0 }; + float totalChi2{ 0 }; + }; + + std::vector pvinfo; + const auto pvtracks = pvs.tracks.scalar(); + for ( const auto& vertex : pvs ) { + PVExtraInfo info; + for ( int trkindex = vertex.begin(); trkindex != vertex.end(); ++trkindex ) { + const auto pvtrack = pvtracks[trkindex]; + const auto weight = pvtrack.weight().cast(); + if ( weight > 0 ) { + info.totalWeight += weight; + info.totalChi2 += pvtrack.ipchi2().cast(); + ++info.nTracks; + } + } + pvinfo.emplace_back( info ); + } + + for ( const LHCb::Track* track : tracks ) { + // keep track of track multiplicity + fillT( tuple, "itrack", itrack++ ); + fillT( tuple, "ntrack", tracks.size() ); + fillT( tuple, "nPV", pvs.size() ); + + // fill lot's of info from the track + const auto& state = track->firstState(); + const auto trkx = state.x(); + const auto trky = state.y(); + const auto trkz = state.z(); + const auto trktx = state.tx(); + const auto trkty = state.ty(); + fillT( tuple, "x", trkx ); + fillT( tuple, "y", trky ); + fillT( tuple, "z", trkz ); + fillT( tuple, "tx", trktx ); + fillT( tuple, "ty", trkty ); + fillT( tuple, "qOverP", state.qOverP() ); + fillT( tuple, "eta", track->pseudoRapidity() ); + fillT( tuple, "phi", track->phi() ); + fillT( tuple, "chi2", track->chi2() ); + fillT( tuple, "ndof", track->nDoF() ); + fillT( tuple, "ghostprob", track->ghostProbability() ); + fillT( tuple, "type", track->type() ); + fillT( tuple, "backward", track->isVeloBackward() ); + const auto& ids = track->lhcbIDs(); + fillT( tuple, "nUThits", + std::count_if( ids.begin(), ids.end(), []( const LHCb::LHCbID& id ) { return id.isUT(); } ) ); + fillT( tuple, "nVPhits", + std::count_if( ids.begin(), ids.end(), []( const LHCb::LHCbID& id ) { return id.isVP(); } ) ); + fillT( tuple, "nFThits", + std::count_if( ids.begin(), ids.end(), []( const LHCb::LHCbID& id ) { return id.isFT(); } ) ); + + // If the fit result it available, fill a bit more + fillT( tuple, "trackWasFitted", track->fitResult() != nullptr ); + if ( track->fitResult() ) { + auto prFitResult = dynamic_cast( track->fitResult() ); + if ( prFitResult ) + fillFitResult( tuple, *prFitResult ); + else { + auto masterFitResult = dynamic_cast( track->fitResult() ); + if ( masterFitResult ) fillFitResult( tuple, *masterFitResult ); + } + } + + // Fill some info from the hitpattern + LHCb::HitPattern hitpattern{ ids }; + fillT( tuple, "numVeloStations", hitpattern.numVeloStations() ); + fillT( tuple, "numVeloStationsOverlap", hitpattern.numVeloStationsOverlap() ); + fillT( tuple, "numVeloHoles", hitpattern.numVeloHoles() ); + fillT( tuple, "numUTLayers", hitpattern.numUT() ); + fillT( tuple, "numFTLayers", hitpattern.numFT() ); + fillT( tuple, "numFTHoles", hitpattern.numFTHoles() ); + fillT( tuple, "numVeloA", hitpattern.numVeloA() ); + fillT( tuple, "numVeloC", hitpattern.numVeloC() ); + + // get the closest PV and unbias. there are various ways to do this. + // the most efficient one is to first find the closest PV, then unbias is. + const auto pvindex = LHCb::bestPVIndex( pvs, trkx, trky, trkz, trktx, trkty ); + if ( pvindex != LHCb::Event::PV::PVIndexInvalid ) { + // extract the velo segment ID needed for unbiasing + auto veloid = LHCb::Event::PV::uniqueVeloSegmentID( ids ); + // put this into an array, and get the unbiased PV + std::array vetoedvelotracks{ veloid }; + const auto pv = LHCb::Event::PV::unbiasedVertex( pvs, pvindex, vetoedvelotracks ); + // fill some info on the PV (without the track) + const auto pvx = pv.position().x(); + const auto pvy = pv.position().y(); + const auto pvz = pv.position().z(); + fillT( tuple, "pvx", pvx ); + fillT( tuple, "pvy", pvy ); + fillT( tuple, "pvz", pvz ); + fillT( tuple, "pvntracks", pvinfo[pvindex].nTracks ); + fillT( tuple, "pvtotalweight", pvinfo[pvindex].totalWeight ); + fillT( tuple, "pvchi2dof", pvinfo[pvindex].totalChi2 / ( 2 * pvinfo[pvindex].nTracks - 3 ) ); + fillT( tuple, "pvxerr", std::sqrt( pv.covMatrix()( 0, 0 ) ) ); + fillT( tuple, "pvyerr", std::sqrt( pv.covMatrix()( 1, 1 ) ) ); + fillT( tuple, "pvzerr", std::sqrt( pv.covMatrix()( 2, 2 ) ) ); + // compute the track state at the z-position of the vertex + auto stateAtPV = state; + stateAtPV.linearTransportTo( pvz ); + const auto dx = stateAtPV.x() - pvx; + const auto dy = stateAtPV.y() - pvy; + // fill ipX and ipY. once this compiles, we also add the errors + fillT( tuple, "x", stateAtPV.x() ); + fillT( tuple, "y", stateAtPV.y() ); + fillT( tuple, "z", stateAtPV.z() ); + fillT( tuple, "xerr", std::sqrt( stateAtPV.covariance()( 0, 0 ) ) ); + fillT( tuple, "yerr", std::sqrt( stateAtPV.covariance()( 1, 1 ) ) ); + fillT( tuple, "ipx", dx ); + fillT( tuple, "ipy", dy ); + + // compute the contribution from the PV to ipX and ipX errors + ROOT::Math::SMatrix H; + H( 0, 0 ) = H( 1, 1 ) = 1; + H( 2, 0 ) = -stateAtPV.tx(); + H( 2, 1 ) = -stateAtPV.ty(); + const Gaudi::SymMatrix2x2 Vpv = ROOT::Math::Similarity( ROOT::Math::Transpose( H ), pv.covMatrix() ); + fillT( tuple, "pvipxerr", std::sqrt( Vpv( 0, 0 ) ) ); + fillT( tuple, "pvipyerr", std::sqrt( Vpv( 1, 1 ) ) ); + fillT( tuple, "ipxerr", std::sqrt( stateAtPV.covariance()( 0, 0 ) + Vpv( 0, 0 ) ) ); + fillT( tuple, "ipyerr", std::sqrt( stateAtPV.covariance()( 1, 1 ) + Vpv( 1, 1 ) ) ); + + // compute the IP chi2 (including the PV error) + const Gaudi::SymMatrix2x2 V = Vpv + stateAtPV.covariance().Sub( 0, 0 ); + const Gaudi::Vector2 residual{ dx, dy }; + Gaudi::SymMatrix2x2 W = V; + W.InvertChol(); + fillT( tuple, "ipchi2", ROOT::Math::Similarity( residual, W ) ); + + // compute the distance to the next PV. PVs should be ordered in z. + // the sign of deltazpv will tell if the next or previous is closer. + std::optional deltazpv; + if ( pvindex > 0 ) deltazpv = pvs[pvindex - 1].position().z() - pvz; + if ( pvindex < pvs.size() - 1 ) { + const auto nextpvdz = pvs[pvindex + 1].position().z() - pvz; + if ( !deltazpv || std::abs( nextpvdz ) < std::abs( *deltazpv ) ) deltazpv = nextpvdz; + } + fillT( tuple, "nextclosestpvdz", deltazpv ? *deltazpv : 9999. ); + } + tuple->write().ignore(); + } + } +} // namespace + +class TrackMonitorNT final : public LHCb::Algorithm::Consumer> { +public: + TrackMonitorNT( const std::string& name, ISvcLocator* pSvcLocator ) + : Consumer( name, pSvcLocator, + { KeyValue{ "TrackContainer", LHCb::TrackLocation::Velo }, + KeyValue{ "PVContainer", LHCb::Event::PV::DefaultLocation } } ) {} + void operator()( const LHCb::Track::Range& tracks, const Vertices& pvs ) const override { + return fillTuple( nTuple( "tracks", "", CLID_ColumnWiseTuple ), tracks, pvs ); + } +}; +DECLARE_COMPONENT( TrackMonitorNT ) diff --git a/Tr/TrackMonitors/src/TrackParticleMonitor.cpp b/Tr/TrackMonitors/src/TrackParticleMonitor.cpp index 4cd79c9fa578c694d222363727858b4cc378935f..3b50a7a20219e4c9e70a5f000d5507c86bffe8ef 100644 --- a/Tr/TrackMonitors/src/TrackParticleMonitor.cpp +++ b/Tr/TrackMonitors/src/TrackParticleMonitor.cpp @@ -124,6 +124,8 @@ private: mutable Gaudi::Accumulators::Histogram<1> track_ipY; mutable Gaudi::Accumulators::Histogram<1> vertex_chi2; mutable Gaudi::Accumulators::Histogram<1> vertex_chi2prob; + mutable Gaudi::Accumulators::Histogram<1> vertex_doca; + mutable Gaudi::Accumulators::Histogram<1> vertex_docapull; mutable Gaudi::Accumulators::Histogram<1> vertex_x_pos; mutable Gaudi::Accumulators::Histogram<1> vertex_y_pos; mutable Gaudi::Accumulators::Histogram<1> vertex_z_pos; @@ -149,8 +151,10 @@ private: { 50, -M_PI, M_PI } } , track_ipX{ owner, "trackIPx", "Track IPx", { 100, -0.1, 0.1, "IPx [mm]" } } , track_ipY{ owner, "trackIPy", "Track IPy", { 100, -0.1, 0.1, "IPy [mm]" } } - , vertex_chi2{ owner, "chi2", "vertex chi2", { 100, 0, 5, "#chi^2" } } + , vertex_chi2{ owner, "chi2", "vertex chi2", { 100, 0, 10, "#chi^2" } } , vertex_chi2prob{ owner, "chi2prob", "vertex chi2prob", { 100, 0, 1, "#chi^2 probability" } } + , vertex_doca{ owner, "doca", "two-prong doca", { 100, -0.15, 0.15 } } + , vertex_docapull{ owner, "docapull", "two-prong doca pull", { 100, -5, 5 } } , vertex_x_pos{ owner, "vtxx", "x position of vertex", { 100, -x_range, x_range, "x [mm]" } } , vertex_y_pos{ owner, "vtxy", "y position of vertex", { 100, -y_range, y_range, "y [mm]" } } , vertex_z_pos{ owner, "vtxz", "z position of vertex", { 100, -z_range, z_range, "z [mm]" } } @@ -187,7 +191,10 @@ private: mutable Gaudi::Accumulators::Histogram<1> matt; mutable Gaudi::Accumulators::Histogram<1> openingAngle; mutable Gaudi::Accumulators::Histogram<1> ipchi2; + mutable Gaudi::Accumulators::Histogram<1> decaylength; mutable Gaudi::Accumulators::Histogram<1> dls; + mutable Gaudi::Accumulators::Histogram<1> ipx; + mutable Gaudi::Accumulators::Histogram<1> ipy; BiasHistos( const TrackParticleMonitor* owner, float const& min_mass, float const& max_mass, float const& pdiff_range, unsigned int const& mass_bins, unsigned int const& pdiff_bins, @@ -229,7 +236,10 @@ private: { bins_angles * 5, 0, M_PI, "phiMatt [rad]" } } , openingAngle{ owner, "openingangle", "Opening angle", { bins_angles * 5, 0, 0.3, "Opening angle [rad]" } } , ipchi2{ owner, "ipchi2", "IP chi2", { 50, 0, 20 } } - , dls{ owner, "dls", "decay length significance", { 50, -5, 5 } } {} + , decaylength{ owner, "decaylength", "decay length [mm]", { 50, -1.0, 1.0 } } + , dls{ owner, "dls", "decay length significance", { 50, -5, 5 } } + , ipx{ owner, "ipx", "mother IP x", { 50, -0.1, 0.1 } } + , ipy{ owner, "ipy", "mother IP y", { 50, -0.1, 0.1 } } {} }; // Profile histograms of mass vs bias variables @@ -462,6 +472,10 @@ void TrackParticleMonitor::operator()( LHCb::Particle::Range const& particles, D double tx = p4.x() / p4.z(); double ty = p4.y() / p4.z(); + // For visualization we would like to have a 'doca-pull' + const auto twoprongdoca = LHCb::StateVertexUtils::doca( *states.front(), *states.back() ); + const auto twoprongchi2 = LHCb::StateVertexUtils::vertexChi2( *states.front(), *states.back() ); + // DECAY PLANE ANGLES // Unit vector perp. to the decay plane @@ -479,6 +493,9 @@ void TrackParticleMonitor::operator()( LHCb::Particle::Range const& particles, D // Vertex histograms ++track_vertex_histos.vertex_chi2[chi2]; ++track_vertex_histos.vertex_chi2prob[chi2prob]; + ++track_vertex_histos.vertex_doca[twoprongdoca]; + ++track_vertex_histos.vertex_docapull[std::sqrt( twoprongchi2 ) * ( twoprongdoca > 0 ? 1 : -1 )]; + ++track_vertex_histos.vertex_x_pos[vertex.position().x()]; ++track_vertex_histos.vertex_y_pos[vertex.position().y()]; ++track_vertex_histos.vertex_z_pos[vertex.position().z()]; @@ -521,12 +538,16 @@ void TrackParticleMonitor::operator()( LHCb::Particle::Range const& particles, D ++bias.matt[phimatt]; ++bias.openingAngle[openangle]; if ( particle->pv() ) { + const auto dx = vertex.position() - particle->pv()->position(); + ++bias.ipx[dx.x() - dx.z() * tx]; + ++bias.ipy[dx.y() - dx.z() * ty]; // FIXME: what a mess this has become. I also don't manage to call it anymore with the output of TrackStateVertex auto [chi2, decaylength, decaylength_err] = LHCb::StateVertexUtils::computeChiSquare( referencePoint( *particle ), threeMomentum( *particle ), LHCb::Event::covMatrix( *particle ), endVertexPos( *( particle->pv() ) ), posCovMatrix( *( particle->pv() ) ) ); if ( decaylength_err > 0 ) { ++bias.ipchi2[chi2.cast()]; + ++bias.decaylength[decaylength.cast()]; ++bias.dls[decaylength.cast() / decaylength_err.cast()]; } } diff --git a/Tr/TrackMonitors/src/TrackVertexMonitor.cpp b/Tr/TrackMonitors/src/TrackVertexMonitor.cpp index 985be3d208ffac4d695950d8bcf534a430d9fc72..8e5a7c33a80104899da49bbad1069400444b15f8 100644 --- a/Tr/TrackMonitors/src/TrackVertexMonitor.cpp +++ b/Tr/TrackMonitors/src/TrackVertexMonitor.cpp @@ -15,7 +15,10 @@ #include "Event/Track.h" #include "Event/TwoProngVertex.h" #include "LHCbDet/InteractionRegion.h" +#include "LHCbMath/StateVertexUtils.h" #include "TrackInterfaces/ITrackVertexer.h" +#include "TrackKernel/PrimaryVertexUtils.h" +#include "TrackKernel/TrackPVFinderUtils.h" #include "TrackKernel/TrackStateVertex.h" #include "GaudiKernel/PhysicalConstants.h" @@ -33,6 +36,74 @@ namespace { for ( const auto& track : tracks ) { states.push_back( &track->firstState() ); } return states; } + + void fillResolutionHistogram3SigmaRMS( const Gaudi::Accumulators::StaticHistogram<2>& hist2d, + Gaudi::Accumulators::ProfileHistogram<1>& prof1d ) { + // reset the target profile + // prof1d.reset() ; + const auto& gaudiAxisX = hist2d.template axis<0>(); + const auto minValueX = gaudiAxisX.minValue(); + const auto maxValueX = gaudiAxisX.maxValue(); + const auto nBinsX = gaudiAxisX.numBins(); + const auto& gaudiAxisY = hist2d.template axis<1>(); + const auto minValueY = gaudiAxisY.minValue(); + const auto maxValueY = gaudiAxisY.maxValue(); + const auto nBinsY = gaudiAxisY.numBins(); + const auto binXSize = ( maxValueX - minValueX ) / nBinsX; + const auto binYSize = ( maxValueY - minValueY ) / nBinsY; + + for ( unsigned int nx = 1; nx <= nBinsX; ++nx ) { + + // the offset of '1' is a bit nagging, so we first extract just the contents + double norm{ 0 }; + std::vector bincontents( nBinsY, 0 ); + for ( unsigned int ny = 1; ny <= nBinsY; ++ny ) { + const auto offset = ny * ( nBinsX + 2 ); + // auto const& [tmp, sumw2] = hist2d.binValue( nx + offset ); + // auto const& [nent, sumw] = tmp; + auto sumw = hist2d.binValue( nx + offset ); + bincontents[ny - 1] += sumw; + norm += sumw; + } + + if ( norm > 0 ) { + // now just start counting + double sumw( 0 ), sumx2w( 0 ); + double xtrunc( 0 ); + unsigned int ny{ 0 }; + for ( ; ny < nBinsY; ++ny ) { + const double x = ny + 0.5; + const double c = bincontents[ny]; + const double up = ny + 1; + const double newsumw = sumw + c; + double newsumx2w = sumx2w + c * x * x; + if ( newsumw > 0.5 * norm ) { + const double newrms = std::sqrt( newsumx2w / newsumw ); + if ( 3 * newrms < up ) { + const double drms = newrms - sqrt( sumx2w / sumw ); + const double frac = ( 3 * drms ); + xtrunc = ( up - ( 1 - frac ) ) / 3; + break; + } + } + sumw = newsumw; + sumx2w = newsumx2w; + } + + // if we ran out of the loop, just return the rms + if ( ny == nBinsY ) xtrunc = std::sqrt( sumx2w / sumw ); + + // get to a number in the right units + const double sigma = xtrunc * binYSize; + const double sigmaerr = sigma / std::sqrt( sumw ); + const auto binXValue = minValueX + ( nx - 0.5 ) * binXSize; + // unfortunately, in default error mode it determines the error from the spread of the values. + // the easiest way to fix it is to fill it twice. the factor sqrt(2) we found by trial and error. + prof1d[binXValue] += sigma + std::numbers::sqrt2 * sigmaerr; + prof1d[binXValue] += sigma - std::numbers::sqrt2 * sigmaerr; + } + } + } } // namespace class TrackVertexMonitor @@ -50,6 +121,12 @@ public: StatusCode initialize() override; + StatusCode stop() override { + fillResolutionHistogram3SigmaRMS( m_trackIPXInvPt, m_trackIPXResolutionVsInvPt ); + fillResolutionHistogram3SigmaRMS( m_trackIPYInvPt, m_trackIPYResolutionVsInvPt ); + return StatusCode::SUCCESS; + } + private: Gaudi::Property m_ipmax{ this, "MaxIP", 0.5 * Gaudi::Units::mm }; Gaudi::Property m_ipmaxprof{ this, "MaxIPProfile", 0.1 * Gaudi::Units::mm }; @@ -61,10 +138,21 @@ private: Gaudi::Property m_zpvmin_wide{ this, "MinZPV_Wide", -150 * Gaudi::Units::cm, "Wide z window for PV plot" }; Gaudi::Property m_zpvmax_wide{ this, "MaxZPV_Wide", 150 * Gaudi::Units::cm, "Wide z window for PV plot" }; Gaudi::Property m_maxLongTrackChisqPerDof{ this, "MaxLongTrackChisqPerDof", 5 }; - Gaudi::Property m_minLongTrackMomentum{ this, "MinLongTrackMomentum", 5 }; + Gaudi::Property m_minLongTrackMomentum{ this, "MinLongTrackMomentum", 5 * Gaudi::Units::GeV }; + Gaudi::Property m_minFastTrackMomentum{ this, "MinFastTrackMomentum", 50 * Gaudi::Units::GeV }; + Gaudi::Property m_minFastTrackPt{ this, "MinFastTrackPt", 2 * Gaudi::Units::GeV }; Gaudi::Property m_nprbins{ this, "NumProfileBins", 20 }; Gaudi::Property m_ntracksPV{ this, "NumTracksPV", 2 }; - Gaudi::Property m_produceHistogram{ this, "produceHistogram", false }; // producing IP 1/pt histograms + + Gaudi::Property m_minpvtotalweight{ this, "MinPVTotalWeight", 10. }; + Gaudi::Property m_maxpviperr{ this, "MaxPVIPError", 0.02 }; + Gaudi::Property m_maxiperr{ this, "MaxIPError", 0.10 }; + Gaudi::Property m_maxipchi2{ this, "MaxIPChi2", 50 }; + Gaudi::Property m_minTwoProngMass{ this, "MinTwoProngMass", 2 * Gaudi::Units::GeV }; + Gaudi::Property m_maxTwoProngDoca{ this, "MaxTwoProngDoca", 0.3 * Gaudi::Units::mm }; + Gaudi::Property m_maxTwoProngChi2{ this, "MaxTwoProngChi2", 16 }; + Gaudi::Property m_maxTwoProngIPChi2{ this, "MixTwoProngIPChi2", 10 }; + Gaudi::Property m_produceHistogram{ this, "produceHistogram", false }; // producing IP 1/pt histograms ToolHandle m_vertexer{ this, "TrackVertexer", "TrackVertexer" }; @@ -76,12 +164,22 @@ private: this, "NumBackTracksPerPV", "NumBackTracksPerPV", { 50, -0.5, 99.5 } }; mutable Gaudi::Accumulators::Histogram<1> m_pvChisquarePerDof{ this, "PV chisquare per dof", "PV chisquare per dof", { 150, 0., 3. } }; + mutable Gaudi::Accumulators::StaticHistogram<1> m_pvAverageWeight{ + this, "PV average weight", "PV average weight", { 50, 0., 1.0 } }; + mutable Gaudi::Accumulators::StaticProfileHistogram<1> m_pvAverageWeightVsNumTracksPr{ + this, "PV average weight vs num tracks", "PV average weight vs num tracks", { 75, 0.5, 150.5 } }; mutable Gaudi::Accumulators::Histogram<1> m_pvXPosition{ this, "PV x position", "PV x position", { 200, -m_xpvmax, m_xpvmax } }; mutable Gaudi::Accumulators::Histogram<1> m_pvYPosition{ this, "PV y position", "PV y position", { 200, -m_ypvmax, m_ypvmax } }; mutable Gaudi::Accumulators::Histogram<1> m_pvZPosition{ this, "PV z position", "PV z position", { 200, m_zpvmin, m_zpvmax } }; + mutable Gaudi::Accumulators::Histogram<1> m_pvXError{ + this, "PV x error", "PV x error", { 50, 0, 0.05 * Gaudi::Units::mm } }; + mutable Gaudi::Accumulators::Histogram<1> m_pvYError{ + this, "PV y error", "PV y error", { 50, 0, 0.05 * Gaudi::Units::mm } }; + mutable Gaudi::Accumulators::Histogram<1> m_pvZError{ + this, "PV z error", "PV z error", { 50, 0, 0.5 * Gaudi::Units::mm } }; mutable Gaudi::Accumulators::Histogram<1> m_pvbeamlineDeltaX{ this, "PV-beamline delta x", "PV-beamline delta x", { 200, -1, 1 } }; mutable Gaudi::Accumulators::Histogram<1> m_pvbeamlineDeltaY{ @@ -90,7 +188,6 @@ private: this, "PV-beamline delta x versus z", "PV-beamline delta x versus z", { m_nprbins, m_zpvmin, m_zpvmax } }; mutable Gaudi::Accumulators::ProfileHistogram<1> m_pvbeamlineDeltaYvsZ{ this, "PV-beamline delta y versus z", "PV-beamline delta y versus z", { m_nprbins, m_zpvmin, m_zpvmax } }; - mutable Gaudi::Accumulators::Histogram<1> m_pvZPositionWide{ this, "PV z position (wide)", "PV z position (wide)", { 200, m_zpvmin_wide, m_zpvmax_wide } }; mutable Gaudi::Accumulators::Histogram<1> m_pvLongChisquarePerDof{ @@ -137,24 +234,33 @@ private: this, "PV forward chisquare per dof", "PV forward chisquare per dof", { 50, 0, 10 } }; mutable Gaudi::Accumulators::Histogram<1> m_pvBackwardChisquareDof{ this, "PV backward chisquare per dof", "PV backward chisquare per dof", { 50, 0, 10 } }; - mutable Gaudi::Accumulators::Histogram<1> m_trackIPX{ - this, "track IP X", "track IP X (biased)", { 50, -m_ipmax, m_ipmax } }; - mutable Gaudi::Accumulators::Histogram<1> m_trackIPY{ - this, "track IP Y", "track IP Y (biased)", { 50, -m_ipmax, m_ipmax } }; - mutable Gaudi::Accumulators::Histogram<1> m_trackTransverseIP{ - this, "fast track transverse IP", "fast track transverse IP", { 50, -m_ipmax, m_ipmax } }; - mutable Gaudi::Accumulators::Histogram<1> m_trackLongitudinalIP{ - this, "fast track longitudinal IP", "fast track longitudinal IP", { 50, -m_ipmax, m_ipmax } }; + mutable Gaudi::Accumulators::Histogram<1> m_trackIPChi2{ this, "track IP chi2", "track IP chi2", { 50, 0, 25 } }; + mutable Gaudi::Accumulators::Histogram<1> m_trackIPX{ this, "track IP X", "track IP X", { 50, -m_ipmax, m_ipmax } }; + mutable Gaudi::Accumulators::Histogram<1> m_trackIPY{ this, "track IP Y", "track IP Y", { 50, -m_ipmax, m_ipmax } }; + mutable Gaudi::Accumulators::Histogram<1> m_trackIPXPull{ this, "track IP X pull", "track IP X pull", { 50, -5, 5 } }; + mutable Gaudi::Accumulators::Histogram<1> m_trackIPYPull{ this, "track IP Y pull", "track IP Y pull", { 50, -5, 5 } }; + mutable Gaudi::Accumulators::Histogram<1> m_fastTrackTransverseIP{ + this, "fast track transverse IP", "fast track transverse IP", { 50, -0.2, 0.2 } }; + mutable Gaudi::Accumulators::Histogram<1> m_fastTrackLongitudinalIP{ + this, "fast track longitudinal IP", "fast track longitudinal IP", { 50, -0.2, 0.2 } }; mutable Gaudi::Accumulators::Histogram<1> m_fastTrackIPX{ - this, "fast track IP X", "fast track IP X", { 50, -m_ipmax, m_ipmax } }; + this, "fast track IP X", "fast track IP X", { 50, -0.2, 0.2 } }; mutable Gaudi::Accumulators::Histogram<1> m_fastTrackIPY{ - this, "fast track IP Y", "fast track IP Y", { 50, -m_ipmax, m_ipmax } }; + this, "fast track IP Y", "fast track IP Y", { 50, -0.2, 0.2 } }; + mutable Gaudi::Accumulators::Histogram<1> m_fastTrackIPXPull{ + this, "fast track IP X pull", "fast track IP X pull", { 50, -5, 5 } }; + mutable Gaudi::Accumulators::Histogram<1> m_fastTrackIPYPull{ + this, "fast track IP Y pull", "fast track IP Y pull", { 50, -5, 5 } }; mutable Gaudi::Accumulators::Histogram<1> m_twoProngMass{ this, "twoprong mass (GeV)", "twoprong mass (GeV)", { 50, 0, 10 } }; mutable Gaudi::Accumulators::Histogram<1> m_twoProngMomentum{ this, "twoprong momentum (GeV)", "twoprong momentum (GeV)", { 50, 0, 200 } }; + mutable Gaudi::Accumulators::Histogram<1> m_twoProngIPX{ + this, "twoprong IP X", "twoprong IP X (mm)", { 50, -0.2, 0.2 } }; + mutable Gaudi::Accumulators::Histogram<1> m_twoProngIPY{ + this, "twoprong IP Y", "twoprong IP Y (mm)", { 50, -0.2, 0.2 } }; mutable Gaudi::Accumulators::Histogram<1> m_twoProngDoca{ - this, "twoprong doca (mm)", "twoprong doca (mm)", { 50, -5, 5 } }; + this, "twoprong doca (mm)", "twoprong doca (mm)", { 50, -0.3, 0.3 } }; mutable Gaudi::Accumulators::Histogram<1> m_twoProngDocaPull{ this, "twoprong doca pull", "twoprong doca pull", { 50, -5, 5 } }; mutable Gaudi::Accumulators::Histogram<1> m_twoProngDecayLength{ @@ -162,10 +268,8 @@ private: mutable Gaudi::Accumulators::Histogram<1> m_twoProngDecayLengthSig{ this, "twoprong decaylength significance", "twoprong decaylength significance", { 50, -5, 5 } }; mutable Gaudi::Accumulators::Histogram<1> m_twoProngCTau{ this, "twoprong ctau", "twoprong ctau", { 50, -0.1, 0.1 } }; - mutable Gaudi::Accumulators::Histogram<1> m_twoProngProperLifetime{ - this, "twoprong proper lifetime (ps)", "twoprong proper lifetime (ps)", { 50, -0.2, 0.2 } }; mutable Gaudi::Accumulators::Histogram<1> m_twoProngIPChi2PerDof{ - this, "twoprong IP chi2 per dof", "twoprong IP chi2 per dof", { 50, 0, 10 } }; + this, "twoprong IP chi2 per dof", "twoprong IP chi2 per dof", { 50, 0, m_maxTwoProngIPChi2 / 2. } }; mutable Gaudi::Accumulators::Histogram<1> m_numPrimaryVertices{ this, "NumPrimaryVertices", "NumPrimaryVertices", { 16, -0.5, 15.5 } }; @@ -192,47 +296,36 @@ private: "PV forward-backward delta x versus z", { m_nprbins, m_zpvmin, m_zpvmax } }; mutable Gaudi::Accumulators::ProfileHistogram<1> m_trackIPXvsPhi{ - this, "track IP X vs phi", "track IP X vs phi (biased)", { m_nprbins, -Gaudi::Units::pi, Gaudi::Units::pi } }; + this, "track IP X vs phi", "track IP X vs phi", { m_nprbins, -Gaudi::Units::pi, Gaudi::Units::pi } }; mutable Gaudi::Accumulators::ProfileHistogram<1> m_trackIPXvsEta{ - this, "track IP X vs eta", "track IP X vs eta (biased)", { m_nprbins, 2.0, 5.0 } }; - - // pt plots - mutable Gaudi::Accumulators::ProfileHistogram<1> m_trackIPXvsPt{ - this, "track IP X vs pt profile", "track IP X vs pt (GeV) (biased)", { 30, 0.0, 15.0 } }; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_trackIPXvsInversePt{ - this, "track IP X vs inverse pt profile", "track IP X vs 1/pt (1/GeV) (biased)", { m_nprbins, 0.0, 3.0 } }; - mutable Gaudi::Accumulators::Histogram<2> m_trackIPXInvPt{ - this, "track IP X vs inverse pt", "(biased) track IP X in 1/pt", { 20, 0.0, 3.0 }, { 50, -m_ipmax, m_ipmax } }; + this, "track IP X vs eta", "track IP X vs eta", { m_nprbins, 2.0, 5.0 } }; + mutable Gaudi::Accumulators::StaticHistogram<2> m_trackIPXInvPt{ + this, "track IP X vs inverse pt", "track IP X in 1/pt", { 20, 0.0, 3.0 }, { 150, 0, 0.3 } }; + mutable Gaudi::Accumulators::ProfileHistogram<1> m_trackIPXResolutionVsInvPt{ + this, "track IP X resolution vs inverse pt", "track IP X resolution vs inverse pt", { 20, 0.0, 3.0 } }; mutable Gaudi::Accumulators::ProfileHistogram<1> m_trackIPYvsPhi{ - this, "track IP Y vs phi", "track IP Y vs phi (biased)", { m_nprbins, -Gaudi::Units::pi, Gaudi::Units::pi } }; + this, "track IP Y vs phi", "track IP Y vs phi", { m_nprbins, -Gaudi::Units::pi, Gaudi::Units::pi } }; mutable Gaudi::Accumulators::ProfileHistogram<1> m_trackIPYvsEta{ - this, "track IP Y vs eta", "track IP Y vs eta (biased)", { m_nprbins, 2.0, 5.0 } }; - - // pt plots - mutable Gaudi::Accumulators::ProfileHistogram<1> m_trackIPYvsPt{ - this, "track IP Y vs pt profile", "track IP Y vs pt (GeV) (biased)", { 30, 0.0, 15.0 } }; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_trackIPYvsInversePt{ - this, - "track IP Y vs inverse pt profile", - "track IP Y vs in 1/pt range (1/GeV) (biased)", - { m_nprbins, 0.0, 3.0 } }; - mutable Gaudi::Accumulators::Histogram<2> m_trackIPYInvPt{ - this, "track IP Y vs inverse pt", "(biased) track IP Y in 1/pt", { 20, 0.0, 3.0 }, { 50, -m_ipmax, m_ipmax } }; + this, "track IP Y vs eta", "track IP Y vs eta", { m_nprbins, 2.0, 5.0 } }; + mutable Gaudi::Accumulators::StaticHistogram<2> m_trackIPYInvPt{ + this, "track IP Y vs inverse pt", "track IP Y in 1/pt", { 20, 0.0, 3.0 }, { 150, 0, 0.3 } }; + mutable Gaudi::Accumulators::ProfileHistogram<1> m_trackIPYResolutionVsInvPt{ + this, "track IP Y resolution vs inverse pt", "track IP Y resolution vs inverse pt", { 20, 0.0, 3.0 } }; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_trackTransverseIPvsPhi{ + mutable Gaudi::Accumulators::ProfileHistogram<1> m_fastTrackTransverseIPvsPhi{ this, "fast track transverse IP vs phi", "fast track transverse IP vs phi", { m_nprbins, -Gaudi::Units::pi, Gaudi::Units::pi } }; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_trackTransverseIPvsEta{ + mutable Gaudi::Accumulators::ProfileHistogram<1> m_fastTrackTransverseIPvsEta{ this, "fast track transverse IP vs eta", "fast track transverse IP vs eta", { m_nprbins, 2.0, 5.0 } }; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_trackLongitudinalIPvsPhi{ + mutable Gaudi::Accumulators::ProfileHistogram<1> m_fastTrackLongitudinalIPvsPhi{ this, "fast track longitudinal IP vs phi", "fast track longitudinal IP vs phi", { m_nprbins, -Gaudi::Units::pi, Gaudi::Units::pi } }; - mutable Gaudi::Accumulators::ProfileHistogram<1> m_trackLongitudinalIPvsEta{ + mutable Gaudi::Accumulators::ProfileHistogram<1> m_fastTrackLongitudinalIPvsEta{ this, "fast track longitudinal IP vs eta", "fast track longitudinal IP vs eta", { m_nprbins, 2.0, 5.0 } }; mutable Gaudi::Accumulators::ProfileHistogram<1> m_fastTrackIPXvsPhi{ this, "fast track IP X vs phi", "fast track IP X vs phi", { m_nprbins, -Gaudi::Units::pi, Gaudi::Units::pi } }; @@ -330,9 +423,8 @@ namespace { } // namespace -void TrackVertexMonitor::operator()( LHCb::RecVertex::Range const& pvcontainer, LHCb::Track::Range const& alltracks, - DetectorElement const& lhcb, - LHCb::Conditions::InteractionRegion const& ir ) const { +void TrackVertexMonitor::operator()( LHCb::RecVertex::Range const& recpvcontainer, LHCb::Track::Range const& alltracks, + DetectorElement const&, LHCb::Conditions::InteractionRegion const& ir ) const { const auto isLong = TrackTypePredicate( LHCb::Track::Types::Long ); const auto isBackward = TrackBackwardPredicate(); const auto isForward = TrackForwardPredicate(); @@ -345,9 +437,12 @@ void TrackVertexMonitor::operator()( LHCb::RecVertex::Range const& pvcontainer, // for now I'll just create the track lists from the Best container // number of primary vertices - ++m_numPrimaryVertices[pvcontainer.size()]; + ++m_numPrimaryVertices[recpvcontainer.size()]; - for ( const LHCb::RecVertex* pv : pvcontainer ) { + // First fill some information for the PVs. + // Also make a selection of the PVs suitable for IP computation. + std::vector isgoodpv; + for ( const auto* pv : recpvcontainer ) { auto tracks = myconvert( pv->tracks() ); auto forwardtracks = myselect( tracks, isForward ); auto backwardtracks = myselect( tracks, isBackward ); @@ -381,11 +476,26 @@ void TrackVertexMonitor::operator()( LHCb::RecVertex::Range const& pvcontainer, ++m_pvbeamlineDeltaY[pvblDeltaY]; if ( std::abs( pvblDeltaX ) < m_xpvmax ) m_pvbeamlineDeltaXvsZ[pv->position().z()] += pvblDeltaX; if ( std::abs( pvblDeltaY ) < m_ypvmax ) m_pvbeamlineDeltaYvsZ[pv->position().z()] += pvblDeltaY; + ++m_pvXError[std::sqrt( pv->covMatrix()( 0, 0 ) )]; + ++m_pvYError[std::sqrt( pv->covMatrix()( 1, 1 ) )]; + ++m_pvZError[std::sqrt( pv->covMatrix()( 2, 2 ) )]; } + bool thisisgoodpv = + std::sqrt( pv->covMatrix()( 0, 0 ) ) < m_maxpviperr && std::sqrt( pv->covMatrix()( 1, 1 ) ) < m_maxpviperr; + // average weight and total weight + if ( pv->weights().size() > 0 ) { + const double sumw = std::accumulate( pv->weights().begin(), pv->weights().end(), 0.0 ); + const double avw = sumw / pv->weights().size(); + ++m_pvAverageWeight[avw]; + m_pvAverageWeightVsNumTracksPr[pv->weights().size()] += avw; + thisisgoodpv = thisisgoodpv && sumw > m_minpvtotalweight; + } + isgoodpv.emplace_back( thisisgoodpv ); + // refit the primary vertex with only the long tracks if ( longtracks.size() >= 2 ) { - auto longvertex = m_vertexer->fit( firstStates( longtracks ), *lhcb.geometry() ); + auto longvertex = m_vertexer->fit( firstStates( longtracks ) ); if ( longvertex ) ++m_pvLongChisquarePerDof[longvertex->chi2() / longvertex->nDoF()]; } @@ -394,7 +504,7 @@ void TrackVertexMonitor::operator()( LHCb::RecVertex::Range const& pvcontainer, auto righttracks = myselect( tracks, TrackVeloSidePredicate( -1 ) ); if ( lefttracks.size() >= m_ntracksPV && righttracks.size() >= m_ntracksPV ) { // fit two vertices - auto leftvertex = m_vertexer->fit( firstStates( lefttracks ), *lhcb.geometry() ); + auto leftvertex = m_vertexer->fit( firstStates( lefttracks ) ); if ( leftvertex ) { ++m_pvLeftX[leftvertex->position().x()]; @@ -409,7 +519,7 @@ void TrackVertexMonitor::operator()( LHCb::RecVertex::Range const& pvcontainer, } */ } - auto rightvertex = m_vertexer->fit( firstStates( righttracks ), *lhcb.geometry() ); + auto rightvertex = m_vertexer->fit( firstStates( righttracks ) ); if ( rightvertex ) { ++m_pvRightX[rightvertex->position().x()]; ++m_pvRightY[rightvertex->position().y()]; @@ -475,8 +585,8 @@ void TrackVertexMonitor::operator()( LHCb::RecVertex::Range const& pvcontainer, if ( forwardtracks.size() >= 2 && backwardtracks.size() >= 2 ) { // fit two vertices - auto forwardvertex = m_vertexer->fit( firstStates( forwardtracks ), *lhcb.geometry() ); - auto backwardvertex = m_vertexer->fit( firstStates( backwardtracks ), *lhcb.geometry() ); + auto forwardvertex = m_vertexer->fit( firstStates( forwardtracks ) ); + auto backwardvertex = m_vertexer->fit( firstStates( backwardtracks ) ); if ( forwardvertex && backwardvertex ) { Gaudi::XYZVector dx = forwardvertex->position() - backwardvertex->position(); @@ -521,137 +631,204 @@ void TrackVertexMonitor::operator()( LHCb::RecVertex::Range const& pvcontainer, } } } + } - // for events with a single vertex, do something with IP of - // highest momentum track, as function of phi and eta. - if ( pvcontainer.size() == 1 && tracks.size() >= 10 ) { - - // now get all good long tracks from the best container: - auto goodlongtracks = myselect( alltracks, [&]( const LHCb::Track* tr ) { - return isLong( tr ) && tr->chi2PerDoF() < m_maxLongTrackChisqPerDof && tr->p() > m_minLongTrackMomentum; - } ); - - for ( const LHCb::Track* tr : goodlongtracks ) { - const LHCb::State& firststate = tr->firstState(); - double dz = pv->position().z() - firststate.z(); - double dx = firststate.x() + dz * firststate.tx() - pv->position().x(); - double dy = firststate.y() + dz * firststate.ty() - pv->position().y(); - Gaudi::XYZVector p3 = firststate.momentum(); - double pt = ( firststate.pt() / Gaudi::Units::GeV ); - double invPt = 1.0 / pt; - ++m_trackIPX[dx]; - ++m_trackIPY[dy]; - // apply a cut for the profiles - if ( std::abs( dx ) < m_ipmaxprof && std::abs( dy ) < m_ipmaxprof ) { - double phi = p3.phi(); - double eta = p3.eta(); - m_trackIPXvsEta[eta] += dx; - m_trackIPXvsPhi[phi] += dx; - m_trackIPYvsEta[eta] += dy; - m_trackIPYvsPhi[phi] += dy; - } - // profiles with no IP cut - m_trackIPXvsPt[pt] += dx; - m_trackIPXvsInversePt[invPt] += dx; - m_trackIPYvsPt[pt] += dy; - m_trackIPYvsInversePt[invPt] += dy; - - // single plots for 1/pt - ++m_trackIPXInvPt[{ invPt, dx }]; - ++m_trackIPYInvPt[{ invPt, dy }]; - } - - if ( goodlongtracks.size() >= 2 ) { - - std::sort( goodlongtracks.begin(), goodlongtracks.end(), []( const LHCb::Track* lhs, const LHCb::Track* rhs ) { - return lhs->firstState().pt() < rhs->firstState().pt(); - } ); - - const LHCb::Track* firsttrack = goodlongtracks.back(); - goodlongtracks.pop_back(); - - // now pick a 2nd track that makes the highest possible invariant mass with this one - double highestmass2( 0 ); - const LHCb::Track* secondtrack = nullptr; - Gaudi::XYZVector firstp3 = firsttrack->firstState().momentum(); - for ( const auto& t : goodlongtracks ) { - Gaudi::XYZVector p3 = t->firstState().momentum(); - double mass2 = p3.r() * firstp3.r() - p3.Dot( firstp3 ); - if ( secondtrack == 0 || highestmass2 < mass2 ) { - highestmass2 = mass2; - secondtrack = t; - } - } + // Compute track IP with respect to their best vertex. + // Because we need 'unbiasing' functionality later, let's convert the input container + LHCb::Event::PV::PrimaryVertexContainer pvcontainer; + LHCb::TrackPVFinder::populateFromRecVertices( + pvcontainer, LHCb::make_span( std::to_address( recpvcontainer.begin() ), recpvcontainer.size() ) ); + LHCb::TrackPVFinder::fitAdaptive( pvcontainer, LHCb::Event::PV::AdaptiveFitConfig{} ); - // recompute the vertex without these tracks - auto newend = tracks.end(); - newend = std::remove( tracks.begin(), newend, firsttrack ); - newend = std::remove( tracks.begin(), newend, secondtrack ); - tracks.erase( newend, tracks.end() ); - auto restvertex = m_vertexer->fit( firstStates( tracks ), *lhcb.geometry() ); - if ( restvertex && firsttrack->nStates() != 0 ) { - const LHCb::State& firststate = firsttrack->firstState(); - double dz = restvertex->position().z() - firststate.z(); - double dx = firststate.x() + dz * firststate.tx() - restvertex->position().x(); - double dy = firststate.y() + dz * firststate.ty() - restvertex->position().y(); - double nt = std::sqrt( firststate.tx() * firststate.tx() + firststate.ty() * firststate.ty() ); - // transverse and longitudinal impact parameter - double iptrans = ( dx * firststate.ty() - dy * firststate.tx() ) / nt; - double iplong = ( dx * firststate.tx() + dy * firststate.ty() ) / nt; - Gaudi::XYZVector p3 = firststate.momentum(); - double phi = p3.phi(); - double eta = p3.eta(); - - ++m_trackTransverseIP[iptrans]; - ++m_trackLongitudinalIP[iplong]; - ++m_fastTrackIPX[dx]; - ++m_fastTrackIPY[dy]; - // apply a cut for the profiles - if ( std::abs( iptrans ) < m_ipmaxprof && std::abs( iplong ) < m_ipmaxprof ) { - m_trackTransverseIPvsEta[eta] += iptrans; - m_trackTransverseIPvsPhi[phi] += iptrans; - m_trackLongitudinalIPvsEta[eta] += iplong; - m_trackLongitudinalIPvsPhi[phi] += iplong; - } - if ( std::abs( dx ) < m_ipmaxprof && std::abs( dy ) < m_ipmaxprof ) { - m_fastTrackIPXvsEta[eta] += dx; - m_fastTrackIPXvsPhi[phi] += dx; - m_fastTrackIPYvsEta[eta] += dy; - m_fastTrackIPYvsPhi[phi] += dy; - } + // get all good long tracks + auto goodlongtracks = myselect( alltracks, [&]( const LHCb::Track* tr ) { + return isLong( tr ) && tr->chi2PerDoF() < m_maxLongTrackChisqPerDof && tr->p() > m_minLongTrackMomentum; + } ); - // The two-track cuts we only make for relatively heavy objects - double mass = std::sqrt( highestmass2 ); - ++m_twoProngMass[mass / Gaudi::Units::GeV]; - if ( mass > 1 * Gaudi::Units::GeV ) { - // compute doca of two tracks - Gaudi::XYZVector dx3 = firsttrack->firstState().position() - secondtrack->firstState().position(); - Gaudi::XYZVector n3 = firsttrack->firstState().slopes().Cross( secondtrack->firstState().slopes() ); - double doca = dx3.Dot( n3 ) / n3.R(); - ++m_twoProngDoca[doca]; - if ( std::abs( doca ) < 200 ) { - m_twoProngDocavsEta[firstp3.eta()] += doca; - m_twoProngDocavsPhi[firstp3.phi()] += doca; - } - // the easiest way to compute the pull is with a vertex fit - auto twoprong = m_vertexer->fit( firsttrack->firstState(), secondtrack->firstState(), *lhcb.geometry() ); - if ( twoprong ) { - double pc = twoprong->p3().R(); - ++m_twoProngMomentum[pc / Gaudi::Units::GeV]; - ++m_twoProngDocaPull[std::sqrt( twoprong->chi2() ) * ( doca > 0 ? 1 : -1 )]; - double chi2, decaylength, decaylengtherr; - m_vertexer->computeDecayLength( *twoprong, *restvertex, chi2, decaylength, decaylengtherr ); - ++m_twoProngDecayLength[decaylength]; - ++m_twoProngDecayLengthSig[decaylength / decaylengtherr]; - ++m_twoProngIPChi2PerDof[chi2 / 2]; - ++m_twoProngCTau[decaylength * mass / pc]; - ++m_twoProngIPChi2PerDof[decaylength * mass / ( pc * Gaudi::Units::c_light * Gaudi::Units::picosecond )]; + const bool m_unbias_pvs_for_all_tracks{ false }; + std::vector bestpvindices; + bestpvindices.reserve( goodlongtracks.size() ); + for ( const LHCb::Track* tr : goodlongtracks ) { + const LHCb::State& state = tr->firstState(); + // find the best PV + const auto pvindex = LHCb::bestPVIndex( pvcontainer, state.x(), state.y(), state.z(), state.tx(), state.ty() ); + bestpvindices.push_back( pvindex ); + if ( pvindex != LHCb::Event::PV::PVIndexInvalid && isgoodpv[pvindex] ) { + + auto fillhistos = [&]( const auto& pv, bool fasttrack ) { + auto stateAtPV = state; + stateAtPV.linearTransportTo( pv.position().z() ); + const auto tx = stateAtPV.tx(); + const auto ty = stateAtPV.ty(); + ROOT::Math::SMatrix H; + H( 0, 0 ) = H( 1, 1 ) = 1; + H( 2, 0 ) = -tx; + H( 2, 1 ) = -ty; + const Gaudi::SymMatrix2x2 Vpv = ROOT::Math::Similarity( ROOT::Math::Transpose( H ), pv.covMatrix() ); + const auto pvipxerr = std::sqrt( Vpv( 0, 0 ) ); + const auto pvipyerr = std::sqrt( Vpv( 1, 1 ) ); + const auto trkipxerr = std::sqrt( stateAtPV.covariance()( 0, 0 ) ); + const auto trkipyerr = std::sqrt( stateAtPV.covariance()( 1, 1 ) ); + + if ( pvipxerr < m_maxpviperr && pvipyerr < m_maxpviperr ) { + + // compute the track state at the z-position of the vertex + const auto dx = stateAtPV.x() - pv.position().x(); + const auto dy = stateAtPV.y() - pv.position().y(); + const auto p3 = stateAtPV.momentum(); + const auto pt = ( stateAtPV.pt() / Gaudi::Units::GeV ); + const auto invPt = 1.0 / pt; + + // compute the IP chi2 + const Gaudi::SymMatrix2x2 V = Vpv + stateAtPV.covariance().Sub( 0, 0 ); + const Gaudi::Vector2 residual{ dx, dy }; + Gaudi::SymMatrix2x2 W = V; + W.InvertChol(); + const double ipchi2 = ROOT::Math::Similarity( residual, W ); + ++m_trackIPChi2[ipchi2]; + + // Make a very loose ipchi2 cut + if ( ipchi2 < m_maxipchi2 ) { + // single plots for 1/pt + // make a cut on the error such that the PV contribution is reasonably small. + const auto absdx = std::abs( dx ); + const auto absdy = std::abs( dy ); + if ( pvipxerr < 0.5 * trkipxerr ) ++m_trackIPXInvPt[{ invPt, absdx }]; + if ( pvipyerr < 0.5 * trkipyerr ) ++m_trackIPYInvPt[{ invPt, absdy }]; + + // making the next plots with a cut on ipx and ipy error, including the track contribution + const auto ipxerr = std::sqrt( V( 0, 0 ) ); + const auto ipyerr = std::sqrt( V( 1, 1 ) ); + const auto absdxpull = absdx / ipxerr; + const auto absdypull = absdy / ipyerr; + + if ( ipxerr < m_maxiperr && ipyerr < m_maxiperr ) { + if ( absdypull < 3 ) { + ++m_trackIPX[dx]; + ++m_trackIPXPull[dx / ipxerr]; + } + if ( absdxpull < 3 ) { + ++m_trackIPY[dy]; + ++m_trackIPYPull[dy / ipyerr]; + } + + // apply a cut for the profiles + const auto phi = p3.phi(); + const auto eta = p3.eta(); + if ( absdx < m_ipmaxprof && absdy < m_ipmaxprof ) { + m_trackIPXvsEta[eta] += dx; + m_trackIPXvsPhi[phi] += dx; + m_trackIPYvsEta[eta] += dy; + m_trackIPYvsPhi[phi] += dy; + } + + if ( fasttrack ) { + if ( absdypull < 3 ) { + ++m_fastTrackIPX[dx]; + ++m_fastTrackIPXPull[dx / ipxerr]; + } + if ( absdxpull < 3 ) { + ++m_fastTrackIPY[dy]; + ++m_fastTrackIPYPull[dy / ipyerr]; + } + const auto nt = std::sqrt( tx * tx + ty * ty ); + const auto iptrans = ( dx * ty - dy * tx ) / nt; + const auto iplong = ( dx * tx + dy * ty ) / nt; + ++m_fastTrackTransverseIP[iptrans]; + ++m_fastTrackLongitudinalIP[iplong]; + if ( std::abs( iptrans ) < m_ipmaxprof && std::abs( iplong ) < m_ipmaxprof ) { + m_fastTrackTransverseIPvsEta[eta] += iptrans; + m_fastTrackTransverseIPvsPhi[phi] += iptrans; + m_fastTrackLongitudinalIPvsEta[eta] += iplong; + m_fastTrackLongitudinalIPvsPhi[phi] += iplong; + } + if ( absdx < m_ipmaxprof && absdy < m_ipmaxprof ) { + m_fastTrackIPXvsEta[eta] += dx; + m_fastTrackIPXvsPhi[phi] += dx; + m_fastTrackIPYvsEta[eta] += dy; + m_fastTrackIPYvsPhi[phi] += dy; + } + } } } } + }; + // It is a bit too expensive to unbias PVs for all tracks. Therefore, we will only do that for the 'fast' tracks + const bool fasttrack = tr->firstState().p() > m_minFastTrackMomentum && tr->firstState().pt() > m_minFastTrackPt; + if ( m_unbias_pvs_for_all_tracks || fasttrack ) { + // extract the velo segment ID needed for unbiasing + const auto veloid = LHCb::Event::PV::uniqueVeloSegmentID( tr->lhcbIDs() ); + // put this into an array, and get the unbiased PV + std::array vetoedvelotracks{ veloid }; + const auto pv = LHCb::Event::PV::unbiasedVertex( pvcontainer, pvindex, vetoedvelotracks ); + fillhistos( pv, fasttrack ); + } else { + const auto& pv = pvcontainer[pvindex]; + fillhistos( pv, fasttrack ); } } } + + // Finally, fill some information on good two-prongs. To speed things up, we use every track at most once. + std::vector trackisused( goodlongtracks.size(), false ); + for ( unsigned itrk = 0; itrk < goodlongtracks.size(); ++itrk ) + if ( !trackisused[itrk] && bestpvindices[itrk] != LHCb::Event::PV::PVIndexInvalid && + isgoodpv[bestpvindices[itrk]] ) { + const LHCb::Track* trackA = goodlongtracks[itrk]; + const auto& stateA = trackA->firstState(); + const auto p3A = stateA.momentum(); + // pick a second track from the same PV + for ( unsigned jtrk = 0; jtrk < itrk; ++jtrk ) + if ( !trackisused[jtrk] && bestpvindices[itrk] == bestpvindices[jtrk] ) { + const LHCb::Track* trackB = goodlongtracks[jtrk]; + const auto& stateB = trackB->firstState(); + const auto p3B = stateB.momentum(); + const auto mass2 = p3A.r() * p3B.r() - p3A.Dot( p3B ); + if ( mass2 > m_minTwoProngMass * m_minTwoProngMass ) { + const auto doca = LHCb::StateVertexUtils::doca( stateA, stateB ); + if ( doca < m_maxTwoProngDoca ) { + const double chi2 = LHCb::StateVertexUtils::vertexChi2( stateA, stateB ); + if ( chi2 < m_maxTwoProngChi2 ) { + // the easiest way to compute the pull is with a vertex fit + auto twoprong = m_vertexer->fit( stateA, stateB ); + if ( twoprong ) { + // unbias the PV + auto veloid1 = LHCb::Event::PV::uniqueVeloSegmentID( trackA->lhcbIDs() ); + auto veloid2 = LHCb::Event::PV::uniqueVeloSegmentID( trackB->lhcbIDs() ); + // put this into an array, and get the unbiased PV + std::array vetoedvelotracks{ veloid1, veloid2 }; + const auto restvertex = + LHCb::Event::PV::unbiasedVertex( pvcontainer, bestpvindices[itrk], vetoedvelotracks ); + const auto deltapos = twoprong->position() - restvertex.position(); + const auto mom = p3A + p3B; + ++m_twoProngIPX[deltapos.x() - deltapos.z() * mom.x() / mom.z()]; + ++m_twoProngIPY[deltapos.y() - deltapos.z() * mom.y() / mom.z()]; + double ipchi2, decaylength, decaylengtherr; + m_vertexer->computeDecayLength( *twoprong, restvertex, ipchi2, decaylength, decaylengtherr ); + if ( ipchi2 < m_maxTwoProngIPChi2 ) { + const auto mass = std::sqrt( mass2 ); + ++m_twoProngMass[mass / Gaudi::Units::GeV]; + ++m_twoProngDoca[doca]; + const auto p3 = p3A + p3B; + const auto pc = p3.R(); + if ( std::abs( doca ) < m_ipmaxprof ) { + m_twoProngDocavsEta[p3.eta()] += doca; + m_twoProngDocavsPhi[p3.phi()] += doca; + } + ++m_twoProngMomentum[pc / Gaudi::Units::GeV]; + ++m_twoProngIPChi2PerDof[ipchi2 / 2]; + ++m_twoProngDocaPull[std::sqrt( twoprong->chi2() ) * ( doca > 0 ? 1 : -1 )]; + ++m_twoProngDecayLength[decaylength]; + ++m_twoProngDecayLengthSig[decaylength / decaylengtherr]; + ++m_twoProngCTau[decaylength * mass / pc]; + trackisused[itrk] = trackisused[jtrk] = true; + } + } + } + } + } + } + } } StatusCode TrackVertexMonitor::initialize() { diff --git a/Tr/TrackTools/src/TrackMuonMatching.cpp b/Tr/TrackTools/src/TrackMuonMatching.cpp index 7d19440e168bc4ba82227b12a06688029784480d..f6256f45e341c08bc98b23d2a87eb017a1c35ca7 100644 --- a/Tr/TrackTools/src/TrackMuonMatching.cpp +++ b/Tr/TrackTools/src/TrackMuonMatching.cpp @@ -69,12 +69,12 @@ private: Gaudi::Property m_matchAtZ{ this, "MatchAtZ", 12500 * Gaudi::Units::mm }; Gaudi::Property m_matchAtFirstMuonHit{ this, "MatchAtFirstMuonHit", false }; Gaudi::Property m_matchChi2Cut{ this, "MatchChi2Cut", 100.0 }; - Gaudi::Property m_minTrackSegmentChi2{ this, "TrackSegmentChi2", 5.0 }; + Gaudi::Property m_maxTrackSegmentChi2{ this, "TrackSegmentChi2", 5.0 }; Gaudi::Property m_allCombinations{ this, "AllCombinations", true }; Gaudi::Property m_fitTracks{ this, "FitTracks", true }; using Type = LHCb::Event::Enum::Track::Type; - Gaudi::Property m_tracktype{ this, "trackType", Type::Long }; + Gaudi::Property m_tracktype{ this, "trackType", Type::Long, "OBSOLETE" }; ToolHandle m_trackFitter{ this, "Fitter", "TrackMasterFitter/Fitter" }; ToolHandle m_extrapolator{ this, "Extrapolator", "TrackLinearExtrapolator" }; @@ -115,9 +115,9 @@ auto TrackMuonMatching::createMatchedTrack( const LHCb::Track& input_track, cons // Add LastMeasurement from muon track // and set the momentum of state LastFTHit from long track if ( muonTrack.hasStateAt( LHCb::State::Location::LastMeasurement ) ) { - matchedTrack->addToStates( *( muonTrack.stateAt( State::Location::LastMeasurement ) ) ); - matchedTrack->stateAt( State::Location::LastMeasurement ) - ->setQOverP( matchedTrack->stateAt( State::Location::LastFTHit )->qOverP() ); + auto s = *( muonTrack.stateAt( State::Location::LastMeasurement ) ); + s.setQOverP( matchedTrack->stateAt( State::Location::LastFTHit )->qOverP() ); + matchedTrack->addToStates( s ); } /// Add muon ids to copied T track for ( LHCbID id : muonTrack.lhcbIDs() ) matchedTrack->addToLhcbIDs( id ); @@ -151,7 +151,7 @@ LHCb::Tracks TrackMuonMatching::operator()( const InputTracks& tracks, const Inp /// Now match this T track to muon tracks for ( InputTracks::const_iterator m = muonTracks.begin(), mEnd = muonTracks.end(); m != mEnd; ++m ) { if ( msgLevel( MSG::DEBUG ) ) debug() << " MuonTrack chi2 " << ( *m )->chi2PerDoF() << endmsg; - if ( ( *m )->chi2PerDoF() > m_minTrackSegmentChi2.value() ) continue; + if ( ( *m )->chi2PerDoF() > m_maxTrackSegmentChi2.value() ) continue; if ( flaglongT ) { @@ -170,11 +170,8 @@ LHCb::Tracks TrackMuonMatching::operator()( const InputTracks& tracks, const Inp /// Matched Muon-T track auto best_matchedTrack = std::make_unique(); for ( InputTracks::const_iterator t = tracks.begin(), tEnd = tracks.end(); t != tEnd; ++t ) { - if ( !( *t )->checkType( m_tracktype ) ) continue; if ( !( *t )->hasT() ) continue; - if ( ( *t )->chi2PerDoF() > m_minTrackSegmentChi2.value() ) continue; - if ( !( *t )->checkType( LHCb::Track::Types::Long ) && !( *t )->checkType( LHCb::Track::Types::Ttrack ) ) - continue; + if ( ( *t )->chi2PerDoF() > m_maxTrackSegmentChi2.value() ) continue; double chi2 = -9999.0; /// Get the Track state closest to this z. Make a copy such that it can be changed. diff --git a/Tr/TrackTools/src/TrackVertexer.cpp b/Tr/TrackTools/src/TrackVertexer.cpp index ff78fa285bf796ce1a48da66cf8e738a6ab0abc6..6f979912f272cafb0946c922c466d70b8feab8fc 100644 --- a/Tr/TrackTools/src/TrackVertexer.cpp +++ b/Tr/TrackTools/src/TrackVertexer.cpp @@ -44,26 +44,24 @@ public: #endif /// Create a vertex from two track states - std::unique_ptr fit( LHCb::State const& stateA, LHCb::State const& stateB, - IGeometryInfo const& geometry ) const override; + std::unique_ptr fit( LHCb::State const& stateA, LHCb::State const& stateB ) const override; /// Create a veretx from a set of states - std::unique_ptr fit( LHCb::span states, - IGeometryInfo const& geometry ) const override; + std::unique_ptr fit( LHCb::span states ) const override; /// Create a vertex from a set of tracks. std::unique_ptr fit( LHCb::span tracks, IGeometryInfo const& geometry ) const override; /// Compute decaylength and IP chi2 wrt to PV. returns true if successful - bool computeDecayLength( const LHCb::TwoProngVertex& vertex, const LHCb::RecVertex& pv, double& chi2, + bool computeDecayLength( const LHCb::TwoProngVertex& vertex, const LHCb::VertexBase& pv, double& chi2, double& decaylength, double& decaylengtherr ) const override; /// Return the ip chi2 for a track (uses stateprovider) - double ipchi2( LHCb::Track const& track, LHCb::RecVertex const& pv, IGeometryInfo const& geometry ) const override; + double ipchi2( LHCb::Track const& track, LHCb::VertexBase const& pv, IGeometryInfo const& geometry ) const override; /// Return the ip chi2 for a track state - double ipchi2( LHCb::State const& state, LHCb::RecVertex const& pv, IGeometryInfo const& geometry ) const override; + double ipchi2( LHCb::State const& state, LHCb::VertexBase const& pv ) const override; private: ToolHandle m_stateprovider{ this, "StateProvider", "TrackStateProvider" }; @@ -79,8 +77,7 @@ private: // Declaration of the Tool Factory DECLARE_COMPONENT( TrackVertexer ) -std::unique_ptr TrackVertexer::fit( const LHCb::State& stateA, const LHCb::State& stateB, - IGeometryInfo const& ) const { +std::unique_ptr TrackVertexer::fit( const LHCb::State& stateA, const LHCb::State& stateB ) const { std::unique_ptr rc; std::array states{ &stateA, &stateB }; LHCb::TrackStateVertex vertex( states, m_maxDChisq, m_maxNumIter ); @@ -103,8 +100,7 @@ std::unique_ptr TrackVertexer::fit( const LHCb::State& sta return rc; } -std::unique_ptr TrackVertexer::fit( LHCb::span tracks, - IGeometryInfo const& ) const { +std::unique_ptr TrackVertexer::fit( LHCb::span tracks ) const { std::unique_ptr recvertex; if ( tracks.size() >= 2 ) { @@ -161,7 +157,7 @@ std::unique_ptr TrackVertexer::fit( LHCb::spanstate( z ) ); // fit the states - auto recvertex = fit( states, geometry ); + auto recvertex = fit( states ); // add the tracks if ( recvertex ) @@ -177,7 +173,7 @@ namespace { inline Gaudi::Vector3 transform( const Gaudi::XYZVector& vec ) { return Gaudi::Vector3( vec.X(), vec.Y(), vec.Z() ); } } // namespace -bool TrackVertexer::computeDecayLength( const LHCb::TwoProngVertex& vertex, const LHCb::RecVertex& pv, double& chi2, +bool TrackVertexer::computeDecayLength( const LHCb::TwoProngVertex& vertex, const LHCb::VertexBase& pv, double& chi2, double& decaylength, double& decaylengtherr ) const { // This calculation is basically a 1-iteration beamspot fit. The // constraint is @@ -229,13 +225,13 @@ bool TrackVertexer::computeDecayLength( const LHCb::TwoProngVertex& vertex, cons return ( OK != 0 ); } -double TrackVertexer::ipchi2( LHCb::Track const& track, LHCb::RecVertex const& pv, +double TrackVertexer::ipchi2( LHCb::Track const& track, LHCb::VertexBase const& pv, IGeometryInfo const& geometry ) const { const LHCb::TrackTraj* traj = m_stateprovider->trajectory( track, geometry ); - return traj ? ipchi2( traj->state( pv.position().z() ), pv, geometry ) : ipchi2( track.firstState(), pv, geometry ); + return traj ? ipchi2( traj->state( pv.position().z() ), pv ) : ipchi2( track.firstState(), pv ); } -double TrackVertexer::ipchi2( const LHCb::State& state, const LHCb::RecVertex& pv, IGeometryInfo const& ) const { +double TrackVertexer::ipchi2( const LHCb::State& state, const LHCb::VertexBase& pv ) const { double tx = state.tx(); double ty = state.ty(); double dz = pv.position().z() - state.z(); diff --git a/Tr/TrackUtils/src/TrackV0Finder.cpp b/Tr/TrackUtils/src/TrackV0Finder.cpp index b4f7b5897e17c31a0e9abb19ef38289d8c3adbd1..bee1c41be41121d43bd462490c019dc00752857f 100644 --- a/Tr/TrackUtils/src/TrackV0Finder.cpp +++ b/Tr/TrackUtils/src/TrackV0Finder.cpp @@ -231,7 +231,7 @@ OutContainer LHCb::TrackV0Finder::operator()( LHCb::Event::PV::PrimaryVertexCont State posstate = postraj.state( z ); State negstate = negtraj.state( z ); - std::unique_ptr vertex{ m_vertexer->fit( posstate, negstate, geometry ) }; + std::unique_ptr vertex{ m_vertexer->fit( posstate, negstate ) }; double decaylength; if ( !vertex ) { @@ -277,7 +277,7 @@ OutContainer LHCb::TrackV0Finder::operator()( LHCb::Event::PV::PrimaryVertexCont ++m_extrapolationFailure; negstate = negtraj.state( z ); } - std::unique_ptr newvertex{ m_vertexer->fit( posstate, negstate, geometry ) }; + std::unique_ptr newvertex{ m_vertexer->fit( posstate, negstate ) }; if ( newvertex ) vertex = std::move( newvertex ); }