diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp index bb48156a703..f05f0929bea 100644 --- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp +++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp @@ -111,10 +111,11 @@ class AdaptiveMultiVertexFinder { bool doRealMultiVertex = true; // Decides if you want to use the ```vertexCompatibility``` of the - // track (set to true) or the ```chi2Track``` (set to false) as an + // track (set to true) or the ```chi2``` of the TrackAtVertex (set to + // false) as an // estimate for a track being an outlier or not. // In case the track refitting is switched on in the AMVFitter, you - // may want to use the refitted ```chi2Track```. + // may want to use the refitted ```chi2```. bool useFastCompatibility = true; // Maximum significance on the distance between two vertices diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.ipp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.ipp index b506afd55d7..e02475d5386 100644 --- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.ipp +++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.ipp @@ -247,7 +247,7 @@ auto Acts::AdaptiveMultiVertexFinder:: if (ipSig < m_cfg.tracksMaxSignificance) { // Create TrackAtVertex objects, unique for each (track, vertex) pair fitterState.tracksAtVerticesMap.emplace(std::make_pair(trk, &vtx), - TrackAtVertex(params, trk)); + TrackAtVertex(trk)); // Add the original track parameters to the list for vtx fitterState.vtxInfoMap[&vtx].trackLinks.push_back(trk); @@ -361,9 +361,8 @@ auto Acts::AdaptiveMultiVertexFinder:: fitterState.tracksAtVerticesMap.at(std::make_pair(trk, &vtx)); if ((trkAtVtx.vertexCompatibility < m_cfg.maxVertexChi2 && m_cfg.useFastCompatibility) || - (trkAtVtx.trackWeight > m_cfg.minWeight && - trkAtVtx.chi2Track < m_cfg.maxVertexChi2 && - !m_cfg.useFastCompatibility)) { + (trkAtVtx.weight > m_cfg.minWeight && + trkAtVtx.chi2 < m_cfg.maxVertexChi2 && !m_cfg.useFastCompatibility)) { // TODO: Understand why looking for compatible tracks only in seed tracks // and not also in all tracks auto foundIter = @@ -399,9 +398,8 @@ auto Acts::AdaptiveMultiVertexFinder:: fitterState.tracksAtVerticesMap.at(std::make_pair(trk, &vtx)); if ((trkAtVtx.vertexCompatibility < m_cfg.maxVertexChi2 && m_cfg.useFastCompatibility) || - (trkAtVtx.trackWeight > m_cfg.minWeight && - trkAtVtx.chi2Track < m_cfg.maxVertexChi2 && - !m_cfg.useFastCompatibility)) { + (trkAtVtx.weight > m_cfg.minWeight && + trkAtVtx.chi2 < m_cfg.maxVertexChi2 && !m_cfg.useFastCompatibility)) { // Find and remove track from seedTracks auto foundSeedIter = std::find_if(seedTracks.begin(), seedTracks.end(), @@ -483,7 +481,7 @@ auto Acts::AdaptiveMultiVertexFinder::keepNewVertex( for (const auto& trk : fitterState.vtxInfoMap[&vtx].trackLinks) { const auto& trkAtVtx = fitterState.tracksAtVerticesMap.at(std::make_pair(trk, &vtx)); - double trackWeight = trkAtVtx.trackWeight; + double trackWeight = trkAtVtx.weight; contaminationNum += trackWeight * (1. - trackWeight); contaminationDeNom += trackWeight * trackWeight; } diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp index 2c2a3969180..71f141491ea 100644 --- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp +++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp @@ -287,12 +287,12 @@ Acts::Result Acts:: for (const auto& trk : vtxInfo.trackLinks) { auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx)); - // Set trackWeight for current track - trkAtVtx.trackWeight = m_cfg.annealingTool.getWeight( + // Set weight for current track + trkAtVtx.weight = m_cfg.annealingTool.getWeight( state.annealingState, trkAtVtx.vertexCompatibility, collectTrackToVertexCompatibilities(state, trk)); - if (trkAtVtx.trackWeight > m_cfg.minWeight) { + if (trkAtVtx.weight > m_cfg.minWeight) { // Check if track is already linearized and whether we need to // relinearize if (!trkAtVtx.isLinearized || vtxInfo.relinearize) { @@ -368,7 +368,7 @@ void Acts::AdaptiveMultiVertexFitter< for (const auto vtx : state.vertexCollection) { for (const auto trk : state.vtxInfoMap[vtx].trackLinks) { auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx)); - if (trkAtVtx.trackWeight > m_cfg.minWeight) { + if (trkAtVtx.weight > m_cfg.minWeight) { KalmanVertexTrackUpdater::update(trkAtVtx, *vtx); } } diff --git a/Core/include/Acts/Vertexing/FullBilloirVertexFitter.ipp b/Core/include/Acts/Vertexing/FullBilloirVertexFitter.ipp index be254c7910c..c8f20b9aa6b 100644 --- a/Core/include/Acts/Vertexing/FullBilloirVertexFitter.ipp +++ b/Core/include/Acts/Vertexing/FullBilloirVertexFitter.ipp @@ -214,11 +214,13 @@ Acts::FullBilloirVertexFitter::fit( Vector4 deltaV = covV * deltaVFac; //-------------------------------------------------------------------------------------- // start momentum related calculations - - std::vector> covDeltaP(nTracks); + // Cross-covariance matrix of the vertex and the momentum + std::vector> crossCovVP(nTracks); + // Covariance matrix of the momentum + std::vector> covP(nTracks); // Update track momenta and calculate the covariance of the track parameters - // after the fit (TODO: parameters -> momenta). + // after the fit. for (std::size_t iTrack = 0; iTrack < billoirTracks.size(); ++iTrack) { auto& billoirTrack = billoirTracks[iTrack]; @@ -241,43 +243,15 @@ Acts::FullBilloirVertexFitter::fit( // ... and add it to the total chi2 value newChi2 += billoirTrack.chi2; - // calculate the 6x6 cross-covariance matrix - // TODO: replace fittedParams with fittedMomentum since d0 and z0 are - // ill-defined. At the moment, only the submatrix of momentum covariances - // is correct. - - // coordinate transformation matrix, i.e., - // d(d0,z0,phi,theta,qOverP,t)/d(x,y,z,t,phi,theta,qOverP) - // TODO: This is not correct since the (x, y, z, t) in the derivatives in - // the D matrix correspond to the global track position at the PCA rather - // than the 4D vertex position. - ActsMatrix transMat; - transMat.setZero(); - transMat.block<2, 2>(0, 0) = billoirTrack.D.template block<2, 2>(0, 0); - transMat(1, 2) = 1.; - transMat(2, 4) = 1.; - transMat(3, 5) = 1.; - transMat(4, 6) = 1.; - transMat(5, 3) = 1.; - - // cov(V,P) - // TODO: This is incorrect (see Ref. (2)), but it will not be needed - // anyway once we replace fittedParams with fittedMomentum - ActsMatrix<4, 3> covVP = billoirTrack.B; - - // cov(P,P), 3x3 matrix - ActsSquareMatrix<3> covP = - billoirTrack.Cinv + - billoirTrack.BCinv.transpose() * covV * billoirTrack.BCinv; - - ActsSquareMatrix<7> cov; - cov.setZero(); - cov.block<4, 4>(0, 0) = covV; - cov.block<4, 3>(0, 4) = covVP; - cov.block<3, 4>(4, 0) = covVP.transpose(); - cov.block<3, 3>(4, 4) = covP; - - covDeltaP[iTrack] = transMat * cov * transMat.transpose(); + // Cross-covariance matrix between the vertex position and the refitted + // momentum, see Eq. 8.22 in Ref. (2). Be mindful of the different choice + // of notation! + crossCovVP[iTrack] = -covV * billoirTrack.BCinv; + + // Covariance matrix of the refitted momentum, see Eq. 8.23 in Ref. (2). + // Be mindful of the different choice of notation! + covP[iTrack] = billoirTrack.Cinv - + billoirTrack.BCinv.transpose() * crossCovVP[iTrack]; } // assign new linearization point (= new vertex position in global frame) @@ -317,23 +291,12 @@ Acts::FullBilloirVertexFitter::fit( for (std::size_t iTrack = 0; iTrack < billoirTracks.size(); ++iTrack) { const auto& billoirTrack = billoirTracks[iTrack]; - - // TODO: replace refittedParams with refittedMomenta since d0 and z0 - // after the vertex fit are ill-defined (FullBilloir only outputs the - // updated track momenta) - - // new refitted trackparameters - BoundVector paramVec = BoundVector::Zero(); - paramVec[eBoundPhi] = trackMomenta[iTrack](0); - paramVec[eBoundTheta] = trackMomenta[iTrack](1); - paramVec[eBoundQOverP] = trackMomenta[iTrack](2); - paramVec[eBoundTime] = linPoint[FreeIndices::eFreeTime]; - BoundTrackParameters refittedParams( - perigee, paramVec, covDeltaP[iTrack], - extractParameters(*billoirTrack.originalTrack) - .particleHypothesis()); - TrackAtVertex trackAtVertex( - billoirTrack.chi2, refittedParams, billoirTrack.originalTrack); + // new refitted track momentum + FittedMomentum fittedMomentum(trackMomenta[iTrack], covP[iTrack], + crossCovVP[iTrack]); + TrackAtVertex trackAtVertex(billoirTrack.originalTrack, + std::move(fittedMomentum), + billoirTrack.chi2); tracksAtVertex.push_back(std::move(trackAtVertex)); } diff --git a/Core/include/Acts/Vertexing/IterativeVertexFinder.ipp b/Core/include/Acts/Vertexing/IterativeVertexFinder.ipp index b1caa4df227..581b6b43838 100644 --- a/Core/include/Acts/Vertexing/IterativeVertexFinder.ipp +++ b/Core/include/Acts/Vertexing/IterativeVertexFinder.ipp @@ -260,7 +260,7 @@ Acts::IterativeVertexFinder::removeUsedCompatibleTracks( for (const auto& trackAtVtx : tracksAtVertex) { // Check compatibility - if (trackAtVtx.trackWeight < m_cfg.cutOffTrackWeight) { + if (trackAtVtx.weight < m_cfg.cutOffTrackWeight) { // Do not remove track here, since it is not compatible with the vertex continue; } @@ -449,7 +449,7 @@ Acts::IterativeVertexFinder::reassignTracksToNewVertex( for (auto tracksIter = tracksBegin; tracksIter != tracksEnd;) { // consider only tracks that are not too tightly assigned to other // vertex - if (tracksIter->trackWeight > m_cfg.cutOffTrackWeightReassign) { + if (tracksIter->weight > m_cfg.cutOffTrackWeightReassign) { tracksIter++; continue; } @@ -560,6 +560,6 @@ int Acts::IterativeVertexFinder::countSignificantTracks( const Vertex& vtx) const { return std::count_if(vtx.tracks().begin(), vtx.tracks().end(), [this](TrackAtVertex trk) { - return trk.trackWeight > m_cfg.cutOffTrackWeight; + return trk.weight > m_cfg.cutOffTrackWeight; }); } diff --git a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp index a501f22a84e..7c94899f1d1 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp +++ b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp @@ -28,24 +28,6 @@ template void update(TrackAtVertex& track, const Vertex& vtx); -namespace detail { - -/// @brief reates a new covariance matrix for the -/// refitted track parameters -/// -/// @param sMat Track ovariance in momentum space -/// @param newTrkCov New track covariance matrixs -/// @param vtxWeight Vertex weight matrix -/// @param vtxCov Vertex covariance matrix -/// @param newTrkParams New track parameter -inline BoundMatrix createFullTrackCovariance(const SquareMatrix3& sMat, - const ActsMatrix<4, 3>& newTrkCov, - const SquareMatrix4& vtxWeight, - const SquareMatrix4& vtxCov, - const BoundVector& newTrkParams); - -} // Namespace detail - } // Namespace KalmanVertexTrackUpdater } // Namespace Acts diff --git a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp index 1f9873986ae..21f51e9224a 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp @@ -43,29 +43,21 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, Vector3 newTrkMomentum = sMat * momJac.transpose() * trkParamWeight * (trkParams - residual - posJac * vtxPos); - // Refit track parameters - BoundVector newTrkParams(BoundVector::Zero()); - // Get phi and theta and correct for possible periodicity changes const auto correctedPhiTheta = Acts::detail::normalizePhiTheta(newTrkMomentum(0), newTrkMomentum(1)); - newTrkParams(BoundIndices::eBoundPhi) = correctedPhiTheta.first; // phi - newTrkParams(BoundIndices::eBoundTheta) = correctedPhiTheta.second; // theta - newTrkParams(BoundIndices::eBoundQOverP) = newTrkMomentum(2); // qOverP + newTrkMomentum(0) = correctedPhiTheta.first; + newTrkMomentum(1) = correctedPhiTheta.second; // Vertex covariance and weight matrices const SquareMatrix3 vtxCov = vtx.fullCovariance().template block<3, 3>(0, 0); const SquareMatrix3 vtxWeight = vtxCov.inverse(); - // New track covariance matrix - const SquareMatrix3 newTrkCov = - -vtxCov * posJac.transpose() * trkParamWeight * momJac * sMat; - KalmanVertexUpdater::MatrixCache matrixCache; // Now determine the smoothed chi2 of the track in the following KalmanVertexUpdater::updatePosition( - vtx, linTrack, track.trackWeight, -1, matrixCache); + vtx, linTrack, track.weight, -1, matrixCache); // Corresponding weight matrix const SquareMatrix3& reducedVtxWeight = matrixCache.newVertexWeight; @@ -82,74 +74,26 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex& track, double chi2 = posDiff.dot(reducedVtxWeight * posDiff) + smParams.dot(trkParamWeight * smParams); - // Not yet 4d ready. This can be removed together will all head<> statements, - // once time is consistently introduced to vertexing - ActsMatrix<4, 3> newFullTrkCov(ActsMatrix<4, 3>::Zero()); - newFullTrkCov.block<3, 3>(0, 0) = newTrkCov; - - SquareMatrix4 vtxFullWeight(SquareMatrix4::Zero()); - vtxFullWeight.block<3, 3>(0, 0) = vtxWeight; - - SquareMatrix4 vtxFullCov(SquareMatrix4::Zero()); - vtxFullCov.block<3, 3>(0, 0) = vtxCov; - - Acts::BoundMatrix fullPerTrackCov = detail::createFullTrackCovariance( - sMat, newFullTrkCov, vtxFullWeight, vtxFullCov, newTrkParams); + // Cross covariance matrix between the vertex position and the refitted track + // momentum + ActsMatrix<4, 3> crossCovVP = ActsMatrix<4, 3>::Zero(); + const SquareMatrix3 spatialCrossCovVP = + -vtxCov * posJac.transpose() * trkParamWeight * momJac * sMat; + crossCovVP.topLeftCorner<3, 3>() = spatialCrossCovVP; - // Create new refitted parameters - std::shared_ptr perigeeSurface = - Surface::makeShared(vtx.position()); + // Covariance matrix of the refitted track momentum + const SquareMatrix3 covP = + sMat + (spatialCrossCovVP).transpose() * + (vtxWeight.block<3, 3>(0, 0) * spatialCrossCovVP); - BoundTrackParameters refittedPerigee = BoundTrackParameters( - perigeeSurface, newTrkParams, std::move(fullPerTrackCov), - track.fittedParams.particleHypothesis()); + // Struct comprising the refitted track momentum and the corresponding + // (cross-) covariances + FittedMomentum fittedMom(newTrkMomentum, covP, crossCovVP); // Set new properties - track.fittedParams = refittedPerigee; - track.chi2Track = chi2; - track.ndf = 2 * track.trackWeight; + track.fittedMomentum = fittedMom; + track.chi2 = chi2; + track.ndf = 2 * track.weight; return; } - -inline Acts::BoundMatrix -Acts::KalmanVertexTrackUpdater::detail::createFullTrackCovariance( - const SquareMatrix3& sMat, const ActsMatrix<4, 3>& newTrkCov, - const SquareMatrix4& vtxWeight, const SquareMatrix4& vtxCov, - const BoundVector& newTrkParams) { - // Now new momentum covariance - ActsSquareMatrix<3> momCov = - sMat + (newTrkCov.block<3, 3>(0, 0)).transpose() * - (vtxWeight.block<3, 3>(0, 0) * newTrkCov.block<3, 3>(0, 0)); - - // Full (x,y,z,phi, theta, q/p) covariance matrix - // To be made 7d again after switching to (x,y,z,phi, theta, q/p, t) - ActsSquareMatrix<6> fullTrkCov(ActsSquareMatrix<6>::Zero()); - - fullTrkCov.block<3, 3>(0, 0) = vtxCov.block<3, 3>(0, 0); - fullTrkCov.block<3, 3>(0, 3) = newTrkCov.block<3, 3>(0, 0); - fullTrkCov.block<3, 3>(3, 0) = (newTrkCov.block<3, 3>(0, 0)).transpose(); - fullTrkCov.block<3, 3>(3, 3) = momCov; - - // Combined track jacobian - ActsMatrix<5, 6> trkJac(ActsMatrix<5, 6>::Zero()); - - // First row - trkJac(0, 0) = -std::sin(newTrkParams[2]); - trkJac(0, 1) = std::cos(newTrkParams[2]); - - double tanTheta = std::tan(newTrkParams[3]); - - // Second row - trkJac(1, 0) = -trkJac(0, 1) / tanTheta; - trkJac(1, 1) = trkJac(0, 0) / tanTheta; - - trkJac.block<4, 4>(1, 2) = ActsMatrix<4, 4>::Identity(); - - // Full perigee track covariance - BoundMatrix fullPerTrackCov(BoundMatrix::Identity()); - fullPerTrackCov.block<5, 5>(0, 0) = - (trkJac * (fullTrkCov * trkJac.transpose())); - - return fullPerTrackCov; -} diff --git a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp b/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp index 03a2af70710..6de225fd367 100644 --- a/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp +++ b/Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp @@ -19,7 +19,7 @@ void Acts::KalmanVertexUpdater::updateVertexWithTrack( template void Acts::KalmanVertexUpdater::detail::update( Vertex& vtx, TrackAtVertex& trk, int sign) { - double trackWeight = trk.trackWeight; + double trackWeight = trk.weight; MatrixCache matrixCache; @@ -51,12 +51,12 @@ void Acts::KalmanVertexUpdater::detail::update( // Otherwise just adds track to existing list of tracks at vertex if (sign > 0) { // Update track - trk.chi2Track = trkChi2; + trk.chi2 = trkChi2; trk.ndf = 2 * trackWeight; } // Remove trk from current vertex if (sign < 0) { - trk.trackWeight = 0; + trk.weight = 0; } } diff --git a/Core/include/Acts/Vertexing/TrackAtVertex.hpp b/Core/include/Acts/Vertexing/TrackAtVertex.hpp index 42f88df8bee..605ff02bf20 100644 --- a/Core/include/Acts/Vertexing/TrackAtVertex.hpp +++ b/Core/include/Acts/Vertexing/TrackAtVertex.hpp @@ -1,6 +1,6 @@ // This file is part of the Acts project. // -// Copyright (C) 2019 CERN for the benefit of the Acts project +// Copyright (C) 2019-2023 CERN for the benefit of the Acts project // // This Source Code Form is subject to the terms of the Mozilla Public // License, v. 2.0. If a copy of the MPL was not distributed with this @@ -15,49 +15,86 @@ namespace Acts { -/// @class TrackAtVertex +/// @struct FittedMomentum +/// +/// @brief Vertex fitters return a vertex position and updated track momenta at said position +/// (+corresponding covariances). The updated track momenta and their +/// covariances are saved in the following struct. +struct FittedMomentum { + FittedMomentum(const Vector3& mom, const ActsSquareMatrix<3>& cov, + const ActsMatrix<4, 3>& crossCov) + : momentum(mom), covariance(cov), crossCovariance(crossCov) {} + + Vector3 momentum; + // covariance(i, j) = Cov(p_i, p_j), where p_i is the i-th component of the + // momentum + ActsSquareMatrix<3> covariance; + // crossCovariance(i, j) = Cov(v_i, p_j), where v_i is the i-th component of + // the vertex and p_j is the j-th component of the momentum. If we perform the + // vertex fit in 3D, the last row of crossCovariance is set to 0. + ActsMatrix<4, 3> crossCovariance; +}; + +/// @struct TrackAtVertex /// /// @brief Defines a track at vertex object /// /// @tparam input_track_t Track object type - template struct TrackAtVertex { /// Deleted default constructor TrackAtVertex() = delete; - /// @brief Parameterized constructor + /// @brief Constructor used before the vertex fit (i.e., when we don't know the fitted momentum yet) + /// + /// @param chi2Track Chi2 of the track + /// @param originalTrack Original perigee parameter + TrackAtVertex(const input_track_t* originalTrack, double chi2Track) + : originalParams(originalTrack), chi2(chi2Track) {} + + /// @brief Constructor used before the vertex fit (i.e., when we don't know the fitted momentum yet) with default chi2 /// - /// @param chi2perTrack Chi2 of track - /// @param paramsAtVertex Fitted perigee parameter /// @param originalTrack Original perigee parameter - TrackAtVertex(double chi2perTrack, const BoundTrackParameters& paramsAtVertex, - const input_track_t* originalTrack) - : fittedParams(paramsAtVertex), - originalParams(originalTrack), - chi2Track(chi2perTrack) {} + TrackAtVertex(const input_track_t* originalTrack) + : originalParams(originalTrack) {} - /// @brief Constructor with default chi2 + /// @brief Constructor used when we know the momentum after the fit /// - /// @param paramsAtVertex Fitted perigee parameter + /// @param chi2Track Chi2 of the track + /// @param fittedMom updated momentum after the vertex fit /// @param originalTrack Original perigee parameter - TrackAtVertex(const BoundTrackParameters& paramsAtVertex, - const input_track_t* originalTrack) - : fittedParams(paramsAtVertex), originalParams(originalTrack) {} + TrackAtVertex(const input_track_t* originalTrack, + std::optional fittedMom, double chi2Track) + : originalParams(originalTrack), + fittedMomentum(std::move(fittedMom)), + chi2(chi2Track) {} - /// Fitted perigee - BoundTrackParameters fittedParams; + /// @brief Constructor used when we know the momentum after the fit with default chi2 + /// + /// @param fittedMom updated momentum after the vertex fit + /// @param originalTrack Original perigee parameter + TrackAtVertex(const input_track_t* originalTrack, + std::optional fittedMom) + : originalParams(originalTrack), fittedMomentum(std::move(fittedMom)) {} - /// Original input parameters + /// Original track parameters const input_track_t* originalParams; - /// Chi2 of track - double chi2Track = 0; + /// Momentum after the vertex fit + /// The full track parameters after the vertex fit are: + /// d0 = 0, + /// z0 = 0, + /// phi = fittedMomentum.momentum(0), + /// theta = fittedMomentum.momentum(1), + /// qOverP = fittedMomentum.momentum(2). + std::optional fittedMomentum = std::nullopt; + + /// Chi2 of the track + double chi2 = 0; /// Number degrees of freedom - /// Note: Can be different from integer value - /// since annealing can result in effective - /// non-interger values + /// Note: Can be different from integer value since annealing can result in + /// effective non-integer values. double ndf = 0; /// Value of the compatibility of the track to the actual vertex, based @@ -65,7 +102,7 @@ struct TrackAtVertex { double vertexCompatibility = 0; /// Weight of track in fit - double trackWeight = 1; + double weight = 1; /// The linearized state of the track at vertex LinearizedTrack linearizedState; diff --git a/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.cpp b/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.cpp index 99c53c97635..6e111231248 100644 --- a/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.cpp +++ b/Examples/Io/Performance/ActsExamples/Io/Performance/VertexPerformanceWriter.cpp @@ -660,13 +660,11 @@ ActsExamples::ProcessCode ActsExamples::VertexPerformanceWriter::writeT( } // Save track parameters after the vertex fit - const auto paramsAtVtxFitted = propagateToVtx(trk.fittedParams); - if (paramsAtVtxFitted != std::nullopt) { + if (trk.fittedMomentum.has_value()) { Acts::ActsVector<3> recoMomFitted = - paramsAtVtxFitted->parameters().segment(Acts::eBoundPhi, 3); + trk.fittedMomentum->momentum; const Acts::ActsMatrix<3, 3>& momCovFitted = - paramsAtVtxFitted->covariance()->block<3, 3>( - Acts::eBoundPhi, Acts::eBoundPhi); + trk.fittedMomentum->covariance; innerRecoPhiFitted.push_back(recoMomFitted[0]); innerRecoThetaFitted.push_back(recoMomFitted[1]); innerRecoQOverPFitted.push_back(recoMomFitted[2]); @@ -687,9 +685,12 @@ ActsExamples::ProcessCode ActsExamples::VertexPerformanceWriter::writeT( innerPullQOverPFitted.push_back( pull(diffMomFitted[2], momCovFitted(2, 2), "q/p")); - const auto& recoUnitDirFitted = paramsAtVtxFitted->direction(); + const auto& recoUnitDirFitted = Acts::makeDirectionFromPhiTheta( + recoMomFitted[0], recoMomFitted[1]); double overlapFitted = trueUnitDir.dot(recoUnitDirFitted); innerMomOverlapFitted.push_back(overlapFitted); + } else { + ACTS_WARNING("No fitted track momentum found!"); } } } diff --git a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFinderTests.cpp b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFinderTests.cpp index 1b89f8811de..f03fae935aa 100644 --- a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFinderTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFinderTests.cpp @@ -198,8 +198,8 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_finder_test) { CHECK_CLOSE_OR_SMALL(recoVtx.covariance(), expVtx.covariance, relTol, small); BOOST_CHECK_EQUAL(recoVtx.tracks().size(), expVtx.nTracks); - CHECK_CLOSE_OR_SMALL(recoVtx.tracks()[0].trackWeight, expVtx.trk1Weight, - relTol, small); + CHECK_CLOSE_OR_SMALL(recoVtx.tracks()[0].weight, expVtx.trk1Weight, relTol, + small); CHECK_CLOSE_OR_SMALL(recoVtx.tracks()[0].vertexCompatibility, expVtx.trk1Comp, relTol, small); } @@ -360,10 +360,16 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_finder_usertype_test) { CHECK_CLOSE_OR_SMALL(recoVtx.covariance(), expVtx.covariance, relTol, small); BOOST_CHECK_EQUAL(recoVtx.tracks().size(), expVtx.nTracks); +<<<<<<< HEAD + CHECK_CLOSE_ABS(recoVtx.tracks()[0].weight, expVtx.trk1Weight, 0.003); + CHECK_CLOSE_ABS(recoVtx.tracks()[0].vertexCompatibility, expVtx.trk1Comp, + 0.003); +======= CHECK_CLOSE_OR_SMALL(recoVtx.tracks()[0].trackWeight, expVtx.trk1Weight, relTol, small); CHECK_CLOSE_OR_SMALL(recoVtx.tracks()[0].vertexCompatibility, expVtx.trk1Comp, relTol, small); +>>>>>>> 159acb7f2b1c631748c4ea2a37d1ead092bb69e7 } } diff --git a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp index bf1867027dd..e5dd9df8a52 100644 --- a/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/AdaptiveMultiVertexFitterTests.cpp @@ -209,22 +209,18 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_fitter_test) { iTrack++) { // Index of current vertex int vtxIdx = (int)(iTrack / nTracksPerVtx); - state.vtxInfoMap[&(vtxList[vtxIdx])].trackLinks.push_back( - &(allTracks[iTrack])); - state.tracksAtVerticesMap.insert( - std::make_pair(std::make_pair(&(allTracks[iTrack]), &(vtxList[vtxIdx])), - TrackAtVertex( - 1., allTracks[iTrack], &(allTracks[iTrack])))); + state.vtxInfoMap[&vtxList[vtxIdx]].trackLinks.push_back(&allTracks[iTrack]); + state.tracksAtVerticesMap.insert(std::make_pair( + std::make_pair(&allTracks[iTrack], &vtxList[vtxIdx]), + TrackAtVertex(&allTracks[iTrack], 1.))); // Use first track also for second vertex to let vtx1 and vtx2 // share this track if (iTrack == 0) { - state.vtxInfoMap[&(vtxList.at(1))].trackLinks.push_back( - &(allTracks[iTrack])); - state.tracksAtVerticesMap.insert( - std::make_pair(std::make_pair(&(allTracks[iTrack]), &(vtxList.at(1))), - TrackAtVertex( - 1., allTracks[iTrack], &(allTracks[iTrack])))); + state.vtxInfoMap[&vtxList.at(1)].trackLinks.push_back(&allTracks[iTrack]); + state.tracksAtVerticesMap.insert(std::make_pair( + std::make_pair(&allTracks[iTrack], &vtxList.at(1)), + TrackAtVertex(&allTracks[iTrack], 1.))); } } @@ -475,7 +471,7 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_fitter_test_athena) { vtxInfo1.trackLinks.push_back(&trk); state.tracksAtVerticesMap.insert( std::make_pair(std::make_pair(&trk, &vtx1), - TrackAtVertex(1.5, trk, &trk))); + TrackAtVertex(&trk, 1.5))); } // Prepare second vertex @@ -502,7 +498,7 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_fitter_test_athena) { vtxInfo2.trackLinks.push_back(&trk); state.tracksAtVerticesMap.insert( std::make_pair(std::make_pair(&trk, &vtx2), - TrackAtVertex(1.5, trk, &trk))); + TrackAtVertex(&trk, 1.5))); } state.vtxInfoMap[&vtx1] = std::move(vtxInfo1); @@ -532,7 +528,7 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_fitter_test_athena) { for (const auto& trk : trks1) { auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx1Fitted)); - ACTS_DEBUG("\tTrack weight:" << trkAtVtx.trackWeight); + ACTS_DEBUG("\tTrack weight:" << trkAtVtx.weight); } ACTS_DEBUG("Vertex 1, chi2: " << vtx1FQ.first); ACTS_DEBUG("Vertex 1, ndf: " << vtx1FQ.second); @@ -543,7 +539,7 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_fitter_test_athena) { for (const auto& trk : trks2) { auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx2Fitted)); - ACTS_DEBUG("\tTrack weight:" << trkAtVtx.trackWeight); + ACTS_DEBUG("\tTrack weight:" << trkAtVtx.weight); } ACTS_DEBUG("Vertex 2, chi2: " << vtx2FQ.first); ACTS_DEBUG("Vertex 2, ndf: " << vtx2FQ.second); @@ -579,7 +575,7 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_fitter_test_athena) { for (const auto& trk : trks1) { auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx1Fitted)); - CHECK_CLOSE_ABS(trkAtVtx.trackWeight, expVtx1TrkWeights[trkCount], 0.001); + CHECK_CLOSE_ABS(trkAtVtx.weight, expVtx1TrkWeights[trkCount], 0.001); trkCount++; } CHECK_CLOSE_ABS(vtx1FQ.first, expVtx1chi2, 0.001); @@ -592,7 +588,7 @@ BOOST_AUTO_TEST_CASE(adaptive_multi_vertex_fitter_test_athena) { for (const auto& trk : trks2) { auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx2Fitted)); - CHECK_CLOSE_ABS(trkAtVtx.trackWeight, expVtx2TrkWeights[trkCount], 0.001); + CHECK_CLOSE_ABS(trkAtVtx.weight, expVtx2TrkWeights[trkCount], 0.001); trkCount++; } CHECK_CLOSE_ABS(vtx2FQ.first, expVtx2chi2, 0.001); diff --git a/Tests/UnitTests/Core/Vertexing/FullBilloirVertexFitterTests.cpp b/Tests/UnitTests/Core/Vertexing/FullBilloirVertexFitterTests.cpp index 2cb83803e9e..5922e59d3ef 100644 --- a/Tests/UnitTests/Core/Vertexing/FullBilloirVertexFitterTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/FullBilloirVertexFitterTests.cpp @@ -261,8 +261,6 @@ BOOST_AUTO_TEST_CASE(billoir_vertex_fitter_defaulttrack_test) { auto fittedVertex = fitter.fit(trksPtr, lin, vfOpts, vfState).value(); if (!fittedVertex.tracks().empty()) { CHECK_CLOSE_ABS(fittedVertex.position(), trueVertex.head(3), 1_mm); - auto tracksAtVtx = fittedVertex.tracks(); - auto trackAtVtx = tracksAtVtx[0]; CHECK_CLOSE_ABS(fittedVertex.time(), trueVertex[3], 1_ns); } diff --git a/Tests/UnitTests/Core/Vertexing/IterativeVertexFinderTests.cpp b/Tests/UnitTests/Core/Vertexing/IterativeVertexFinderTests.cpp index 3fb399d1d80..ec2f63d0639 100644 --- a/Tests/UnitTests/Core/Vertexing/IterativeVertexFinderTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/IterativeVertexFinderTests.cpp @@ -242,8 +242,7 @@ BOOST_AUTO_TEST_CASE(iterative_finder_test) { tracks.push_back(std::make_unique(params)); - TrackAtVertex trAtVt(0., params, - tracks.back().get()); + TrackAtVertex trAtVt(tracks.back().get(), 0.); tracksAtTrueVtx.push_back(trAtVt); } @@ -460,7 +459,7 @@ BOOST_AUTO_TEST_CASE(iterative_finder_test_user_track_type) { auto params = extractParameters(paramsUT); - TrackAtVertex trAtVt(0., params, tracks.back().get()); + TrackAtVertex trAtVt(tracks.back().get(), 0.); tracksAtTrueVtx.push_back(trAtVt); } @@ -642,7 +641,7 @@ BOOST_AUTO_TEST_CASE(iterative_finder_test_athena_reference) { // CHECK_CLOSE_ABS(recoVtx.position(), expVtx.position, 0.001_mm); // CHECK_CLOSE_ABS(recoVtx.covariance(), expVtx.covariance, 0.001_mm); // BOOST_CHECK_EQUAL(recoVtx.tracks().size(), expVtx.nTracks); - // CHECK_CLOSE_ABS(recoVtx.tracks()[0].trackWeight, expVtx.trk1Weight, + // CHECK_CLOSE_ABS(recoVtx.tracks()[0].weight, expVtx.trk1Weight, // 0.003); CHECK_CLOSE_ABS(recoVtx.tracks()[0].vertexCompatibility, // expVtx.trk1Comp, // 0.003); diff --git a/Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp b/Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp index 141196193d0..e546d2a1acd 100644 --- a/Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp @@ -164,14 +164,11 @@ BOOST_AUTO_TEST_CASE(Kalman_Vertex_TrackUpdater) { .value(); // Create TrackAtVertex - TrackAtVertex trkAtVtx(0., params, ¶ms); + TrackAtVertex trkAtVtx(¶ms, 0.); // Set linearized state of trackAtVertex trkAtVtx.linearizedState = linTrack; - // Copy parameters for later comparison of old and new version - auto fittedParamsCopy = trkAtVtx.fittedParams; - // Create a vertex Vector3 vtxPos(vXYDist(gen), vXYDist(gen), vZDist(gen)); Vertex vtx(vtxPos); @@ -179,27 +176,9 @@ BOOST_AUTO_TEST_CASE(Kalman_Vertex_TrackUpdater) { // Update trkAtVertex with assumption of originating from vtx KalmanVertexTrackUpdater::update(trkAtVtx, vtx); - // The old distance - double oldDistance = - ip3dEst.calculateDistance(geoContext, fittedParamsCopy, vtxPos, state) - .value(); - - // The new distance after update - double newDistance = - ip3dEst - .calculateDistance(geoContext, trkAtVtx.fittedParams, vtxPos, state) - .value(); - if (debug) { - std::cout << "Old distance: " << oldDistance << std::endl; - std::cout << "New distance: " << newDistance << std::endl; - } - - // Parameters should have changed - BOOST_CHECK_NE(fittedParamsCopy, trkAtVtx.fittedParams); - - // After update, track should be closer to the vertex - BOOST_CHECK_LT(newDistance, oldDistance); - + // Momentum should have changed + BOOST_CHECK_NE(paramVec.segment<3>(eBoundPhi), + trkAtVtx.fittedMomentum.value().momentum); } // end for loop } // end test case diff --git a/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp b/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp index 6f3aa39b4fd..224ef739b51 100644 --- a/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp +++ b/Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp @@ -160,7 +160,7 @@ BOOST_AUTO_TEST_CASE(Kalman_Vertex_Updater) { .value(); // Create TrackAtVertex - TrackAtVertex trkAtVtx(0., params, ¶ms); + TrackAtVertex trkAtVtx(¶ms, 0.); // Set linearized state of trackAtVertex trkAtVtx.linearizedState = linTrack;