src/Vertexing/GaussianTrackDensity.cpp

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