src/Vertexing/GaussianTrackDensity.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/GaussianTrackDensity.hpp"
0010
0011 #include "Acts/Vertexing/VertexingError.hpp"
0012
0013 #include <cmath>
0014 #include <numbers>
0015
0016 namespace Acts {
0017
0018 Result<std::optional<std::pair<double, double>>>
0019 Acts::GaussianTrackDensity::globalMaximumWithWidth(
0020     State& state, const std::vector<InputTrack>& trackList) const {
0021   auto result = addTracks(state, trackList);
0022   if (!result.ok()) {
0023     return result.error();
0024   }
0025
0026   double maxPosition = 0.;
0027   double maxDensity = 0.;
0028   double maxSecondDerivative = 0.;
0029
0030   for (const auto& track : state.trackEntries) {
0031     double trialZ = track.z;
0032
0033     auto [density, firstDerivative, secondDerivative] =
0034         trackDensityAndDerivatives(state, trialZ);
0035     if (secondDerivative >= 0. || density <= 0.) {
0036       continue;
0037     }
0038     std::tie(maxPosition, maxDensity, maxSecondDerivative) =
0039         updateMaximum(trialZ, density, secondDerivative, maxPosition,
0040                       maxDensity, maxSecondDerivative);
0041
0042     trialZ += stepSize(density, firstDerivative, secondDerivative);
0043     std::tie(density, firstDerivative, secondDerivative) =
0044         trackDensityAndDerivatives(state, trialZ);
0045
0046     if (secondDerivative >= 0. || density <= 0.) {
0047       continue;
0048     }
0049     std::tie(maxPosition, maxDensity, maxSecondDerivative) =
0050         updateMaximum(trialZ, density, secondDerivative, maxPosition,
0051                       maxDensity, maxSecondDerivative);
0052     trialZ += stepSize(density, firstDerivative, secondDerivative);
0053     std::tie(density, firstDerivative, secondDerivative) =
0054         trackDensityAndDerivatives(state, trialZ);
0055     if (secondDerivative >= 0. || density <= 0.) {
0056       continue;
0057     }
0058     std::tie(maxPosition, maxDensity, maxSecondDerivative) =
0059         updateMaximum(trialZ, density, secondDerivative, maxPosition,
0060                       maxDensity, maxSecondDerivative);
0061   }
0062
0063   if (maxSecondDerivative == 0.) {
0064     return std::nullopt;
0065   }
0066
0067   return std::pair{maxPosition, std::sqrt(-(maxDensity / maxSecondDerivative))};
0068 }
0069
0070 Result<std::optional<double>> Acts::GaussianTrackDensity::globalMaximum(
0071     State& state, const std::vector<InputTrack>& trackList) const {
0072   auto maxRes = globalMaximumWithWidth(state, trackList);
0073   if (!maxRes.ok()) {
0074     return maxRes.error();
0075   }
0076   const auto& maxOpt = *maxRes;
0077   if (!maxOpt.has_value()) {
0078     return std::nullopt;
0079   }
0080   return maxOpt->first;
0081 }
0082
0083 Result<void> Acts::GaussianTrackDensity::addTracks(
0084     State& state, const std::vector<InputTrack>& trackList) const {
0085   for (auto trk : trackList) {
0086     const BoundTrackParameters& boundParams = m_cfg.extractParameters(trk);
0087     // Get required track parameters
0088     const double d0 = boundParams.parameters()[BoundIndices::eBoundLoc0];
0089     const double z0 = boundParams.parameters()[BoundIndices::eBoundLoc1];
0090     // Get track covariance
0091     if (!boundParams.covariance().has_value()) {
0092       return VertexingError::NoCovariance;
0093     }
0094     const auto perigeeCov = *(boundParams.covariance());
0095     const double covDD =
0096         perigeeCov(BoundIndices::eBoundLoc0, BoundIndices::eBoundLoc0);
0097     const double covZZ =
0098         perigeeCov(BoundIndices::eBoundLoc1, BoundIndices::eBoundLoc1);
0099     const double covDZ =
0100         perigeeCov(BoundIndices::eBoundLoc0, BoundIndices::eBoundLoc1);
0101     const double covDeterminant = (perigeeCov.block<2, 2>(0, 0)).determinant();
0102
0103     // Do track selection based on track cov matrix and m_cfg.d0SignificanceCut
0104     if ((covDD <= 0) || (d0 * d0 / covDD > m_cfg.d0SignificanceCut) ||
0105         (covZZ <= 0) || (covDeterminant <= 0)) {
0106       continue;
0107     }
0108
0109     // Calculate track density quantities
0110     double constantTerm =
0111         -(d0 * d0 * covZZ + z0 * z0 * covDD + 2. * d0 * z0 * covDZ) /
0112         (2. * covDeterminant);
0113     const double linearTerm =
0114         (d0 * covDZ + z0 * covDD) /
0115         covDeterminant;  // minus signs and factors of 2 cancel...
0116     const double quadraticTerm = -covDD / (2. * covDeterminant);
0117     double discriminant =
0118         linearTerm * linearTerm -
0119         4. * quadraticTerm * (constantTerm + 2. * m_cfg.z0SignificanceCut);
0120     if (discriminant < 0) {
0121       continue;
0122     }
0123
0124     // Add the track to the current maps in the state
0125     discriminant = std::sqrt(discriminant);
0126     const double zMax = (-linearTerm - discriminant) / (2. * quadraticTerm);
0127     const double zMin = (-linearTerm + discriminant) / (2. * quadraticTerm);
0128     constantTerm -= std::log(2. * std::numbers::pi * std::sqrt(covDeterminant));
0129
0130     state.trackEntries.emplace_back(z0, constantTerm, linearTerm, quadraticTerm,
0131                                     zMin, zMax);
0132   }
0133   return Result<void>::success();
0134 }
0135
0136 std::tuple<double, double, double>
0137 Acts::GaussianTrackDensity::trackDensityAndDerivatives(State& state,
0138                                                        double z) const {
0139   GaussianTrackDensityStore densityResult(z);
0140   for (const auto& trackEntry : state.trackEntries) {
0141     densityResult.addTrackToDensity(trackEntry);
0142   }
0143   return densityResult.densityAndDerivatives();
0144 }
0145
0146 std::tuple<double, double, double> Acts::GaussianTrackDensity::updateMaximum(
0147     double newZ, double newValue, double newSecondDerivative, double maxZ,
0148     double maxValue, double maxSecondDerivative) const {
0149   if (newValue > maxValue) {
0150     maxZ = newZ;
0151     maxValue = newValue;
0152     maxSecondDerivative = newSecondDerivative;
0153   }
0154   return {maxZ, maxValue, maxSecondDerivative};
0155 }
0156
0157 double Acts::GaussianTrackDensity::stepSize(double y, double dy,
0158                                             double ddy) const {
0159   return (m_cfg.isGaussianShaped ? (y * dy) / (dy * dy - y * ddy) : -dy / ddy);
0160 }
0161
0162 void Acts::GaussianTrackDensity::GaussianTrackDensityStore::addTrackToDensity(
0163     const TrackEntry& entry) {
0164   // Take track only if it's within bounds
0165   if (entry.lowerBound < m_z && m_z < entry.upperBound) {
0166     double delta = std::exp(entry.c0 + m_z * (entry.c1 + m_z * entry.c2));
0167     double qPrime = entry.c1 + 2. * m_z * entry.c2;
0168     double deltaPrime = delta * qPrime;
0169     m_density += delta;
0170     m_firstDerivative += deltaPrime;
0171     m_secondDerivative += 2. * entry.c2 * delta + qPrime * deltaPrime;
0172   }
0173 }
0174
0175 }  // namespace Acts