From a086ffd6b7b086909cdaf9d4e37e1740e586c1b4 Mon Sep 17 00:00:00 2001
From: Wouter Hulsbergen <wouter.hulsbergen@nikhef.nl>
Date: Sat, 17 Feb 2024 00:04:08 +0100
Subject: [PATCH] added function to vertex container of states
---
.../RecEvent/include/Event/PrimaryVertices.h | 12 +--
Event/RecEvent/src/PrimaryVertices.cpp | 6 +-
.../include/Event/TrackVertexUtils.h | 34 +++++++-
Event/TrackEvent/src/TrackVertexUtils.cpp | 85 +++++++++----------
4 files changed, 84 insertions(+), 53 deletions(-)
diff --git a/Event/RecEvent/include/Event/PrimaryVertices.h b/Event/RecEvent/include/Event/PrimaryVertices.h
index e30b9ddf8c6..c278d201fed 100644
--- a/Event/RecEvent/include/Event/PrimaryVertices.h
+++ b/Event/RecEvent/include/Event/PrimaryVertices.h
@@ -169,12 +169,12 @@ namespace LHCb::Event::PV {
private:
using SymMatrix3x3 = LHCb::LinAlg::MatSym<double, 3>;
using Vector3 = LHCb::LinAlg::Vec<double, 3>;
- uint16_t m_begin{0}; /// Index of the first track of this PV in the track list
- uint16_t m_end{0}; /// Last+1 index of the first track of this PV in the track list
- SymMatrix3x3 m_halfD2Chi2DX2; /// Cached weight matrix, for unbiasing.
- Vector3 m_halfDChi2DX; /// Cached first derivative, for unbiasing.
- double m_sumOfWeights{0}; /// Cached sum of weights
- VeloSegmentIDMap m_veloidmap; /// veloids
+ uint16_t m_begin{0}; /// Index of the first track of this PV in the track list
+ uint16_t m_end{0}; /// Last+1 index of the first track of this PV in the track list
+ SymMatrix3x3 m_halfD2Chi2DX2; /// Cached weight matrix, for unbiasing.
+ Vector3 m_halfDChi2DX; /// Cached first derivative, for unbiasing.
+ double m_sumOfWeights{0}; /// Cached sum of weights
+ VeloSegmentIDMap m_veloidmap; /// veloids
public:
using PrimaryVertex::PrimaryVertex;
/// Constructor with a selected key (only for testing)
diff --git a/Event/RecEvent/src/PrimaryVertices.cpp b/Event/RecEvent/src/PrimaryVertices.cpp
index 4d2ab958f0d..3a4712fa735 100644
--- a/Event/RecEvent/src/PrimaryVertices.cpp
+++ b/Event/RecEvent/src/PrimaryVertices.cpp
@@ -88,7 +88,7 @@ namespace LHCb::Event::PV {
pv.setHalfD2Chi2DX2( accumulator.halfD2Chi2DX2 );
pv.setHalfDChi2DX( LinAlg::initialize_with_zeros<LHCb::LinAlg::Vec<double, 3>>() );
pv.setCovMatrix( convertToSMatrix<double>( cov ) );
- pv.setSumOfWeights( accumulator.sumw ) ;
+ pv.setSumOfWeights( accumulator.sumw );
return std::abs( delta( 2 ) ) < maxDeltaZConverged && abs( deltachi2 ) < maxDeltaChi2Converged;
};
@@ -222,7 +222,7 @@ namespace LHCb::Event::PV {
Chi2Accumulator accumulator = accumulate<true, true>( pv, pvtracks, chi2max );
pv.setHalfDChi2DX( accumulator.halfDChi2DX );
pv.setHalfD2Chi2DX2( accumulator.halfD2Chi2DX2 );
- pv.setSumOfWeights( accumulator.sumw ) ;
+ pv.setSumOfWeights( accumulator.sumw );
}
/// Fill the track fields that we do not read from the Velo tracks or the persistency. This routine is called by the
@@ -330,7 +330,7 @@ namespace LHCb::Event::PV {
weight * pvtrack.halfD2Chi2DX2( PVTrackTag::PosCovMatrixElement::v22 ).cast();
accumulator.chi2 -= weight * pvtrack.ipchi2().cast();
accumulator.entries -= 1;
- accumulator.sumw -= weight ;
+ accumulator.sumw -= weight;
}
}
}
diff --git a/Event/TrackEvent/include/Event/TrackVertexUtils.h b/Event/TrackEvent/include/Event/TrackVertexUtils.h
index d6568f924b7..6153eb0d203 100644
--- a/Event/TrackEvent/include/Event/TrackVertexUtils.h
+++ b/Event/TrackEvent/include/Event/TrackVertexUtils.h
@@ -30,11 +30,20 @@ namespace LHCb::TrackVertexUtils {
/////////////////////////////////////////////////////////////////////////
/// Add a track to a vertex represented by a position and a weight
/// matrix. Both the weight matrix and the covariance matrix are
- /// computed. Returns the chi2
+ /// computed. Make sure to initialize vertexweight with zeros. Returns the chi2.
/////////////////////////////////////////////////////////////////////////
double addToVertex( const LHCb::State& state, Gaudi::XYZPoint& vertexpos, Gaudi::SymMatrix3x3& vertexweight,
Gaudi::SymMatrix3x3& vertexcov );
+ /////////////////////////////////////////////////////////////////////////
+ /// Add a range of tracks to a vertex represented by a position and a weight
+ /// matrix. Both the weight matrix and the covariance matrix are
+ /// computed. Make sure to initialize vertexweight with zeros. Returns the chi2.
+ /////////////////////////////////////////////////////////////////////////
+ template <typename StateContainer>
+ double addToVertex( const StateContainer& states, Gaudi::XYZPoint& vertexpos, Gaudi::SymMatrix3x3& vertexweight,
+ Gaudi::SymMatrix3x3& vertexcov );
+
/////////////////////////////////////////////////////////////////////////
/// Computes a vertex from two track states.
/////////////////////////////////////////////////////////////////////////
@@ -48,3 +57,26 @@ namespace LHCb::TrackVertexUtils {
Gaudi::SymMatrix3x3& vertexcov );
} // namespace LHCb::TrackVertexUtils
+
+// Below are implementations.
+namespace LHCb::TrackVertexUtils {
+
+ /// Helper function to compute chi2 derivatives for vertex fit. See .cpp file.
+ double addToDerivatives( const LHCb::State& state, const Gaudi::XYZPoint& vertexpos, Gaudi::Vector3& halfDChi2DX,
+ Gaudi::SymMatrix3x3& halfD2Chi2DX2 );
+
+ /// Helper function to compute vertex position from derivatives. See .cpp file.
+ double solve( const Gaudi::Vector3& halfDChi2DX, const Gaudi::SymMatrix3x3& halfD2Chi2DX2, Gaudi::XYZPoint& vertexpos,
+ Gaudi::SymMatrix3x3& vertexcov );
+
+ template <typename StateContainer>
+ double addToVertex( const StateContainer& states, Gaudi::XYZPoint& vertexpos, Gaudi::SymMatrix3x3& vertexweight,
+ Gaudi::SymMatrix3x3& vertexcov ) {
+ Gaudi::Vector3 halfDChi2DX{}; // does default initializer initialize to zero? not for libEIGEN!
+ double chi2{0};
+ for ( const auto& state : states ) chi2 += addToDerivatives( state, vertexpos, halfDChi2DX, vertexweight );
+ // compute the vertex
+ chi2 += solve( halfDChi2DX, vertexweight, vertexpos, vertexcov );
+ return chi2;
+ }
+} // namespace LHCb::TrackVertexUtils
diff --git a/Event/TrackEvent/src/TrackVertexUtils.cpp b/Event/TrackEvent/src/TrackVertexUtils.cpp
index 1585502f293..ab391bf508a 100644
--- a/Event/TrackEvent/src/TrackVertexUtils.cpp
+++ b/Event/TrackEvent/src/TrackVertexUtils.cpp
@@ -18,51 +18,50 @@ namespace LHCb::TrackVertexUtils {
/// to a vertex. This should probably go into LHCb math.
/////////////////////////////////////////////////////////////////////////
Gaudi::Vector3 transform( const Gaudi::XYZVector& vec ) { return Gaudi::Vector3( vec.X(), vec.Y(), vec.Z() ); }
- namespace {
- double addToDerivatives( const LHCb::State& state, const Gaudi::XYZPoint& vertexpos, Gaudi::Vector3& halfDChi2DX,
- Gaudi::SymMatrix3x3& halfD2Chi2DX2 ) {
- // compute residual
- Gaudi::Vector2 res;
- double dz = vertexpos.z() - state.z();
- res( 0 ) = state.x() + dz * state.tx() - vertexpos.x();
- res( 1 ) = state.y() + dz * state.ty() - vertexpos.y();
- // compute the weight matrix
- const auto& trkcov = state.covariance();
- Gaudi::SymMatrix2x2 invcov = trkcov.Sub<Gaudi::SymMatrix2x2>( 0, 0 );
- // extrapolate to the vertex
- invcov( 0, 0 ) += dz * dz * trkcov( 2, 2 ) + 2 * dz * trkcov( 2, 0 );
- invcov( 1, 0 ) += dz * dz * trkcov( 3, 2 ) + dz * ( trkcov( 3, 0 ) + trkcov( 2, 1 ) );
- invcov( 1, 1 ) += dz * dz * trkcov( 3, 3 ) + 2 * dz * trkcov( 3, 1 );
- invcov.Invert();
- // I tried to make this faster by writing it out, but H does
- // not contain sufficiently manby zeroes. Better to
- // parallelize.
- ROOT::Math::SMatrix<double, 3, 2> H;
- H( 0, 0 ) = H( 1, 1 ) = 1;
- H( 2, 0 ) = -state.tx();
- H( 2, 1 ) = -state.ty();
- halfD2Chi2DX2 += Similarity( H, invcov );
- halfDChi2DX += ( H * invcov ) * res;
- // You could potentially save time by reusing HW. However, it
- // does not help enough. Perhaps better when parallelized.
- // ROOT::Math::SMatrix<double,3,2> HW = H*invcov ;
- // Gaudi::SymMatrix3x3 HWH ;
- // ROOT::Math::AssignSym::Evaluate( HWH, HW*Transpose(H) ) ;
- // halfD2Chi2DX2 += HWH ;
- // halfDChi2DX += HW * res ;
- return Similarity( res, invcov );
- }
+ double addToDerivatives( const LHCb::State& state, const Gaudi::XYZPoint& vertexpos, Gaudi::Vector3& halfDChi2DX,
+ Gaudi::SymMatrix3x3& halfD2Chi2DX2 ) {
+ // compute residual
+ Gaudi::Vector2 res;
+ double dz = vertexpos.z() - state.z();
+ res( 0 ) = state.x() + dz * state.tx() - vertexpos.x();
+ res( 1 ) = state.y() + dz * state.ty() - vertexpos.y();
+ // compute the weight matrix
+ const auto& trkcov = state.covariance();
+ Gaudi::SymMatrix2x2 invcov = trkcov.Sub<Gaudi::SymMatrix2x2>( 0, 0 );
+ // extrapolate to the vertex
+ invcov( 0, 0 ) += dz * dz * trkcov( 2, 2 ) + 2 * dz * trkcov( 2, 0 );
+ invcov( 1, 0 ) += dz * dz * trkcov( 3, 2 ) + dz * ( trkcov( 3, 0 ) + trkcov( 2, 1 ) );
+ invcov( 1, 1 ) += dz * dz * trkcov( 3, 3 ) + 2 * dz * trkcov( 3, 1 );
+ invcov.Invert();
- double solve( const Gaudi::Vector3& halfDChi2DX, const Gaudi::SymMatrix3x3& halfD2Chi2DX2,
- Gaudi::XYZPoint& vertexpos, Gaudi::SymMatrix3x3& vertexcov ) {
- vertexcov = halfD2Chi2DX2;
- vertexcov.InvertChol();
- Gaudi::Vector3 delta = vertexcov * halfDChi2DX;
- vertexpos = {vertexpos.x() + delta( 0 ), vertexpos.y() + delta( 1 ), vertexpos.z() + delta( 2 )};
- return -1 * Dot( delta, halfDChi2DX );
- }
- } // namespace
+ // I tried to make this faster by writing it out, but H does
+ // not contain sufficiently manby zeroes. Better to
+ // parallelize.
+ ROOT::Math::SMatrix<double, 3, 2> H;
+ H( 0, 0 ) = H( 1, 1 ) = 1;
+ H( 2, 0 ) = -state.tx();
+ H( 2, 1 ) = -state.ty();
+ halfD2Chi2DX2 += Similarity( H, invcov );
+ halfDChi2DX += ( H * invcov ) * res;
+ // You could potentially save time by reusing HW. However, it
+ // does not help enough. Perhaps better when parallelized.
+ // ROOT::Math::SMatrix<double,3,2> HW = H*invcov ;
+ // Gaudi::SymMatrix3x3 HWH ;
+ // ROOT::Math::AssignSym::Evaluate( HWH, HW*Transpose(H) ) ;
+ // halfD2Chi2DX2 += HWH ;
+ // halfDChi2DX += HW * res ;
+ return Similarity( res, invcov );
+ }
+
+ double solve( const Gaudi::Vector3& halfDChi2DX, const Gaudi::SymMatrix3x3& halfD2Chi2DX2, Gaudi::XYZPoint& vertexpos,
+ Gaudi::SymMatrix3x3& vertexcov ) {
+ vertexcov = halfD2Chi2DX2;
+ vertexcov.InvertChol();
+ Gaudi::Vector3 delta = vertexcov * halfDChi2DX;
+ vertexpos = {vertexpos.x() + delta( 0 ), vertexpos.y() + delta( 1 ), vertexpos.z() + delta( 2 )};
+ return -1 * Dot( delta, halfDChi2DX );
+ }
double vertex( const LHCb::State& stateA, const LHCb::State& stateB, Gaudi::XYZPoint& vertexpos,
Gaudi::SymMatrix3x3& vertexweight, Gaudi::SymMatrix3x3& vertexcov ) {
--
GitLab