ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
KalmanVertexUpdater.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file KalmanVertexUpdater.ipp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019 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 
9 #include <algorithm>
11 
12 template <typename input_track_t>
15  detail::update<input_track_t>(vtx, trk, 1);
16 }
17 
18 template <typename input_track_t>
21  double trackWeight = trk.trackWeight;
22 
23  MatrixCache matrixCache;
24 
25  updatePosition(vtx, trk.linearizedState, trackWeight, sign, matrixCache);
26 
27  // Get fit quality parameters wrt to old vertex
28  std::pair fitQuality = vtx.fitQuality();
29  double chi2 = fitQuality.first;
30  double ndf = fitQuality.second;
31 
32  // Chi2 wrt to track parameters
33  double trkChi2 = detail::trackParametersChi2<input_track_t>(
34  trk.linearizedState, matrixCache);
35 
36  // Calculate new chi2
37  chi2 += sign * (detail::vertexPositionChi2<input_track_t>(vtx, matrixCache) +
38  trackWeight * trkChi2);
39 
40  // Calculate ndf
41  ndf += sign * trackWeight * 2.;
42 
43  // Updating the vertex
44  vtx.setPosition(matrixCache.newVertexPos);
45  vtx.setCovariance(matrixCache.newVertexCov);
46  vtx.setFitQuality(chi2, ndf);
47 
48  // Updates track at vertex if already there
49  // by removing it first and adding new one.
50  // Otherwise just adds track to existing list of tracks at vertex
51  if (sign > 0) {
52  // Update track
53  trk.chi2Track = trkChi2;
54  trk.ndf = 2 * trackWeight;
55  }
56  // Remove trk from current vertex
57  if (sign < 0) {
58  trk.trackWeight = 0;
59  }
60 }
61 
62 template <typename input_track_t>
64  const Acts::Vertex<input_track_t>& vtx,
65  const Acts::LinearizedTrack& linTrack, double trackWeight, int sign,
66  MatrixCache& matrixCache) {
67  // Retrieve linTrack information
68  // TODO: To make 4-D compatible, remove block<> and head<> statements
69  const auto posJac = linTrack.positionJacobian.block<5, 3>(0, 0);
70  const auto momJac =
71  linTrack.momentumJacobian.block<5, 3>(0, 0); // B_k in comments below
72  const auto trkParams = linTrack.parametersAtPCA.head<5>();
73  const auto constTerm = linTrack.constantTerm.head<5>();
74  const auto trkParamWeight = linTrack.weightAtPCA.block<5, 5>(0, 0);
75 
76  // Vertex to be updated
77  const auto& oldVtxPos = vtx.position();
78  matrixCache.oldVertexWeight = (vtx.covariance()).inverse();
79 
80  // W_k matrix
81  matrixCache.momWeightInv =
82  (momJac.transpose() * (trkParamWeight * momJac)).inverse();
83 
84  // G_b = G_k - G_k*B_k*W_k*B_k^(T)*G_k^T
85  auto gBmat = trkParamWeight -
86  trkParamWeight *
87  (momJac * (matrixCache.momWeightInv * momJac.transpose())) *
88  trkParamWeight.transpose();
89 
90  // New vertex cov matrix
91  matrixCache.newVertexWeight =
92  matrixCache.oldVertexWeight +
93  trackWeight * sign * posJac.transpose() * (gBmat * posJac);
94  matrixCache.newVertexCov = matrixCache.newVertexWeight.inverse();
95 
96  // New vertex position
97  matrixCache.newVertexPos =
98  matrixCache.newVertexCov * (matrixCache.oldVertexWeight * oldVtxPos +
99  trackWeight * sign * posJac.transpose() *
100  gBmat * (trkParams - constTerm));
101 }
102 
103 template <typename input_track_t>
105  const Vertex<input_track_t>& oldVtx, const MatrixCache& matrixCache) {
106  Vector3D posDiff = matrixCache.newVertexPos - oldVtx.position();
107 
108  // Calculate and return corresponding chi2
109  return posDiff.transpose() * (matrixCache.oldVertexWeight * posDiff);
110 }
111 
112 template <typename input_track_t>
114  const LinearizedTrack& linTrack, const MatrixCache& matrixCache) {
115  // Track properties
116  const auto posJac = linTrack.positionJacobian.block<5, 3>(0, 0);
117  const auto momJac = linTrack.momentumJacobian.block<5, 3>(0, 0);
118  const auto trkParams = linTrack.parametersAtPCA.head<5>();
119  const auto constTerm = linTrack.constantTerm.head<5>();
120  const auto trkParamWeight = linTrack.weightAtPCA.block<5, 5>(0, 0);
121 
122  const auto jacVtx = posJac * matrixCache.newVertexPos;
123 
124  // Refitted track momentum
125  Vector3D newTrackMomentum = matrixCache.momWeightInv * momJac.transpose() *
126  trkParamWeight * (trkParams - constTerm - jacVtx);
127 
128  // Refitted track parameters
129  auto newTrkParams = constTerm + jacVtx + momJac * newTrackMomentum;
130 
131  // Parameter difference
132  auto paramDiff = trkParams - newTrkParams;
133 
134  // Return chi2
135  return paramDiff.transpose() * (trkParamWeight * paramDiff);
136 }