ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AdaptiveMultiVertexFitter.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file AdaptiveMultiVertexFitter.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 
12 
13 template <typename input_track_t, typename linearizer_t>
16  State& state, const std::vector<Vertex<input_track_t>*>& verticesToFit,
17  const linearizer_t& linearizer,
18  const VertexingOptions<input_track_t>& vertexingOptions) const {
19  // Set all vertices to fit in the current state
20  state.vertexCollection = verticesToFit;
21 
22  // Perform fit on all vertices simultaneously
23  auto fitRes = fitImpl(state, linearizer, vertexingOptions);
24 
25  if (!fitRes.ok()) {
26  return fitRes.error();
27  }
28 
29  return {};
30 }
31 
32 template <typename input_track_t, typename linearizer_t>
35  State& state, const linearizer_t& linearizer,
36  const VertexingOptions<input_track_t>& vertexingOptions) const {
37  auto& geoContext = vertexingOptions.geoContext;
38 
39  // Reset annealing tool
41 
42  // Indicates how much the vertex positions have shifted
43  // in last fit iteration. Will be false if vertex position
44  // shift was too big. Needed if equilibrium is reached in
45  // annealing procedure but fitter has not fully converged
46  // yet and needs some more iterations until vertex position
47  // shifts between iterations are small (converged).
48  bool isSmallShift = true;
49 
50  // Number of iterations counter
51  unsigned int nIter = 0;
52 
53  // Start iterating
54  while (nIter < m_cfg.maxIterations &&
55  (!state.annealingState.equilibriumReached || !isSmallShift)) {
56  // Initial loop over all vertices in state.vertexCollection
57 
58  for (auto currentVtx : state.vertexCollection) {
59  VertexInfo<input_track_t>& currentVtxInfo = state.vtxInfoMap[currentVtx];
60  currentVtxInfo.relinearize = false;
61  // Store old position of vertex, i.e. seed position
62  // in case of first iteration or position determined
63  // in previous iteration afterwards
64  currentVtxInfo.oldPosition = currentVtx->fullPosition();
65 
66  SpacePointVector dist =
67  currentVtxInfo.oldPosition - currentVtxInfo.linPoint;
68  double perpDist = std::sqrt(dist[0] * dist[0] + dist[1] * dist[1]);
69  // Determine if relinearization is needed
70  if (perpDist > m_cfg.maxDistToLinPoint) {
71  // Relinearization needed, distance too big
72  currentVtxInfo.relinearize = true;
73  // Prepare for fit with new vertex position
74  prepareVertexForFit(state, currentVtx, vertexingOptions);
75  }
76  // Determine if constraint vertex exist
77  if (state.vtxInfoMap[currentVtx].constraintVertex.fullCovariance() !=
78  SpacePointSymMatrix::Zero()) {
79  currentVtx->setFullPosition(
80  state.vtxInfoMap[currentVtx].constraintVertex.fullPosition());
81  currentVtx->setFitQuality(
82  state.vtxInfoMap[currentVtx].constraintVertex.fitQuality());
83  currentVtx->setFullCovariance(
84  state.vtxInfoMap[currentVtx].constraintVertex.fullCovariance());
85  }
86 
87  else if (currentVtx->fullCovariance() == SpacePointSymMatrix::Zero()) {
88  return VertexingError::NoCovariance;
89  }
90  double weight =
91  1. / m_cfg.annealingTool.getWeight(state.annealingState, 1.);
92 
93  currentVtx->setFullCovariance(currentVtx->fullCovariance() * weight);
94 
95  // Set vertexCompatibility for all TrackAtVertex objects
96  // at current vertex
97  setAllVertexCompatibilities(state, currentVtx, vertexingOptions);
98  } // End loop over vertex collection
99 
100  // Now after having estimated all compatibilities of all tracks at
101  // all vertices, run again over all vertices to set track weights
102  // and update the vertex
103  setWeightsAndUpdate(state, linearizer, vertexingOptions);
104  if (!state.annealingState.equilibriumReached) {
105  m_cfg.annealingTool.anneal(state.annealingState);
106  }
107  isSmallShift = checkSmallShift(state);
108 
109  ++nIter;
110  }
111  // Multivertex fit is finished
112 
113  // Check if smoothing is required
114  if (m_cfg.doSmoothing) {
115  doVertexSmoothing(state, geoContext);
116  }
117 
118  return {};
119 }
120 
121 template <typename input_track_t, typename linearizer_t>
124  State& state, Vertex<input_track_t>& newVertex,
125  const linearizer_t& linearizer,
126  const VertexingOptions<input_track_t>& vertexingOptions) const {
127  if (state.vtxInfoMap[&newVertex].trackLinks.empty()) {
128  return VertexingError::EmptyInput;
129  }
130 
131  std::vector<Vertex<input_track_t>*> verticesToFit;
132 
133  // Prepares vtx and tracks for fast estimation method of their
134  // compatibility with vertex
135  auto res = prepareVertexForFit(state, &newVertex, vertexingOptions);
136  if (!res.ok()) {
137  return res.error();
138  }
139  // List of vertices added in last iteration
140  std::vector<Vertex<input_track_t>*> lastIterAddedVertices = {&newVertex};
141  // List of vertices added in current iteration
142  std::vector<Vertex<input_track_t>*> currentIterAddedVertices;
143 
144  // Loop as long as new vertices are found that share tracks with
145  // previously added vertices
146  while (!lastIterAddedVertices.empty()) {
147  for (auto& lastVtxIter : lastIterAddedVertices) {
148  // Loop over all track at current lastVtxIter
149  const std::vector<const input_track_t*>& trks =
150  state.vtxInfoMap[lastVtxIter].trackLinks;
151  for (const auto& trk : trks) {
152  // Retrieve list of links to all vertices that currently use the current
153  // track
154  auto range = state.trackToVerticesMultiMap.equal_range(trk);
155 
156  // Loop over all attached vertices and add those to vertex fit
157  // which are not already in `verticesToFit`
158  for (auto vtxIter = range.first; vtxIter != range.second; ++vtxIter) {
159  auto newVtxIter = vtxIter->second;
160  if (!isAlreadyInList(newVtxIter, verticesToFit)) {
161  // Add newVtxIter to verticesToFit
162  verticesToFit.push_back(newVtxIter);
163 
164  // Add newVtxIter vertex to currentIterAddedVertices
165  // if vertex != lastVtxIter
166  if (newVtxIter != lastVtxIter) {
167  currentIterAddedVertices.push_back(newVtxIter);
168  }
169  }
170  } // End for loop over linksToVertices
171  }
172  } // End loop over lastIterAddedVertices
173 
174  lastIterAddedVertices = currentIterAddedVertices;
175  currentIterAddedVertices.clear();
176  } // End while loop
177 
178  state.vertexCollection = verticesToFit;
179 
180  // Perform fit on all added vertices
181  auto fitRes = fitImpl(state, linearizer, vertexingOptions);
182  if (!fitRes.ok()) {
183  return fitRes.error();
184  }
185 
186  return {};
187 }
188 
189 template <typename input_track_t, typename linearizer_t>
193  const std::vector<Vertex<input_track_t>*>& verticesVec) const {
194  return std::find(verticesVec.begin(), verticesVec.end(), vtx) !=
195  verticesVec.end();
196 }
197 
198 template <typename input_track_t, typename linearizer_t>
201  State& state, Vertex<input_track_t>* vtx,
202  const VertexingOptions<input_track_t>& vertexingOptions) const {
203  // The current vertex info object
204  auto& currentVtxInfo = state.vtxInfoMap[vtx];
205  // The seed position
206  const Vector3D& seedPos = currentVtxInfo.seedPosition.template head<3>();
207 
208  // Loop over all tracks at current vertex
209  for (const auto& trk : currentVtxInfo.trackLinks) {
210  auto res = m_cfg.ipEst.estimate3DImpactParameters(
211  vertexingOptions.geoContext, vertexingOptions.magFieldContext,
212  m_extractParameters(*trk), seedPos);
213  if (!res.ok()) {
214  return res.error();
215  }
216  // Set ip3dParams for current trackAtVertex
217  currentVtxInfo.ip3dParams.emplace(trk, *(res.value()));
218  }
219  return {};
220 }
221 
222 template <typename input_track_t, typename linearizer_t>
226  State& state, Vertex<input_track_t>* currentVtx,
227  const VertexingOptions<input_track_t>& vertexingOptions) const {
228  VertexInfo<input_track_t>& currentVtxInfo = state.vtxInfoMap[currentVtx];
229 
230  // Loop over tracks at current vertex and
231  // estimate compatibility with vertex
232  for (const auto& trk : currentVtxInfo.trackLinks) {
233  auto& trkAtVtx =
234  state.tracksAtVerticesMap.at(std::make_pair(trk, currentVtx));
235  // Recover from cases where linearization point != 0 but
236  // more tracks were added later on
237  if (currentVtxInfo.ip3dParams.find(trk) ==
238  currentVtxInfo.ip3dParams.end()) {
239  auto res = m_cfg.ipEst.estimate3DImpactParameters(
240  vertexingOptions.geoContext, vertexingOptions.magFieldContext,
241  m_extractParameters(*trk),
242  VectorHelpers::position(currentVtxInfo.linPoint));
243  if (!res.ok()) {
244  return res.error();
245  }
246  // Set ip3dParams for current trackAtVertex
247  currentVtxInfo.ip3dParams.emplace(trk, *(res.value()));
248  }
249  // Set compatibility with current vertex
250  auto compRes = m_cfg.ipEst.get3dVertexCompatibility(
251  vertexingOptions.geoContext, &(currentVtxInfo.ip3dParams.at(trk)),
252  VectorHelpers::position(currentVtxInfo.oldPosition));
253  if (!compRes.ok()) {
254  return compRes.error();
255  }
256  trkAtVtx.vertexCompatibility = *compRes;
257  }
258  return {};
259 }
260 
261 template <typename input_track_t, typename linearizer_t>
264  State& state, const linearizer_t& linearizer,
265  const VertexingOptions<input_track_t>& vertexingOptions) const {
266  for (auto vtx : state.vertexCollection) {
267  VertexInfo<input_track_t>& currentVtxInfo = state.vtxInfoMap[vtx];
268 
269  for (const auto& trk : currentVtxInfo.trackLinks) {
270  auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx));
271 
272  // Set trackWeight for current track
273  double currentTrkWeight = m_cfg.annealingTool.getWeight(
274  state.annealingState, trkAtVtx.vertexCompatibility,
275  collectTrackToVertexCompatibilities(state, trk));
276  trkAtVtx.trackWeight = currentTrkWeight;
277 
278  if (trkAtVtx.trackWeight > m_cfg.minWeight) {
279  // Check if linearization state exists or need to be relinearized
280  if (trkAtVtx.linearizedState.covarianceAtPCA ==
281  BoundSymMatrix::Zero() ||
282  state.vtxInfoMap[vtx].relinearize) {
283  auto result = linearizer.linearizeTrack(
284  m_extractParameters(*trk), state.vtxInfoMap[vtx].oldPosition,
285  vertexingOptions.geoContext, vertexingOptions.magFieldContext);
286  if (!result.ok()) {
287  return result.error();
288  }
289  trkAtVtx.linearizedState = *result;
290  state.vtxInfoMap[vtx].linPoint = state.vtxInfoMap[vtx].oldPosition;
291  }
292  // Update the vertex with the new track
293  KalmanVertexUpdater::updateVertexWithTrack<input_track_t>(*vtx,
294  trkAtVtx);
295  } else {
296  ACTS_VERBOSE("Track weight too low. Skip track.");
297  }
298  } // End loop over tracks at vertex
299 
300  ACTS_VERBOSE("New vertex position: " << vtx->fullPosition());
301  } // End loop over vertex collection
302 
303  return {};
304 }
305 
306 template <typename input_track_t, typename linearizer_t>
307 std::vector<double>
310  const input_track_t* trk) const {
311  std::vector<double> trkToVtxCompatibilities;
312  trkToVtxCompatibilities.reserve(state.vertexCollection.size());
313  auto range = state.trackToVerticesMultiMap.equal_range(trk);
314 
315  for (auto vtxIter = range.first; vtxIter != range.second; ++vtxIter) {
316  trkToVtxCompatibilities.push_back(
317  state.tracksAtVerticesMap.at(std::make_pair(trk, vtxIter->second))
318  .vertexCompatibility);
319  }
320 
321  return trkToVtxCompatibilities;
322 }
323 
324 template <typename input_track_t, typename linearizer_t>
326  input_track_t, linearizer_t>::checkSmallShift(State& state) const {
327  for (auto vtx : state.vertexCollection) {
328  auto diff = state.vtxInfoMap[vtx].oldPosition.template head<3>() -
329  vtx->fullPosition().template head<3>();
330  const auto& vtxWgt =
331  (vtx->fullCovariance().template block<3, 3>(0, 0)).inverse();
332  double relativeShift = diff.dot(vtxWgt * diff);
333  if (relativeShift > m_cfg.maxRelativeShift) {
334  return false;
335  }
336  }
337  return true;
338 }
339 
340 template <typename input_track_t, typename linearizer_t>
343  for (const auto vtx : state.vertexCollection) {
344  for (const auto trk : state.vtxInfoMap[vtx].trackLinks) {
345  KalmanVertexTrackUpdater::update<input_track_t>(
346  geoContext, state.tracksAtVerticesMap.at(std::make_pair(trk, vtx)),
347  *vtx);
348  }
349  }
350 }