ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FullBilloirVertexFitter.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FullBilloirVertexFitter.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 
14 
15 namespace {
16 
20 template <typename input_track_t>
21 struct BilloirTrack {
23 
24  BilloirTrack(const input_track_t* params, Acts::LinearizedTrack lTrack)
25  : originalTrack(params), linTrack(std::move(lTrack)) {}
26 
27  BilloirTrack(const BilloirTrack& arg) = default;
28 
29  const input_track_t* originalTrack;
30  Acts::LinearizedTrack linTrack;
31  double chi2;
32  Jacobian DiMat; // position jacobian
33  Acts::ActsMatrixD<Acts::eBoundParametersSize, 3> EiMat; // momentum jacobian
34  Acts::ActsSymMatrixD<3> CiMat; // = EtWmat * Emat (see below)
35  Acts::ActsMatrixD<4, 3> BiMat; // = DiMat^T * Wi * EiMat
36  Acts::ActsSymMatrixD<3> CiInv; // = (EiMat^T * Wi * EiMat)^-1
37  Acts::Vector3D UiVec; // = EiMat^T * Wi * dqi
38  Acts::ActsMatrixD<4, 3> BCiMat; // = BiMat * Ci^-1
39  Acts::BoundVector deltaQ;
40 };
41 
45 struct BilloirVertex {
46  BilloirVertex() = default;
47 
49  Acts::SpacePointSymMatrix::Zero()}; // Amat = sum{DiMat^T * Wi * DiMat}
51  Acts::SpacePointVector::Zero()}; // Tvec = sum{DiMat^T * Wi * dqi}
53  Acts::SpacePointSymMatrix::Zero()}; // BCBmat =
54  // sum{BiMat
55  // * Ci^-1 *
56  // BiMat^T}
58  Acts::SpacePointVector::Zero()}; // BCUvec = sum{BiMat * Ci^-1 * UiVec}
59 };
60 
61 } // end anonymous namespace
62 
63 template <typename input_track_t, typename linearizer_t>
66  const std::vector<const input_track_t*>& paramVector,
67  const linearizer_t& linearizer,
68  const VertexingOptions<input_track_t>& vertexingOptions) const {
69  double chi2 = std::numeric_limits<double>::max();
70  double newChi2 = 0;
71  unsigned int nTracks = paramVector.size();
72 
73  if (nTracks == 0) {
74  return Vertex<input_track_t>(Vector3D(0., 0., 0.));
75  }
76 
77  // Set number of degrees of freedom
78  // ndf = (5-3) * nTracks - 3;
79  int ndf = 2 * nTracks - 3;
80  if (nTracks < 2) {
81  ndf = 1;
82  }
83 
84  // Determine if we do contraint fit or not by checking if an
85  // invertible non-zero constraint vertex covariance is given
86  bool isConstraintFit = false;
87  if (vertexingOptions.vertexConstraint.covariance().determinant() != 0) {
88  isConstraintFit = true;
89  ndf += 3;
90  }
91 
92  std::vector<BilloirTrack<input_track_t>> billoirTracks;
93 
94  std::vector<Vector3D> trackMomenta;
95 
96  SpacePointVector linPoint(vertexingOptions.vertexConstraint.fullPosition());
97 
98  Vertex<input_track_t> fittedVertex;
99 
100  for (int nIter = 0; nIter < m_cfg.maxIterations; ++nIter) {
101  billoirTracks.clear();
102 
103  newChi2 = 0;
104 
105  BilloirVertex billoirVertex;
106  int iTrack = 0;
107  // iterate over all tracks
108  for (const input_track_t* trackContainer : paramVector) {
109  const auto& trackParams = extractParameters(*trackContainer);
110  if (nIter == 0) {
111  double phi = trackParams.parameters()[ParID_t::ePHI];
112  double theta = trackParams.parameters()[ParID_t::eTHETA];
113  double qop = trackParams.parameters()[ParID_t::eQOP];
114  trackMomenta.push_back(Vector3D(phi, theta, qop));
115  }
116 
117  auto result = linearizer.linearizeTrack(trackParams, linPoint,
118  vertexingOptions.geoContext,
119  vertexingOptions.magFieldContext);
120  if (result.ok()) {
121  const auto& linTrack = *result;
122  const auto& parametersAtPCA = linTrack.parametersAtPCA;
123  double d0 = parametersAtPCA[ParID_t::eLOC_D0];
124  double z0 = parametersAtPCA[ParID_t::eLOC_Z0];
125  double phi = parametersAtPCA[ParID_t::ePHI];
126  double theta = parametersAtPCA[ParID_t::eTHETA];
127  double qOverP = parametersAtPCA[ParID_t::eQOP];
128 
129  // calculate f(V_0,p_0) f_d0 = f_z0 = 0
130  double fPhi = trackMomenta[iTrack][0];
131  double fTheta = trackMomenta[iTrack][1];
132  double fQOvP = trackMomenta[iTrack][2];
133  BilloirTrack<input_track_t> currentBilloirTrack(trackContainer,
134  linTrack);
135 
136  currentBilloirTrack.deltaQ << d0, z0, phi - fPhi, theta - fTheta,
137  qOverP - fQOvP, 0;
138 
139  // position jacobian (D matrix)
141  Dmat = linTrack.positionJacobian;
142 
143  // momentum jacobian (E matrix)
145  Emat = linTrack.momentumJacobian;
146  // cache some matrix multiplications
149  BoundSymMatrix Wi = linTrack.weightAtPCA;
150 
151  DtWmat = Dmat.transpose() * Wi;
152  EtWmat = Emat.transpose() * Wi;
153 
154  // compute billoir tracks
155  currentBilloirTrack.DiMat = Dmat;
156  currentBilloirTrack.EiMat = Emat;
157  currentBilloirTrack.CiMat = EtWmat * Emat;
158  currentBilloirTrack.BiMat = DtWmat * Emat; // DiMat^T * Wi * EiMat
159  currentBilloirTrack.UiVec =
160  EtWmat * currentBilloirTrack.deltaQ; // EiMat^T * Wi * dqi
161  currentBilloirTrack.CiInv =
162  (EtWmat * Emat).inverse(); // (EiMat^T * Wi * EiMat)^-1
163 
164  // sum up over all tracks
165  billoirVertex.Tvec +=
166  DtWmat * currentBilloirTrack.deltaQ; // sum{DiMat^T * Wi * dqi}
167  billoirVertex.Amat += DtWmat * Dmat; // sum{DiMat^T * Wi * DiMat}
168 
169  // remember those results for all tracks
170  currentBilloirTrack.BCiMat =
171  currentBilloirTrack.BiMat *
172  currentBilloirTrack.CiInv; // BCi = BiMat * Ci^-1
173 
174  // and some summed results
175  billoirVertex.BCUvec +=
176  currentBilloirTrack.BCiMat *
177  currentBilloirTrack.UiVec; // sum{BiMat * Ci^-1 * UiVec}
178  billoirVertex.BCBmat +=
179  currentBilloirTrack.BCiMat *
180  currentBilloirTrack.BiMat
181  .transpose(); // sum{BiMat * Ci^-1 * BiMat^T}
182 
183  billoirTracks.push_back(currentBilloirTrack);
184  ++iTrack;
185  } else {
186  return result.error();
187  }
188  } // end loop tracks
189 
190  // calculate delta (billoirFrameOrigin-position), might be changed by the
191  // beam-const
192  SpacePointVector Vdel =
193  billoirVertex.Tvec -
194  billoirVertex.BCUvec; // Vdel = Tvec-sum{BiMat*Ci^-1*UiVec}
195  SpacePointSymMatrix VwgtMat =
196  billoirVertex.Amat -
197  billoirVertex.BCBmat; // VwgtMat = Amat-sum{BiMat*Ci^-1*BiMat^T}
198  if (isConstraintFit) {
199  // this will be 0 for first iteration but != 0 from second on
200  SpacePointVector posInBilloirFrame =
201  vertexingOptions.vertexConstraint.fullPosition() - linPoint;
202 
203  Vdel += vertexingOptions.vertexConstraint.fullCovariance().inverse() *
204  posInBilloirFrame;
205  VwgtMat += vertexingOptions.vertexConstraint.fullCovariance().inverse();
206  }
207 
208  // cov(deltaV) = VwgtMat^-1
209  SpacePointSymMatrix covDeltaVmat = VwgtMat.inverse();
210  // deltaV = cov_(deltaV) * Vdel;
211  SpacePointVector deltaV = covDeltaVmat * Vdel;
212  //--------------------------------------------------------------------------------------
213  // start momentum related calculations
214 
215  std::vector<std::optional<BoundSymMatrix>> covDeltaPmat(nTracks);
216 
217  iTrack = 0;
218  for (auto& bTrack : billoirTracks) {
219  Vector3D deltaP =
220  (bTrack.CiInv) * (bTrack.UiVec - bTrack.BiMat.transpose() * deltaV);
221 
222  // update track momenta
223  trackMomenta[iTrack] += deltaP;
224 
225  // correct for 2PI / PI periodicity
226  auto correctedPhiTheta = detail::ensureThetaBounds(
227  trackMomenta[iTrack][0], trackMomenta[iTrack][1]);
228 
229  trackMomenta[iTrack][0] = correctedPhiTheta.first;
230  trackMomenta[iTrack][1] = correctedPhiTheta.second;
231 
232  // calculate 5x5 covdelta_P matrix
233  // d(d0,z0,phi,theta,qOverP, t)/d(x,y,z,phi,theta,qOverP,
234  // t)-transformation matrix
236  transMat.setZero();
237  transMat(0, 0) = bTrack.DiMat(0, 0);
238  transMat(0, 1) = bTrack.DiMat(0, 1);
239  transMat(1, 0) = bTrack.DiMat(1, 0);
240  transMat(1, 1) = bTrack.DiMat(1, 1);
241  transMat(1, 2) = 1.;
242  transMat(2, 3) = 1.;
243  transMat(3, 4) = 1.;
244  transMat(4, 5) = 1.;
245  transMat(5, 6) = 1.;
246 
247  // some intermediate calculations to get 5x5 matrix
248  // cov(V,V), 4x4 matrix
249  SpacePointSymMatrix VVmat = covDeltaVmat;
250 
251  // cov(V,P)
252  ActsMatrixD<4, 3> VPmat = bTrack.BiMat;
253 
254  // cov(P,P), 3x3 matrix
255  ActsSymMatrixD<3> PPmat;
256  PPmat = bTrack.CiInv +
257  bTrack.BCiMat.transpose() * covDeltaVmat * bTrack.BCiMat;
258 
259  ActsSymMatrixD<7> covMat;
260  covMat.setZero();
261  covMat.block<4, 4>(0, 0) = VVmat;
262  covMat.block<4, 3>(0, 4) = VPmat;
263  covMat.block<3, 4>(4, 0) = VPmat.transpose();
264 
265  covMat.block<3, 3>(4, 4) = PPmat;
266 
267  // covdelta_P calculation
268  covDeltaPmat[iTrack] = transMat * covMat * transMat.transpose();
269  // Calculate chi2 per track.
270  bTrack.chi2 =
271  ((bTrack.deltaQ - bTrack.DiMat * deltaV - bTrack.EiMat * deltaP)
272  .transpose())
273  .dot(bTrack.linTrack.weightAtPCA *
274  (bTrack.deltaQ - bTrack.DiMat * deltaV -
275  bTrack.EiMat * deltaP));
276  newChi2 += bTrack.chi2;
277 
278  ++iTrack;
279  }
280 
281  if (isConstraintFit) {
282  // last term will also be 0 again but only in the first iteration
283  // = calc. vtx in billoir frame - ( isConstraintFit pos. in billoir
284  // frame )
285 
286  SpacePointVector deltaTrk =
287  deltaV -
288  (vertexingOptions.vertexConstraint.fullPosition() - linPoint);
289 
290  newChi2 +=
291  (deltaTrk.transpose())
292  .dot(
293  vertexingOptions.vertexConstraint.fullCovariance().inverse() *
294  deltaTrk);
295  }
296 
297  if (!std::isnormal(newChi2)) {
298  return VertexingError::NumericFailure;
299  }
300 
301  // assign new linearization point (= new vertex position in global frame)
302  linPoint += deltaV;
303  if (newChi2 < chi2) {
304  chi2 = newChi2;
305 
306  SpacePointVector vertexPos(linPoint);
307 
308  fittedVertex.setFullPosition(vertexPos);
309  fittedVertex.setFullCovariance(covDeltaVmat);
310  fittedVertex.setFitQuality(chi2, ndf);
311 
312  std::vector<TrackAtVertex<input_track_t>> tracksAtVertex;
313 
314  std::shared_ptr<PerigeeSurface> perigee =
315  Surface::makeShared<PerigeeSurface>(
316  VectorHelpers::position(vertexPos));
317 
318  iTrack = 0;
319  for (auto& bTrack : billoirTracks) {
320  // new refitted trackparameters
321  BoundVector paramVec;
322  paramVec << 0., 0., trackMomenta[iTrack](0), trackMomenta[iTrack](1),
323  trackMomenta[iTrack](2), 0.;
324 
325  BoundParameters refittedParams(vertexingOptions.geoContext,
326  covDeltaPmat[iTrack], paramVec, perigee);
327 
328  TrackAtVertex<input_track_t> trackVx(bTrack.chi2, refittedParams,
329  bTrack.originalTrack);
330  tracksAtVertex.push_back(std::move(trackVx));
331  ++iTrack;
332  }
333  fittedVertex.setTracksAtVertex(tracksAtVertex);
334  }
335  } // end loop iterations
336  return std::move(fittedVertex);
337 }