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