refactor!: Fitted params to fitted momenta

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

Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.ipp

@@ -247,7 +247,7 @@ auto Acts::AdaptiveMultiVertexFinder<vfitter_t, sfinder_t>::
     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<vfitter_t, sfinder_t>::
         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<vfitter_t, sfinder_t>::
         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<vfitter_t, sfinder_t>::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;
   }

Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp

@@ -287,12 +287,12 @@ Acts::Result<void> 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<input_track_t>(trkAtVtx, *vtx);
       }
     }

Core/include/Acts/Vertexing/FullBilloirVertexFitter.ipp

@@ -214,11 +214,13 @@ Acts::FullBilloirVertexFitter<input_track_t, linearizer_t>::fit(
     Vector4 deltaV = covV * deltaVFac;
     //--------------------------------------------------------------------------------------
     // start momentum related calculations
-
-    std::vector<std::optional<BoundSquareMatrix>> covDeltaP(nTracks);
+    // Cross-covariance matrix of the vertex and the momentum
+    std::vector<ActsMatrix<4, 3>> crossCovVP(nTracks);
+    // Covariance matrix of the momentum
+    std::vector<ActsSquareMatrix<3>> 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<input_track_t, linearizer_t>::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<eBoundSize, 7> 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<input_track_t, linearizer_t>::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<input_track_t> trackAtVertex(
-            billoirTrack.chi2, refittedParams, billoirTrack.originalTrack);
+        // new refitted track momentum
+        FittedMomentum fittedMomentum(trackMomenta[iTrack], covP[iTrack],
+                                      crossCovVP[iTrack]);
+        TrackAtVertex<input_track_t> trackAtVertex(billoirTrack.originalTrack,
+                                                   std::move(fittedMomentum),
+                                                   billoirTrack.chi2);
         tracksAtVertex.push_back(std::move(trackAtVertex));
       }
 

Core/include/Acts/Vertexing/IterativeVertexFinder.ipp

@@ -260,7 +260,7 @@ Acts::IterativeVertexFinder<vfitter_t, sfinder_t>::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<vfitter_t, sfinder_t>::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<vfitter_t, sfinder_t>::countSignificantTracks(
     const Vertex<InputTrack_t>& vtx) const {
   return std::count_if(vtx.tracks().begin(), vtx.tracks().end(),
                        [this](TrackAtVertex<InputTrack_t> trk) {
-                         return trk.trackWeight > m_cfg.cutOffTrackWeight;
+                         return trk.weight > m_cfg.cutOffTrackWeight;
                        });
 }

Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.hpp

@@ -28,24 +28,6 @@ template <typename input_track_t>
 void update(TrackAtVertex<input_track_t>& track,
             const Vertex<input_track_t>& 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
 

Core/include/Acts/Vertexing/KalmanVertexTrackUpdater.ipp

@@ -43,29 +43,21 @@ void Acts::KalmanVertexTrackUpdater::update(TrackAtVertex<input_track_t>& 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<input_track_t>(
-      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<input_track_t>& 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> perigeeSurface =
-      Surface::makeShared<PerigeeSurface>(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;
-}

Core/include/Acts/Vertexing/KalmanVertexUpdater.ipp

@@ -19,7 +19,7 @@ void Acts::KalmanVertexUpdater::updateVertexWithTrack(
 template <typename input_track_t>
 void Acts::KalmanVertexUpdater::detail::update(
     Vertex<input_track_t>& vtx, TrackAtVertex<input_track_t>& 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;
   }
 }