src/Vertexing/AdaptiveMultiVertexFinder.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Vertexing/AdaptiveMultiVertexFinder.hpp"
0010
0011 #include "Acts/Utilities/AlgebraHelpers.hpp"
0012 #include "Acts/Vertexing/IVertexFinder.hpp"
0013 #include "Acts/Vertexing/VertexingError.hpp"
0014
0015 #include <algorithm>
0016
0017 namespace Acts {
0018
0019 Result<std::vector<Vertex>> AdaptiveMultiVertexFinder::find(
0020     const std::vector<InputTrack>& allTracks,
0021     const VertexingOptions& vertexingOptions,
0022     IVertexFinder::State& anyState) const {
0023   if (allTracks.empty()) {
0024     ACTS_ERROR("Empty track collection handed to find method");
0025     return VertexingError::EmptyInput;
0026   }
0027
0028   State& state = anyState.template as<State>();
0029   IVertexFinder::State& seedFinderState = state.seedFinderState;
0030   VertexFitterState fitterState(*m_cfg.bField,
0031                                 vertexingOptions.magFieldContext);
0032
0033   const std::vector<InputTrack>& origTracks = allTracks;
0034   std::vector<InputTrack> seedTracks = allTracks;
0035   std::vector<std::unique_ptr<Vertex>> allVertices;
0036   std::vector<Vertex*> allVerticesPtr;
0037
0038   int iteration = 0;
0039   std::vector<InputTrack> removedSeedTracks;
0040   while (!seedTracks.empty() && iteration < m_cfg.maxIterations &&
0041          (m_cfg.addSingleTrackVertices || seedTracks.size() >= 2)) {
0042     Vertex currentConstraint = vertexingOptions.constraint;
0043
0044     // Retrieve seed vertex from all remaining seedTracks
0045     auto seedResult = doSeeding(seedTracks, currentConstraint, vertexingOptions,
0046                                 seedFinderState, removedSeedTracks);
0047     if (!seedResult.ok()) {
0048       return seedResult.error();
0049     }
0050     const auto& seedOptional = *seedResult;
0051
0052     if (!seedOptional.has_value()) {
0053       ACTS_DEBUG(
0054           "No seed found anymore. Break and stop primary vertex finding.");
0055       break;
0056     }
0057     const auto& seedVertex = seedOptional.value();
0058
0059     ACTS_DEBUG("Position of vertex candidate after seeding: "
0060                << seedVertex.fullPosition().transpose());
0061
0062     allVertices.push_back(std::make_unique<Vertex>(seedVertex));
0063     Vertex& vtxCandidate = *allVertices.back();
0064     allVerticesPtr.push_back(&vtxCandidate);
0065
0066     // Clear the seed track collection that has been removed in last iteration
0067     // now after seed finding is done
0068     removedSeedTracks.clear();
0069
0070     // Tracks that are used for searching compatible tracks near a vertex
0071     // candidate
0072     std::vector<InputTrack> searchTracks;
0073     if (m_cfg.doRealMultiVertex) {
0074       searchTracks = origTracks;
0075     } else {
0076       searchTracks = seedTracks;
0077     }
0078
0079     auto prepResult = canPrepareVertexForFit(searchTracks, seedTracks,
0080                                              vtxCandidate, currentConstraint,
0081                                              fitterState, vertexingOptions);
0082
0083     if (!prepResult.ok()) {
0084       return prepResult.error();
0085     }
0086     if (!(*prepResult)) {
0087       ACTS_DEBUG(
0088           "Could not prepare for fit. Discarding the vertex candindate.");
0089       allVertices.pop_back();
0090       allVerticesPtr.pop_back();
0091       if (m_cfg.doNotBreakWhileSeeding) {
0092         continue;
0093       } else {
0094         break;
0095       }
0096     }
0097     // Update fitter state with all vertices
0098     fitterState.addVertexToMultiMap(vtxCandidate);
0099
0100     // Perform the fit
0101     auto fitResult = m_cfg.vertexFitter.addVtxToFit(fitterState, vtxCandidate,
0102                                                     vertexingOptions);
0103     if (!fitResult.ok()) {
0104       return fitResult.error();
0105     }
0106     ACTS_DEBUG("Position of vertex candidate after the fit: "
0107                << vtxCandidate.fullPosition().transpose());
0108     // Check if vertex is good vertex
0109     auto [nCompatibleTracks, isGoodVertex] =
0110         checkVertexAndCompatibleTracks(vtxCandidate, seedTracks, fitterState,
0111                                        vertexingOptions.useConstraintInFit);
0112
0113     ACTS_DEBUG("Vertex is good vertex: " << isGoodVertex);
0114     if (nCompatibleTracks > 0) {
0115       removeCompatibleTracksFromSeedTracks(vtxCandidate, seedTracks,
0116                                            fitterState, removedSeedTracks);
0117     } else {
0118       bool removedIncompatibleTrack = removeTrackIfIncompatible(
0119           vtxCandidate, seedTracks, fitterState, removedSeedTracks,
0120           vertexingOptions.geoContext);
0121       if (!removedIncompatibleTrack) {
0122         ACTS_DEBUG(
0123             "Could not remove any further track from seed tracks. Break.");
0124         allVertices.pop_back();
0125         allVerticesPtr.pop_back();
0126         break;
0127       }
0128     }
0129     auto keepNewVertexResult =
0130         keepNewVertex(vtxCandidate, allVerticesPtr, fitterState);
0131     if (!keepNewVertexResult.ok()) {
0132       return keepNewVertexResult.error();
0133     }
0134     bool keepVertex = isGoodVertex && *keepNewVertexResult;
0135     ACTS_DEBUG("New vertex will be saved: " << keepVertex);
0136
0137     // Delete vertex from allVertices list again if it's not kept
0138     if (!keepVertex) {
0139       auto deleteVertexResult =
0140           deleteLastVertex(vtxCandidate, allVertices, allVerticesPtr,
0141                            fitterState, vertexingOptions);
0142       if (!deleteVertexResult.ok()) {
0143         return deleteVertexResult.error();
0144       }
0145     }
0146     iteration++;
0147   }  // end while loop
0148
0149   return getVertexOutputList(allVerticesPtr, fitterState);
0150 }
0151
0152 auto AdaptiveMultiVertexFinder::doSeeding(
0153     const std::vector<InputTrack>& trackVector, Vertex& currentConstraint,
0154     const VertexingOptions& vertexingOptions,
0155     IVertexFinder::State& seedFinderState,
0156     const std::vector<InputTrack>& removedSeedTracks) const
0157     -> Result<std::optional<Vertex>> {
0158   VertexingOptions seedOptions = vertexingOptions;
0159   seedOptions.constraint = currentConstraint;
0160
0161   m_cfg.seedFinder->setTracksToRemove(seedFinderState, removedSeedTracks);
0162
0163   // Run seed finder
0164   auto seedResult =
0165       m_cfg.seedFinder->find(trackVector, seedOptions, seedFinderState);
0166
0167   if (!seedResult.ok()) {
0168     return seedResult.error();
0169   }
0170   const auto& seedVector = *seedResult;
0171
0172   ACTS_DEBUG("Found " << seedVector.size() << " seeds");
0173
0174   if (seedVector.empty()) {
0175     return std::nullopt;
0176   }
0177   Vertex seedVertex = seedVector.back();
0178
0179   // Update constraints according to seed vertex
0180   setConstraintAfterSeeding(currentConstraint, seedOptions.useConstraintInFit,
0181                             seedVertex);
0182
0183   return seedVertex;
0184 }
0185
0186 void AdaptiveMultiVertexFinder::setConstraintAfterSeeding(
0187     Vertex& currentConstraint, bool useVertexConstraintInFit,
0188     Vertex& seedVertex) const {
0189   if (useVertexConstraintInFit) {
0190     if (!m_cfg.useSeedConstraint) {
0191       // Set seed vertex constraint to old constraint before seeding
0192       seedVertex.setFullCovariance(currentConstraint.fullCovariance());
0193     } else {
0194       // Use the constraint provided by the seed finder
0195       currentConstraint.setFullPosition(seedVertex.fullPosition());
0196       currentConstraint.setFullCovariance(seedVertex.fullCovariance());
0197     }
0198   } else {
0199     currentConstraint.setFullPosition(seedVertex.fullPosition());
0200     currentConstraint.setFullCovariance(m_cfg.initialVariances.asDiagonal());
0201     currentConstraint.setFitQuality(m_cfg.defaultConstrFitQuality);
0202   }
0203 }
0204
0205 Result<double> AdaptiveMultiVertexFinder::getIPSignificance(
0206     const InputTrack& track, const Vertex& vtx,
0207     const VertexingOptions& vertexingOptions) const {
0208   // TODO: In original implementation the covariance of the given vertex is set
0209   // to zero. I did the same here now, but consider removing this and just
0210   // passing the vtx object to the estimator without changing its covariance.
0211   // After all, the vertex seed does have a non-zero convariance in general and
0212   // it probably should be used.
0213   Vertex newVtx = vtx;
0214   if (!m_cfg.useVertexCovForIPEstimation) {
0215     newVtx.setFullCovariance(SquareMatrix4::Zero());
0216   }
0217
0218   auto estRes = m_cfg.ipEstimator.getImpactParameters(
0219       m_cfg.extractParameters(track), newVtx, vertexingOptions.geoContext,
0220       vertexingOptions.magFieldContext, m_cfg.useTime);
0221   if (!estRes.ok()) {
0222     return estRes.error();
0223   }
0224
0225   ImpactParametersAndSigma ipas = *estRes;
0226
0227   // TODO: throw error when encountering negative standard deviations
0228   double chi2Time = 0;
0229   if (m_cfg.useTime) {
0230     if (ipas.sigmaDeltaT.value() > 0) {
0231       chi2Time = std::pow(ipas.deltaT.value() / ipas.sigmaDeltaT.value(), 2);
0232     }
0233   }
0234
0235   double significance = 0.;
0236   if (ipas.sigmaD0 > 0 && ipas.sigmaZ0 > 0) {
0237     significance = std::sqrt(std::pow(ipas.d0 / ipas.sigmaD0, 2) +
0238                              std::pow(ipas.z0 / ipas.sigmaZ0, 2) + chi2Time);
0239   }
0240
0241   return significance;
0242 }
0243
0244 Result<void> AdaptiveMultiVertexFinder::addCompatibleTracksToVertex(
0245     const std::vector<InputTrack>& tracks, Vertex& vtx,
0246     VertexFitterState& fitterState,
0247     const VertexingOptions& vertexingOptions) const {
0248   for (const auto& trk : tracks) {
0249     auto params = m_cfg.extractParameters(trk);
0250     auto pos = params.position(vertexingOptions.geoContext);
0251     // If track is too far away from vertex, do not consider checking the IP
0252     // significance
0253     if (m_cfg.tracksMaxZinterval < std::abs(pos[eZ] - vtx.position()[eZ])) {
0254       continue;
0255     }
0256     auto sigRes = getIPSignificance(trk, vtx, vertexingOptions);
0257     if (!sigRes.ok()) {
0258       return sigRes.error();
0259     }
0260     double ipSig = *sigRes;
0261     if (ipSig < m_cfg.tracksMaxSignificance) {
0262       // Create TrackAtVertex objects, unique for each (track, vertex) pair
0263       fitterState.tracksAtVerticesMap.emplace(std::make_pair(trk, &vtx),
0264                                               TrackAtVertex(params, trk));
0265
0266       // Add the original track parameters to the list for vtx
0267       fitterState.vtxInfoMap[&vtx].trackLinks.push_back(trk);
0268     }
0269   }
0270   return {};
0271 }
0272
0273 Result<bool> AdaptiveMultiVertexFinder::canRecoverFromNoCompatibleTracks(
0274     const std::vector<InputTrack>& allTracks,
0275     const std::vector<InputTrack>& seedTracks, Vertex& vtx,
0276     const Vertex& currentConstraint, VertexFitterState& fitterState,
0277     const VertexingOptions& vertexingOptions) const {
0278   // Recover from cases where no compatible tracks to vertex
0279   // candidate were found
0280   // TODO: This is for now how it's done in athena... this look a bit
0281   // nasty to me
0282   if (fitterState.vtxInfoMap[&vtx].trackLinks.empty()) {
0283     // Find nearest track to vertex candidate
0284     double smallestDeltaZ = std::numeric_limits<double>::max();
0285     double newZ = 0;
0286     bool nearTrackFound = false;
0287     for (const auto& trk : seedTracks) {
0288       auto pos =
0289           m_cfg.extractParameters(trk).position(vertexingOptions.geoContext);
0290       auto zDistance = std::abs(pos[eZ] - vtx.position()[eZ]);
0291       if (zDistance < smallestDeltaZ) {
0292         smallestDeltaZ = zDistance;
0293         nearTrackFound = true;
0294         newZ = pos[eZ];
0295       }
0296     }
0297     if (nearTrackFound) {
0298       vtx.setFullPosition(Vector4(0., 0., newZ, 0.));
0299
0300       // Update vertex info for current vertex
0301       fitterState.vtxInfoMap[&vtx] =
0302           VertexInfo(currentConstraint, vtx.fullPosition());
0303
0304       // Try to add compatible track with adapted vertex position
0305       auto res = addCompatibleTracksToVertex(allTracks, vtx, fitterState,
0306                                              vertexingOptions);
0307       if (!res.ok()) {
0308         return Result<bool>::failure(res.error());
0309       }
0310
0311       if (fitterState.vtxInfoMap[&vtx].trackLinks.empty()) {
0312         ACTS_DEBUG(
0313             "No tracks near seed were found, while at least one was "
0314             "expected. Break.");
0315         return Result<bool>::success(false);
0316       }
0317
0318     } else {
0319       ACTS_DEBUG("No nearest track to seed found. Break.");
0320       return Result<bool>::success(false);
0321     }
0322   }
0323
0324   return Result<bool>::success(true);
0325 }
0326
0327 Result<bool> AdaptiveMultiVertexFinder::canPrepareVertexForFit(
0328     const std::vector<InputTrack>& allTracks,
0329     const std::vector<InputTrack>& seedTracks, Vertex& vtx,
0330     const Vertex& currentConstraint, VertexFitterState& fitterState,
0331     const VertexingOptions& vertexingOptions) const {
0332   // Add vertex info to fitter state
0333   fitterState.vtxInfoMap[&vtx] =
0334       VertexInfo(currentConstraint, vtx.fullPosition());
0335
0336   // Add all compatible tracks to vertex
0337   auto resComp = addCompatibleTracksToVertex(allTracks, vtx, fitterState,
0338                                              vertexingOptions);
0339   if (!resComp.ok()) {
0340     return Result<bool>::failure(resComp.error());
0341   }
0342
0343   // Try to recover from cases where adding compatible track was not possible
0344   auto resRec = canRecoverFromNoCompatibleTracks(allTracks, seedTracks, vtx,
0345                                                  currentConstraint, fitterState,
0346                                                  vertexingOptions);
0347   if (!resRec.ok()) {
0348     return Result<bool>::failure(resRec.error());
0349   }
0350
0351   return Result<bool>::success(*resRec);
0352 }
0353
0354 std::pair<int, bool> AdaptiveMultiVertexFinder::checkVertexAndCompatibleTracks(
0355     Vertex& vtx, const std::vector<InputTrack>& seedTracks,
0356     VertexFitterState& fitterState, bool useVertexConstraintInFit) const {
0357   bool isGoodVertex = false;
0358   int nCompatibleTracks = 0;
0359   for (const auto& trk : fitterState.vtxInfoMap[&vtx].trackLinks) {
0360     const auto& trkAtVtx =
0361         fitterState.tracksAtVerticesMap.at(std::make_pair(trk, &vtx));
0362     if ((trkAtVtx.vertexCompatibility < m_cfg.maxVertexChi2 &&
0363          m_cfg.useFastCompatibility) ||
0364         (trkAtVtx.trackWeight > m_cfg.minWeight &&
0365          trkAtVtx.chi2Track < m_cfg.maxVertexChi2 &&
0366          !m_cfg.useFastCompatibility)) {
0367       // TODO: Understand why looking for compatible tracks only in seed tracks
0368       // and not also in all tracks
0369       if (rangeContainsValue(seedTracks, trk)) {
0370         nCompatibleTracks++;
0371         ACTS_DEBUG("Compatible track found.");
0372
0373         if (m_cfg.addSingleTrackVertices && useVertexConstraintInFit) {
0374           isGoodVertex = true;
0375           break;
0376         }
0377         if (nCompatibleTracks > 1) {
0378           isGoodVertex = true;
0379           break;
0380         }
0381       }
0382     }
0383   }  // end loop over all tracks at vertex
0384
0385   return {nCompatibleTracks, isGoodVertex};
0386 }
0387
0388 auto AdaptiveMultiVertexFinder::removeCompatibleTracksFromSeedTracks(
0389     Vertex& vtx, std::vector<InputTrack>& seedTracks,
0390     VertexFitterState& fitterState,
0391     std::vector<InputTrack>& removedSeedTracks) const -> void {
0392   for (const auto& trk : fitterState.vtxInfoMap[&vtx].trackLinks) {
0393     const auto& trkAtVtx =
0394         fitterState.tracksAtVerticesMap.at(std::make_pair(trk, &vtx));
0395     if ((trkAtVtx.vertexCompatibility < m_cfg.maxVertexChi2 &&
0396          m_cfg.useFastCompatibility) ||
0397         (trkAtVtx.trackWeight > m_cfg.minWeight &&
0398          trkAtVtx.chi2Track < m_cfg.maxVertexChi2 &&
0399          !m_cfg.useFastCompatibility)) {
0400       // Find and remove track from seedTracks
0401       auto foundSeedIter = std::ranges::find(seedTracks, trk);
0402       if (foundSeedIter != seedTracks.end()) {
0403         seedTracks.erase(foundSeedIter);
0404         removedSeedTracks.push_back(trk);
0405       }
0406     }
0407   }
0408 }
0409
0410 bool AdaptiveMultiVertexFinder::removeTrackIfIncompatible(
0411     Vertex& vtx, std::vector<InputTrack>& seedTracks,
0412     VertexFitterState& fitterState, std::vector<InputTrack>& removedSeedTracks,
0413     const GeometryContext& geoCtx) const {
0414   // Try to find the track with highest compatibility
0415   double maxCompatibility = 0;
0416
0417   auto maxCompSeedIt = seedTracks.end();
0418   std::optional<InputTrack> removedTrack = std::nullopt;
0419   for (const auto& trk : fitterState.vtxInfoMap[&vtx].trackLinks) {
0420     const auto& trkAtVtx =
0421         fitterState.tracksAtVerticesMap.at(std::make_pair(trk, &vtx));
0422     double compatibility = trkAtVtx.vertexCompatibility;
0423     if (compatibility > maxCompatibility) {
0424       // Try to find track in seed tracks
0425       auto foundSeedIter = std::ranges::find(seedTracks, trk);
0426       if (foundSeedIter != seedTracks.end()) {
0427         maxCompatibility = compatibility;
0428         maxCompSeedIt = foundSeedIter;
0429         removedTrack = trk;
0430       }
0431     }
0432   }
0433   if (maxCompSeedIt != seedTracks.end()) {
0434     // Remove track with highest compatibility from seed tracks
0435     seedTracks.erase(maxCompSeedIt);
0436     removedSeedTracks.push_back(removedTrack.value());
0437   } else {
0438     // Could not find any seed with compatibility > 0, use alternative
0439     // method to remove a track from seed tracks: Closest track in z to
0440     // vtx candidate
0441     double smallestDeltaZ = std::numeric_limits<double>::max();
0442     auto smallestDzSeedIter = seedTracks.end();
0443     for (unsigned int i = 0; i < seedTracks.size(); i++) {
0444       auto pos = m_cfg.extractParameters(seedTracks[i]).position(geoCtx);
0445       double zDistance = std::abs(pos[eZ] - vtx.position()[eZ]);
0446       if (zDistance < smallestDeltaZ) {
0447         smallestDeltaZ = zDistance;
0448         smallestDzSeedIter = seedTracks.begin() + i;
0449         removedTrack = seedTracks[i];
0450       }
0451     }
0452     if (smallestDzSeedIter != seedTracks.end()) {
0453       seedTracks.erase(smallestDzSeedIter);
0454       removedSeedTracks.push_back(removedTrack.value());
0455     } else {
0456       ACTS_DEBUG("No track found to remove. Stop vertex finding now.");
0457       return false;
0458     }
0459   }
0460   return true;
0461 }
0462
0463 Result<bool> AdaptiveMultiVertexFinder::keepNewVertex(
0464     Vertex& vtx, const std::vector<Vertex*>& allVertices,
0465     VertexFitterState& fitterState) const {
0466   double contamination = 0.;
0467   double contaminationNum = 0;
0468   double contaminationDeNom = 0;
0469   for (const auto& trk : fitterState.vtxInfoMap[&vtx].trackLinks) {
0470     const auto& trkAtVtx =
0471         fitterState.tracksAtVerticesMap.at(std::make_pair(trk, &vtx));
0472     double trackWeight = trkAtVtx.trackWeight;
0473     contaminationNum += trackWeight * (1. - trackWeight);
0474     // MARK: fpeMaskBegin(FLTUND, 1, #2590)
0475     contaminationDeNom += trackWeight * trackWeight;
0476     // MARK: fpeMaskEnd(FLTUND)
0477   }
0478   if (contaminationDeNom != 0) {
0479     contamination = contaminationNum / contaminationDeNom;
0480   }
0481   if (contamination > m_cfg.maximumVertexContamination) {
0482     return Result<bool>::success(false);
0483   }
0484
0485   auto isMergedVertexResult = isMergedVertex(vtx, allVertices);
0486   if (!isMergedVertexResult.ok()) {
0487     return Result<bool>::failure(isMergedVertexResult.error());
0488   }
0489
0490   if (*isMergedVertexResult) {
0491     return Result<bool>::success(false);
0492   }
0493
0494   return Result<bool>::success(true);
0495 }
0496
0497 Result<bool> AdaptiveMultiVertexFinder::isMergedVertex(
0498     const Vertex& vtx, const std::vector<Vertex*>& allVertices) const {
0499   const Vector4& candidatePos = vtx.fullPosition();
0500   const SquareMatrix4& candidateCov = vtx.fullCovariance();
0501
0502   for (const auto otherVtx : allVertices) {
0503     if (&vtx == otherVtx) {
0504       continue;
0505     }
0506     const Vector4& otherPos = otherVtx->fullPosition();
0507     const SquareMatrix4& otherCov = otherVtx->fullCovariance();
0508
0509     double significance = 0;
0510     if (!m_cfg.doFullSplitting) {
0511       if (m_cfg.useTime) {
0512         const Vector2 deltaZT = otherPos.tail<2>() - candidatePos.tail<2>();
0513         SquareMatrix2 sumCovZT = candidateCov.bottomRightCorner<2, 2>() +
0514                                  otherCov.bottomRightCorner<2, 2>();
0515         auto sumCovZTInverse = safeInverse(sumCovZT);
0516         if (!sumCovZTInverse) {
0517           ACTS_ERROR("Vertex z-t covariance matrix is singular.");
0518           ACTS_ERROR("sumCovZT:\n" << sumCovZT);
0519           return Result<bool>::failure(VertexingError::SingularMatrix);
0520         }
0521         significance = std::sqrt(deltaZT.dot(*sumCovZTInverse * deltaZT));
0522       } else {
0523         const double deltaZPos = otherPos[eZ] - candidatePos[eZ];
0524         const double sumVarZ = otherCov(eZ, eZ) + candidateCov(eZ, eZ);
0525         if (sumVarZ <= 0) {
0526           ACTS_ERROR("Variance of the vertex's z-coordinate is not positive.");
0527           ACTS_ERROR("sumVarZ:\n" << sumVarZ);
0528           return Result<bool>::failure(VertexingError::SingularMatrix);
0529         }
0530         // Use only z significance
0531         significance = std::abs(deltaZPos) / std::sqrt(sumVarZ);
0532       }
0533     } else {
0534       if (m_cfg.useTime) {
0535         // Use 4D information for significance
0536         const Vector4 deltaPos = otherPos - candidatePos;
0537         SquareMatrix4 sumCov = candidateCov + otherCov;
0538         auto sumCovInverse = safeInverse(sumCov);
0539         if (!sumCovInverse) {
0540           ACTS_ERROR("Vertex 4D covariance matrix is singular.");
0541           ACTS_ERROR("sumCov:\n" << sumCov);
0542           return Result<bool>::failure(VertexingError::SingularMatrix);
0543         }
0544         significance = std::sqrt(deltaPos.dot(*sumCovInverse * deltaPos));
0545       } else {
0546         // Use 3D information for significance
0547         const Vector3 deltaPos = otherPos.head<3>() - candidatePos.head<3>();
0548         SquareMatrix3 sumCov =
0549             candidateCov.topLeftCorner<3, 3>() + otherCov.topLeftCorner<3, 3>();
0550         auto sumCovInverse = safeInverse(sumCov);
0551         if (!sumCovInverse) {
0552           ACTS_ERROR("Vertex 3D covariance matrix is singular.");
0553           ACTS_ERROR("sumCov:\n" << sumCov);
0554           return Result<bool>::failure(VertexingError::SingularMatrix);
0555         }
0556         significance = std::sqrt(deltaPos.dot(*sumCovInverse * deltaPos));
0557       }
0558     }
0559     if (significance < 0.) {
0560       ACTS_ERROR(
0561           "Found a negative significance (i.e., a negative chi2) when checking "
0562           "if vertices are merged. This should never happen since the vertex "
0563           "covariance should be positive definite, and thus its inverse "
0564           "should be positive definite as well.");
0565       return Result<bool>::failure(VertexingError::MatrixNotPositiveDefinite);
0566     }
0567     if (significance < m_cfg.maxMergeVertexSignificance) {
0568       return Result<bool>::success(true);
0569     }
0570   }
0571   return Result<bool>::success(false);
0572 }
0573
0574 Result<void> AdaptiveMultiVertexFinder::deleteLastVertex(
0575     Vertex& vtx, std::vector<std::unique_ptr<Vertex>>& allVertices,
0576     std::vector<Vertex*>& allVerticesPtr, VertexFitterState& fitterState,
0577     const VertexingOptions& vertexingOptions) const {
0578   allVertices.pop_back();
0579   allVerticesPtr.pop_back();
0580
0581   // Update fitter state with removed vertex candidate
0582   fitterState.removeVertexFromMultiMap(vtx);
0583   // fitterState.vertexCollection contains all vertices that will be fit. When
0584   // we called addVtxToFit, vtx and all vertices that share tracks with vtx were
0585   // added to vertexCollection. Now, we want to refit the same set of vertices
0586   // excluding vx. Therefore, we remove vtx from vertexCollection.
0587   auto removeResult = fitterState.removeVertexFromCollection(vtx, logger());
0588   if (!removeResult.ok()) {
0589     return removeResult.error();
0590   }
0591
0592   for (auto& [key, value] : fitterState.tracksAtVerticesMap) {
0593     // Delete all linearized tracks for current (bad) vertex
0594     if (key.second == &vtx) {
0595       value.isLinearized = false;
0596     }
0597   }
0598
0599   // If no vertices share tracks with vtx we don't need to refit
0600   if (fitterState.vertexCollection.empty()) {
0601     return {};
0602   }
0603
0604   // Do the fit with removed vertex
0605   auto fitResult = m_cfg.vertexFitter.fit(fitterState, vertexingOptions);
0606   if (!fitResult.ok()) {
0607     return fitResult.error();
0608   }
0609
0610   return {};
0611 }
0612
0613 Result<std::vector<Vertex>> AdaptiveMultiVertexFinder::getVertexOutputList(
0614     const std::vector<Vertex*>& allVerticesPtr,
0615     VertexFitterState& fitterState) const {
0616   std::vector<Vertex> outputVec;
0617   for (auto vtx : allVerticesPtr) {
0618     auto& outVtx = *vtx;
0619     std::vector<TrackAtVertex> tracksAtVtx;
0620     for (const auto& trk : fitterState.vtxInfoMap[vtx].trackLinks) {
0621       tracksAtVtx.push_back(
0622           fitterState.tracksAtVerticesMap.at(std::make_pair(trk, vtx)));
0623     }
0624     outVtx.setTracksAtVertex(std::move(tracksAtVtx));
0625     outputVec.push_back(outVtx);
0626   }
0627   return Result<std::vector<Vertex>>(outputVec);
0628 }
0629 }  // namespace Acts