src/Vertexing/AdaptiveMultiVertexFinder.cpp

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