src/Vertexing/GridDensityVertexFinder.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/GridDensityVertexFinder.hpp"
0010
0011 namespace Acts {
0012
0013 auto GridDensityVertexFinder::find(const std::vector<InputTrack>& trackVector,
0014                                    const VertexingOptions& vertexingOptions,
0015                                    IVertexFinder::State& anyState) const
0016     -> Result<std::vector<Vertex>> {
0017   auto& state = anyState.as<State>();
0018   // Remove density contributions from tracks removed from track collection
0019   if (m_cfg.cacheGridStateForTrackRemoval && state.isInitialized &&
0020       !state.tracksToRemove.empty()) {
0021     // Bool to check if removable tracks, that pass selection, still exist
0022     bool couldRemoveTracks = false;
0023     for (auto trk : state.tracksToRemove) {
0024       if (!state.trackSelectionMap.at(trk)) {
0025         // Track was never added to grid, so cannot remove it
0026         continue;
0027       }
0028       couldRemoveTracks = true;
0029       auto binAndTrackGrid = state.binAndTrackGridMap.at(trk);
0030       m_cfg.gridDensity.removeTrackGridFromMainGrid(
0031           binAndTrackGrid.first, binAndTrackGrid.second, state.mainGrid);
0032     }
0033     if (!couldRemoveTracks) {
0034       // No tracks were removed anymore
0035       // Return empty seed
0036       // (Note: Upstream finder should check for this break condition)
0037       return std::vector<Vertex>{};
0038     }
0039   } else {
0040     state.mainGrid =
0041         MainGridVector::Zero(m_cfg.gridDensity.config().mainGridSize);
0042     // Fill with track densities
0043     for (auto trk : trackVector) {
0044       const BoundTrackParameters& trkParams = m_cfg.extractParameters(trk);
0045       // Take only tracks that fulfill selection criteria
0046       if (!doesPassTrackSelection(trkParams)) {
0047         if (m_cfg.cacheGridStateForTrackRemoval) {
0048           state.trackSelectionMap[trk] = false;
0049         }
0050         continue;
0051       }
0052       auto binAndTrackGrid =
0053           m_cfg.gridDensity.addTrack(trkParams, state.mainGrid);
0054       // Cache track density contribution to main grid if enabled
0055       if (m_cfg.cacheGridStateForTrackRemoval) {
0056         state.binAndTrackGridMap[trk] = binAndTrackGrid;
0057         state.trackSelectionMap[trk] = true;
0058       }
0059     }
0060     state.isInitialized = true;
0061   }
0062
0063   double z = 0;
0064   double width = 0;
0065   if (!state.mainGrid.isZero()) {
0066     if (!m_cfg.estimateSeedWidth) {
0067       // Get z value of highest density bin
0068       auto maxZres = m_cfg.gridDensity.getMaxZPosition(state.mainGrid);
0069
0070       if (!maxZres.ok()) {
0071         return maxZres.error();
0072       }
0073       z = *maxZres;
0074     } else {
0075       // Get z value of highest density bin and width
0076       auto maxZres = m_cfg.gridDensity.getMaxZPositionAndWidth(state.mainGrid);
0077
0078       if (!maxZres.ok()) {
0079         return maxZres.error();
0080       }
0081       z = (*maxZres).first;
0082       width = (*maxZres).second;
0083     }
0084   }
0085
0086   // Construct output vertex
0087   Vector3 seedPos = vertexingOptions.constraint.position() + Vector3(0., 0., z);
0088
0089   Vertex returnVertex = Vertex(seedPos);
0090
0091   SquareMatrix4 seedCov = vertexingOptions.constraint.fullCovariance();
0092
0093   if (width != 0.) {
0094     // Use z-constraint from seed width
0095     seedCov(2, 2) = width * width;
0096   }
0097
0098   returnVertex.setFullCovariance(seedCov);
0099
0100   return std::vector<Vertex>{returnVertex};
0101 }
0102
0103 auto GridDensityVertexFinder::doesPassTrackSelection(
0104     const BoundTrackParameters& trk) const -> bool {
0105   // Get required track parameters
0106   const double d0 = trk.parameters()[BoundIndices::eBoundLoc0];
0107   const double z0 = trk.parameters()[BoundIndices::eBoundLoc1];
0108   // Get track covariance
0109   if (!trk.covariance().has_value()) {
0110     return false;
0111   }
0112   const auto perigeeCov = *(trk.covariance());
0113   const double covDD =
0114       perigeeCov(BoundIndices::eBoundLoc0, BoundIndices::eBoundLoc0);
0115   const double covZZ =
0116       perigeeCov(BoundIndices::eBoundLoc1, BoundIndices::eBoundLoc1);
0117   const double covDZ =
0118       perigeeCov(BoundIndices::eBoundLoc0, BoundIndices::eBoundLoc1);
0119   const double covDeterminant = covDD * covZZ - covDZ * covDZ;
0120
0121   // Do track selection based on track cov matrix and d0SignificanceCut
0122   if ((covDD <= 0) || (d0 * d0 / covDD > m_cfg.d0SignificanceCut) ||
0123       (covZZ <= 0) || (covDeterminant <= 0)) {
0124     return false;
0125   }
0126
0127   // Calculate track density quantities
0128   double constantTerm =
0129       -(d0 * d0 * covZZ + z0 * z0 * covDD + 2. * d0 * z0 * covDZ) /
0130       (2. * covDeterminant);
0131   const double linearTerm = (d0 * covDZ + z0 * covDD) / covDeterminant;
0132   const double quadraticTerm = -covDD / (2. * covDeterminant);
0133   double discriminant =
0134       linearTerm * linearTerm -
0135       4. * quadraticTerm * (constantTerm + 2. * m_cfg.z0SignificanceCut);
0136   if (discriminant < 0) {
0137     return false;
0138   }
0139
0140   return true;
0141 }
0142
0143 }  // namespace Acts