Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GaussianTrackDensity.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GaussianTrackDensity.ipp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2020 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
10 
11 #include <math.h>
12 
13 namespace Acts {
14 
15 template <typename input_track_t>
16 std::pair<double, double>
18  State& state, const std::vector<const input_track_t*>& trackList,
19  const std::function<BoundTrackParameters(input_track_t)>& extractParameters)
20  const {
21  auto result = addTracks(state, trackList, extractParameters);
22  if (not result.ok()) {
23  return std::make_pair(0., 0.);
24  }
25 
26  double maxPosition = 0.;
27  double maxDensity = 0.;
28  double maxSecondDerivative = 0.;
29 
30  for (const auto& track : state.trackEntries) {
31  double trialZ = track.z;
32 
33  auto [density, firstDerivative, secondDerivative] =
34  trackDensityAndDerivatives(state, trialZ);
35  if (secondDerivative >= 0. || density <= 0.) {
36  continue;
37  }
38  std::tie(maxPosition, maxDensity, maxSecondDerivative) =
39  updateMaximum(trialZ, density, secondDerivative, maxPosition,
40  maxDensity, maxSecondDerivative);
41 
42  trialZ += stepSize(density, firstDerivative, secondDerivative);
43  std::tie(density, firstDerivative, secondDerivative) =
44  trackDensityAndDerivatives(state, trialZ);
45 
46  if (secondDerivative >= 0. || density <= 0.) {
47  continue;
48  }
49  std::tie(maxPosition, maxDensity, maxSecondDerivative) =
50  updateMaximum(trialZ, density, secondDerivative, maxPosition,
51  maxDensity, maxSecondDerivative);
52  trialZ += stepSize(density, firstDerivative, secondDerivative);
53  std::tie(density, firstDerivative, secondDerivative) =
54  trackDensityAndDerivatives(state, trialZ);
55  if (secondDerivative >= 0. || density <= 0.) {
56  continue;
57  }
58  std::tie(maxPosition, maxDensity, maxSecondDerivative) =
59  updateMaximum(trialZ, density, secondDerivative, maxPosition,
60  maxDensity, maxSecondDerivative);
61  }
62 
63  return (maxSecondDerivative == 0.)
64  ? std::make_pair(0., 0.)
65  : std::make_pair(maxPosition,
66  std::sqrt(-(maxDensity / maxSecondDerivative)));
67 }
68 
69 template <typename input_track_t>
71  State& state, const std::vector<const input_track_t*>& trackList,
72  const std::function<BoundTrackParameters(input_track_t)>& extractParameters)
73  const {
74  return globalMaximumWithWidth(state, trackList, extractParameters).first;
75 }
76 
77 template <typename input_track_t>
79  State& state, const std::vector<const input_track_t*>& trackList,
80  const std::function<BoundTrackParameters(input_track_t)>& extractParameters)
81  const {
82  for (auto trk : trackList) {
83  const BoundTrackParameters& boundParams = extractParameters(*trk);
84  // Get required track parameters
85  const double d0 = boundParams.parameters()[BoundIndices::eBoundLoc0];
86  const double z0 = boundParams.parameters()[BoundIndices::eBoundLoc1];
87  // Get track covariance
88  if (not boundParams.covariance().has_value()) {
89  return VertexingError::NoCovariance;
90  }
91  const auto perigeeCov = *(boundParams.covariance());
92  const double covDD =
94  const double covZZ =
96  const double covDZ =
98  const double covDeterminant = (perigeeCov.block<2, 2>(0, 0)).determinant();
99 
100  // Do track selection based on track cov matrix and m_cfg.d0SignificanceCut
101  if ((covDD <= 0) || (d0 * d0 / covDD > m_cfg.d0SignificanceCut) ||
102  (covZZ <= 0) || (covDeterminant <= 0)) {
103  continue;
104  }
105 
106  // Calculate track density quantities
107  double constantTerm =
108  -(d0 * d0 * covZZ + z0 * z0 * covDD + 2. * d0 * z0 * covDZ) /
109  (2. * covDeterminant);
110  const double linearTerm =
111  (d0 * covDZ + z0 * covDD) /
112  covDeterminant; // minus signs and factors of 2 cancel...
113  const double quadraticTerm = -covDD / (2. * covDeterminant);
114  double discriminant =
115  linearTerm * linearTerm -
116  4. * quadraticTerm * (constantTerm + 2. * m_cfg.z0SignificanceCut);
117  if (discriminant < 0) {
118  continue;
119  }
120 
121  // Add the track to the current maps in the state
122  discriminant = std::sqrt(discriminant);
123  const double zMax = (-linearTerm - discriminant) / (2. * quadraticTerm);
124  const double zMin = (-linearTerm + discriminant) / (2. * quadraticTerm);
125  constantTerm -= std::log(2. * M_PI * std::sqrt(covDeterminant));
126 
127  state.trackEntries.emplace_back(z0, constantTerm, linearTerm, quadraticTerm,
128  zMin, zMax);
129  }
130  return Result<void>::success();
131 }
132 
133 template <typename input_track_t>
134 std::tuple<double, double, double>
136  State& state, double z) const {
137  GaussianTrackDensityStore densityResult(z);
138  for (const auto& trackEntry : state.trackEntries) {
139  densityResult.addTrackToDensity(trackEntry);
140  }
141  return densityResult.densityAndDerivatives();
142 }
143 
144 template <typename input_track_t>
145 std::tuple<double, double, double>
147  double newZ, double newValue, double newSecondDerivative, double maxZ,
148  double maxValue, double maxSecondDerivative) const {
149  if (newValue > maxValue) {
150  maxZ = newZ;
151  maxValue = newValue;
152  maxSecondDerivative = newSecondDerivative;
153  }
154  return {maxZ, maxValue, maxSecondDerivative};
155 }
156 
157 template <typename input_track_t>
159  double ddy) const {
160  return (m_cfg.isGaussianShaped ? (y * dy) / (dy * dy - y * ddy) : -dy / ddy);
161 }
162 
163 template <typename input_track_t>
166  // Take track only if it's within bounds
167  if (entry.lowerBound < m_z && m_z < entry.upperBound) {
168  double delta = std::exp(entry.c0 + m_z * (entry.c1 + m_z * entry.c2));
169  double qPrime = entry.c1 + 2. * m_z * entry.c2;
170  double deltaPrime = delta * qPrime;
171  m_density += delta;
172  m_firstDerivative += deltaPrime;
173  m_secondDerivative += 2. * entry.c2 * delta + qPrime * deltaPrime;
174  }
175 }
176 
177 } // namespace Acts