ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
KalmanVertexTrackUpdater.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file KalmanVertexTrackUpdater.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 
13 
14 template <typename input_track_t>
17  const Vertex<input_track_t>& vtx) {
18  const auto vtxPos = vtx.fullPosition().template head<3>();
19 
20  // Get the linearized track
21  const LinearizedTrack& linTrack = track.linearizedState;
22 
23  // Check if linearized state exists
24  if (linTrack.covarianceAtPCA.determinant() == 0.) {
25  // Track has no linearized state, returning w/o update
26  return;
27  }
28 
29  // Retrieve linTrack information
30  const auto posJac = linTrack.positionJacobian.block<5, 3>(0, 0);
31  const auto momJac = linTrack.momentumJacobian.block<5, 3>(0, 0);
32  const auto trkParams = linTrack.parametersAtPCA.head<5>();
33  const auto trkParamWeight = linTrack.weightAtPCA.block<5, 5>(0, 0);
34 
35  // Calculate S matrix
36  ActsSymMatrixD<3> sMat =
37  (momJac.transpose() * (trkParamWeight * momJac)).inverse();
38 
39  const auto residual = linTrack.constantTerm.head<5>();
40 
41  // Refit track momentum
42  Vector3D newTrkMomentum = sMat * momJac.transpose() * trkParamWeight *
43  (trkParams - residual - posJac * vtxPos);
44 
45  // Refit track parameters
46  BoundVector newTrkParams(BoundVector::Zero());
47 
48  // Get phi and theta and correct for possible periodicity changes
49  auto correctedPhiTheta =
50  Acts::detail::ensureThetaBounds(newTrkMomentum(0), newTrkMomentum(1));
51 
52  newTrkParams(ParID_t::ePHI) = correctedPhiTheta.first; // phi
53  newTrkParams(ParID_t::eTHETA) = correctedPhiTheta.second; // theta
54  newTrkParams(ParID_t::eQOP) = newTrkMomentum(2); // qOverP
55 
56  // Vertex covariance and weight matrices
57  const auto vtxCov = vtx.fullCovariance().template block<3, 3>(0, 0);
58  const auto vtxWeight = vtxCov.inverse();
59 
60  // New track covariance matrix
61  auto newTrkCov =
62  -vtxCov * posJac.transpose() * trkParamWeight * momJac * sMat;
63 
65 
66  // Now determine the smoothed chi2 of the track in the following
67  KalmanVertexUpdater::updatePosition<input_track_t>(
68  vtx, linTrack, track.trackWeight, -1, matrixCache);
69 
70  // Corresponding weight matrix
71  const auto& reducedVtxWeight = matrixCache.newVertexWeight;
72 
73  // Difference in positions
74  auto posDiff = vtx.position() - matrixCache.newVertexPos;
75 
76  // Get smoothed params
77  auto smParams =
78  trkParams - (residual + posJac * vtx.fullPosition().template head<3>() +
79  momJac * newTrkMomentum);
80 
81  // New chi2 to be set later
82  double chi2 = posDiff.dot(reducedVtxWeight * posDiff) +
83  smParams.dot(trkParamWeight * smParams);
84 
85  // Not yet 4d ready. This can be removed together will all head<> statements,
86  // once time is consistently introduced to vertexing
87  ActsMatrixD<eSpacePointSize, 3> newFullTrkCov(
89  newFullTrkCov.block<3, 3>(0, 0) = newTrkCov;
90 
91  SpacePointSymMatrix vtxFullWeight(SpacePointSymMatrix::Zero());
92  vtxFullWeight.block<3, 3>(0, 0) = vtxWeight;
93 
94  SpacePointSymMatrix vtxFullCov(SpacePointSymMatrix::Zero());
95  vtxFullCov.block<3, 3>(0, 0) = vtxCov;
96 
97  const auto fullPerTrackCov = detail::createFullTrackCovariance(
98  sMat, newFullTrkCov, vtxFullWeight, vtxFullCov, newTrkParams);
99 
100  // Create new refitted parameters
101  std::shared_ptr<PerigeeSurface> perigeeSurface =
102  Surface::makeShared<PerigeeSurface>(
103  VectorHelpers::position(vtx.fullPosition()));
104 
105  BoundParameters refittedPerigee = BoundParameters(
106  gctx, std::move(fullPerTrackCov), newTrkParams, perigeeSurface);
107 
108  // Set new properties
109  track.fittedParams = refittedPerigee;
110  track.chi2Track = chi2;
111  track.ndf = 2 * track.trackWeight;
112 
113  return;
114 }
115 
118  const ActsSymMatrixD<3>& sMat,
119  const ActsMatrixD<eSpacePointSize, 3>& newTrkCov,
120  const SpacePointSymMatrix& vtxWeight, const SpacePointSymMatrix& vtxCov,
121  const BoundVector& newTrkParams) {
122  // Now new momentum covariance
123  ActsSymMatrixD<3> momCov =
124  sMat + (newTrkCov.block<3, 3>(0, 0)).transpose() *
125  (vtxWeight.block<3, 3>(0, 0) * newTrkCov.block<3, 3>(0, 0));
126 
127  // Full (x,y,z,phi, theta, q/p) covariance matrix
128  // To be made 7d again after switching to (x,y,z,phi, theta, q/p, t)
130 
131  fullTrkCov.block<3, 3>(0, 0) = vtxCov.block<3, 3>(0, 0);
132  fullTrkCov.block<3, 3>(0, 3) = newTrkCov.block<3, 3>(0, 0);
133  fullTrkCov.block<3, 3>(3, 0) = (newTrkCov.block<3, 3>(0, 0)).transpose();
134  fullTrkCov.block<3, 3>(3, 3) = momCov;
135 
136  // Combined track jacobian
138 
139  // First row
140  trkJac(0, 0) = -std::sin(newTrkParams[2]);
141  trkJac(0, 1) = std::cos(newTrkParams[2]);
142 
143  double tanTheta = std::tan(newTrkParams[3]);
144 
145  // Second row
146  trkJac(1, 0) = -trkJac(0, 1) / tanTheta;
147  trkJac(1, 1) = trkJac(0, 0) / tanTheta;
148 
149  trkJac.block<4, 4>(1, 2) = ActsSymMatrixD<4>::Identity();
150 
151  // Full perigee track covariance
152  BoundMatrix fullPerTrackCov(BoundMatrix::Identity());
153  fullPerTrackCov.block<5, 5>(0, 0) =
154  (trkJac * (fullTrkCov * trkJac.transpose()));
155 
156  return fullPerTrackCov;
157 }