ATLAS Offline Software
DistributedKalmanFilter.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3  */
4 
6 // DistributedKalmanFilter.cxx
7 // Source file for class DistributedKalmanFilter
9 // (c) ATLAS Detector software
11 // Author: Dmitry Emeliyanov, RAL
12 // D.Emeliyanov@rl.ac.uk
14 
15 #include "GaudiKernel/SystemOfUnits.h"
16 
17 #include "TrkTrack/Track.h"
18 #include "TrkTrack/TrackStateOnSurface.h"
19 #include "TrkParameters/TrackParameters.h"
20 
21 #include "TrkRIO_OnTrack/RIO_OnTrack.h"
22 #include "TrkPrepRawData/PrepRawData.h"
23 #include "TrkSurfaces/Surface.h"
24 #include "TrkSurfaces/RectangleBounds.h"
25 #include "TrkSurfaces/DiscBounds.h"
26 #include "TrkEventPrimitives/PropDirection.h"
27 
28 #include "TrkEventPrimitives/ParticleHypothesis.h"
29 
30 #include "TrkDetElementBase/TrkDetElementBase.h"
31 #include "TrkMeasurementBase/MeasurementBase.h"
32 
33 #include "TrkEventPrimitives/FitQuality.h"
34 #include "TrkEventPrimitives/FitQualityOnSurface.h"
35 
36 #include "TrkEventUtils/PrepRawDataComparisonFunction.h"
37 
38 #include "AtlasDetDescr/AtlasDetectorID.h"
39 
40 //DKF stuff
41 
42 #include "TrkDistributedKalmanFilter/DistributedKalmanFilter.h"
43 #include "TrkDistributedKalmanFilter/TrkBaseNode.h"
44 #include "TrkDistributedKalmanFilter/TrkFilteringNodes.h"
45 #include "TrkDistributedKalmanFilter/TrkPlanarSurface.h"
46 #include "TrkDistributedKalmanFilter/TrkTrackState.h"
47 
48 #include "TrkTrack/TrackInfo.h"
49 
50 #include <memory>
51 
52 #include <vector>
53 #include <algorithm>
54 #include <cmath>
55 
56 namespace {
57  std::unique_ptr<Trk::Perigee>
58  createMeasuredPerigee(Trk::TrkTrackState* pTS) {
59  AmgSymMatrix(5) pM;
60  for (int i = 0; i < 5; i++) {
61  for (int j = 0; j < 5; j++) {
62  (pM) (i, j) = pTS->getTrackCovariance(i, j);
63  }
64  }
65  const Trk::PerigeeSurface perSurf {};
66  return std::make_unique<Trk::Perigee>(pTS->getTrackState(0),
67  pTS->getTrackState(1),
68  pTS->getTrackState(2),
69  pTS->getTrackState(3),
70  pTS->getTrackState(4),
71  perSurf,
72  std::move(pM));
73  }
74 
75  void
76  matrixInversion5x5(double a[5][5]) {
77  /**** 5x5 matrix inversion by Gaussian elimination ****/
78  int i, j, k, l;
79  double factor;
80  double temp[5];
81  double b[5][5];
82 
83  // Set b to I
84 
85  memset(&b[0][0], 0, sizeof(b));
86  b[0][0] = 1.0;
87  b[1][1] = 1.0;
88  b[2][2] = 1.0;
89  b[3][3] = 1.0;
90  b[4][4] = 1.0;
91 
92  for (i = 0; i < 5; i++) {
93  for (j = i + 1; j < 5; j++)
94  if (std::abs(a[i][i]) < std::abs(a[j][i])) {
95  for (l = 0; l < 5; l++) temp[l] = a[i][l];
96  for (l = 0; l < 5; l++) a[i][l] = a[j][l];
97  for (l = 0; l < 5; l++) a[j][l] = temp[l];
98  for (l = 0; l < 5; l++) temp[l] = b[i][l];
99  for (l = 0; l < 5; l++) b[i][l] = b[j][l];
100  for (l = 0; l < 5; l++) b[j][l] = temp[l];
101  }
102  factor = a[i][i];
103  for (j = 4; j > -1; j--) {
104  b[i][j] /= factor;
105  a[i][j] /= factor;
106  }
107  for (j = i + 1; j < 5; j++) {
108  factor = -a[j][i];
109  for (k = 0; k < 5; k++) {
110  a[j][k] += a[i][k] * factor;
111  b[j][k] += b[i][k] * factor;
112  }
113  }
114  }
115  for (i = 4; i > 0; i--) {
116  for (j = i - 1; j > -1; j--) {
117  factor = -a[j][i];
118  for (k = 0; k < 5; k++) {
119  a[j][k] += a[i][k] * factor;
120  b[j][k] += b[i][k] * factor;
121  }
122  }
123  }
124  for (i = 0; i < 5; i++) for (j = 0; j < 5; j++) a[i][j] = b[i][j];
125  }
126 }
127 
128 
129 // constructor
130 
131 Trk::DistributedKalmanFilter::DistributedKalmanFilter(const std::string& t,
132  const std::string& n,
133  const IInterface* p)
134  : AthAlgTool(t, n, p)
135  , m_idHelper(nullptr) {
136  // AlgTool stuff
137  declareInterface<ITrackFitter>(this);
138  // ME temporary fix
139  declareProperty("sortingReferencePoint", m_option_sortingRefPoint);
140 }
141 
142 // destructor
143 Trk::DistributedKalmanFilter::~DistributedKalmanFilter() = default;
144 
145 // initialize
146 StatusCode Trk::DistributedKalmanFilter::initialize() {
147  //StatusCode s = AthAlgTool::finalize();
148 
149  ATH_CHECK(detStore()->retrieve(m_idHelper, "AtlasID"));
150 
151  ATH_CHECK(m_fieldCacheCondObjInputKey.initialize());
152  ATH_CHECK(m_ROTcreator.retrieve());
153  ATH_CHECK(m_extrapolator.retrieve());
154 
155  ATH_MSG_VERBOSE("initialize() successful in Trk::DistributedKalmanFilter");
156  return StatusCode::SUCCESS;
157 }
158 
159 // finalize
160 StatusCode Trk::DistributedKalmanFilter::finalize() {
161  // init message stream
162  //StatusCode sc = AthAlgTool::finalize();
163  ATH_MSG_VERBOSE(" finalize() successful in Trk::DistributedKalmanFilter");
164  return StatusCode::SUCCESS;
165 }
166 
167 // helper operator for STL min_element
168 struct minTrkPar {
169  bool operator () (const Trk::TrackParameters* one, const Trk::TrackParameters* two) const {
170  return(one->position().mag() < two->position().mag());
171  };
172 };
173 
174 // refit a track
175 std::unique_ptr<Trk::Track>
176 Trk::DistributedKalmanFilter::fit(const EventContext& ctx,
177  const Trk::Track& inputTrack,
178  const RunOutlierRemoval runOutlier,
179  const ParticleHypothesis matEffects) const {
180  ATH_MSG_VERBOSE("--> enter DistributedKalmanFilter::fit(Track,,)");
181 
182  // protection against track not having any parameters
183  if (inputTrack.trackParameters()->empty()) {
184  ATH_MSG_ERROR("need estimated track parameters near " << "origin, reject fit");
185  return nullptr;
186  }
187  // protection against not having RIO_OnTracks
188  if (inputTrack.measurementsOnTrack()->empty()) {
189  ATH_MSG_ERROR("try to refit track with empty " << "vec<RIO_OnTrack>, reject fit");
190  return nullptr;
191  }
192 
193  // create PrepRawData subset
194  PrepRawDataSet prepRDColl;
195 
196  // collect PrepRawData pointers from input track ROTs
197  DataVector<const MeasurementBase>::const_iterator it = inputTrack.measurementsOnTrack()->begin();
198  DataVector<const MeasurementBase>::const_iterator itEnd = inputTrack.measurementsOnTrack()->end();
199  for (; it != itEnd; ++it) {
200  if (!(*it)) {
201  ATH_MSG_WARNING("This track contains empty MeasurementBase "
202  << "objects! Skipped this MB..");
203  } else {
204  const RIO_OnTrack* testROT = dynamic_cast<const RIO_OnTrack*>(*it);
205  const PrepRawData* prepRD = (testROT) ? (testROT)->prepRawData() : nullptr;
206  if (!prepRD) {
207  ATH_MSG_WARNING("RIO_OnTrack->prepRawRata() "
208  << "returns no object, ignore... ");
209  } else {
210  prepRDColl.push_back(prepRD);
211  }
212  }
213  }
214  const TrackParameters* minPar =
215  *(std::min_element(inputTrack.trackParameters()->begin(),
216  inputTrack.trackParameters()->end(),
217  minTrkPar()));
218 
219  // fit set of PrepRawData using main method, start with first Track Parameter in inputTrack
220  ATH_MSG_VERBOSE("call fit(PRDSet,TP,,)");
221  return fit(ctx, prepRDColl, *minPar, runOutlier, matEffects);
222 }
223 
224 // @TODO remove ? unused :
225 double Trk::DistributedKalmanFilter::integrate(double Rk[5],
226  TrkPlanarSurface* pSB,
227  TrkPlanarSurface* pSE, double* Rf,
228  MagField::AtlasFieldCache& fieldCache) const {
229  const double C = 0.02999975;
230  const int nStepMax = 5;
231 
232  double sint, cost, sinf, cosf;
233  double gP[3], lP[3], gV[3], a, b, c, descr, CQ, Ac, Av;
234  double V[3], P[3], D[4], gB[3];
235  int i, nStep;
236  double sl, ds, path = 0.0;
237 
238  sint = sin(Rk[3]);
239  cosf = cos(Rk[2]);
240  sinf = sin(Rk[2]);
241  cost = cos(Rk[3]);
242  gV[0] = sint * cosf;
243  gV[1] = sint * sinf;
244  gV[2] = cost;
245  CQ = C * Rk[4];
246 
247  if (pSB != nullptr) {
248  lP[0] = Rk[0];
249  lP[1] = Rk[1];
250  lP[2] = 0.0;
251  pSB->transformPointToGlobal(lP, gP);
252  } else {
253  gP[0] = -Rk[0] * sinf;
254  gP[1] = Rk[0] * cosf;
255  gP[2] = Rk[1];
256  }
257  for (i = 0; i < 4; i++) D[i] = pSE->getPar(i);
258 
259  getMagneticField(gP, gB, fieldCache);
260 
261  c = D[0] * gP[0] + D[1] * gP[1] + D[2] * gP[2] + D[3];
262  b = D[0] * gV[0] + D[1] * gV[1] + D[2] * gV[2];
263  a = 0.5 * CQ * (gB[0] * (D[1] * gV[2] - D[2] * gV[1]) +
264  gB[1] * (D[2] * gV[0] - D[0] * gV[2]) +
265  gB[2] * (D[0] * gV[1] - D[1] * gV[0]));
266 
267  descr = b * b - 4.0 * a * c;
268 
269  if (descr < 0.0) {
270  printf("D<0 - extrapolation failed\n");
271  return 0;
272  }
273 
274  nStep = nStepMax;
275  path = 0.0;
276  while (nStep > 0) {
277  c = D[0] * gP[0] + D[1] * gP[1] + D[2] * gP[2] + D[3];
278  b = D[0] * gV[0] + D[1] * gV[1] + D[2] * gV[2];
279  a = 0.5 * CQ *
280  (gB[0] * (D[1] * gV[2] - D[2] * gV[1]) +
281  gB[1] * (D[2] * gV[0] - D[0] * gV[2]) +
282  gB[2] * (D[0] * gV[1] - D[1] * gV[0]));
283  sl = -c / b;
284  sl = sl * (1 - a * sl / b);
285  ds = sl / nStep;
286  path += ds;
287  Av = ds * CQ;
288  Ac = 0.5 * ds * Av;
289 
290  P[0] = gP[0] + gV[0] * ds + Ac * (gV[1] * gB[2] - gV[2] * gB[1]);
291  P[1] = gP[1] + gV[1] * ds + Ac * (gV[2] * gB[0] - gV[0] * gB[2]);
292  P[2] = gP[2] + gV[2] * ds + Ac * (gV[0] * gB[1] - gV[1] * gB[0]);
293  V[0] = gV[0] + Av * (gV[1] * gB[2] - gV[2] * gB[1]);
294  V[1] = gV[1] + Av * (gV[2] * gB[0] - gV[0] * gB[2]);
295  V[2] = gV[2] + Av * (gV[0] * gB[1] - gV[1] * gB[0]);
296  for (i = 0; i < 3; i++) {
297  gV[i] = V[i];
298  gP[i] = P[i];
299  }
300  getMagneticField(gP, gB, fieldCache);
301  nStep--;
302  }
303  pSE->transformPointToLocal(P, lP);
304 
305  Rf[0] = lP[0];
306  Rf[1] = lP[1];
307  Rf[2] = atan2(V[1], V[0]);
308  Rf[3] = acos(V[2]);
309  Rf[4] = Rk[4];
310 
311  return path;
312 }
313 
314 std::unique_ptr<Trk::TrkTrackState>
315 Trk::DistributedKalmanFilter::extrapolate(Trk::TrkTrackState* pTS,
316  Trk::TrkPlanarSurface* pSB,
317  Trk::TrkPlanarSurface* pSE,
318  MagField::AtlasFieldCache& fieldCache) const {
319  const double C = 0.02999975;
320  const double minStep = 30.0;
321 
322  double J[5][5], Rf[5], AG[5][5], Gi[5][5], Gf[5][5], A[5][5];
323  int i, j, m;
324 
325  bool samePlane = false;
326 
327  if (pSB != nullptr) {
328  double diff = 0.0;
329  for (i = 0; i < 4; i++) diff += fabs(pSE->getPar(i) - pSB->getPar(i));
330  if (diff < 1e-5) {
331  samePlane = true;
332  }
333  }
334 
335  if (!samePlane) {
336  double gP[3], gPi[3], lP[3], gV[3], a, b, c, s, J0[7][5], descr, CQ, Ac, Av, Cc;
337  double V[3], P[3], M[3][3], D[4], Jm[7][7],
338  J1[5][7], gB[3], gBi[3], gBf[3], dBds[3], Buf[5][7], DVx, DVy, DVz;
339  int nStep, nStepMax;
340  double sl, ds = 0.0;
341 
342  //numericalJacobian(pTS,pSB,pSE,J, fieldCache);
343  double sint, cost, sinf, cosf;
344  sint = sin(pTS->getTrackState(3));
345  cosf = cos(pTS->getTrackState(2));
346  sinf = sin(pTS->getTrackState(2));
347  cost = cos(pTS->getTrackState(3));
348  gV[0] = sint * cosf;
349  gV[1] = sint * sinf;
350  gV[2] = cost;
351  CQ = C * pTS->getTrackState(4);
352 
353  memset(&J0[0][0], 0, sizeof(J0));
354 
355  if (pSB != nullptr) {
356  double L[3][3];
357  lP[0] = pTS->getTrackState(0);
358  lP[1] = pTS->getTrackState(1);
359  lP[2] = 0.0;
360  pSB->transformPointToGlobal(lP, gP);
361  for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) L[i][j] = pSB->getInvRotMatrix(i, j);
362 
363  J0[0][0] = L[0][0];
364  J0[0][1] = L[0][1];
365  J0[1][0] = L[1][0];
366  J0[1][1] = L[1][1];
367  J0[2][0] = L[2][0];
368  J0[2][1] = L[2][1];
369  J0[3][2] = -sinf * sint;
370  J0[3][3] = cosf * cost;
371  J0[4][2] = cosf * sint;
372  J0[4][3] = sinf * cost;
373  J0[5][3] = -sint;
374  J0[6][4] = 1.0;
375  } else {
376  gP[0] = -pTS->getTrackState(0) * sinf;
377  gP[1] = pTS->getTrackState(0) * cosf;
378  gP[2] = pTS->getTrackState(1);
379  J0[0][0] = -sinf;
380  J0[0][2] = -pTS->getTrackState(0) * cosf;
381  J0[1][0] = cosf;
382  J0[1][2] = -pTS->getTrackState(0) * sinf;
383  J0[2][1] = 1.0;
384  J0[3][2] = -sinf * sint;
385  J0[3][3] = cosf * cost;
386  J0[4][2] = cosf * sint;
387  J0[4][3] = sinf * cost;
388  J0[5][3] = -sint;
389  J0[6][4] = 1.0;
390  }
391  for (i = 0; i < 4; i++) D[i] = pSE->getPar(i);
392  for (i = 0; i < 3; i++) gPi[i] = gP[i];
393 
394  getMagneticField(gP, gB, fieldCache);
395 
396  for (i = 0; i < 3; i++) gBi[i] = gB[i];
397 
398  c = D[0] * gP[0] + D[1] * gP[1] + D[2] * gP[2] + D[3];
399  b = D[0] * gV[0] + D[1] * gV[1] + D[2] * gV[2];
400  a = 0.5 * CQ * (gB[0] * (D[1] * gV[2] - D[2] * gV[1]) +
401  gB[1] * (D[2] * gV[0] - D[0] * gV[2]) +
402  gB[2] * (D[0] * gV[1] - D[1] * gV[0]));
403 
404  descr = b * b - 4.0 * a * c;
405 
406  if (descr < 0.0) {
407  // printf("D<0 - extrapolation failed\n");
408  return nullptr;
409  }
410 
411  bool useExpansion = true;
412  double ratio = 4 * a * c / (b * b);
413 
414  if (fabs(ratio) > 0.1) useExpansion = false;
415 
416  if (useExpansion) {
417  sl = -c / b;
418  sl = sl * (1 - a * sl / b);
419  } else {
420  int signb = (b < 0.0) ? -1 : 1;
421  sl = (-b + signb * sqrt(descr)) / (2 * a);
422  }
423 
424  if (fabs(sl) < minStep) nStepMax = 1;
425  else {
426  nStepMax = (int) (fabs(sl) / minStep) + 1;
427  }
428  if ((nStepMax < 0) || (nStepMax > 1000)) {
429  return nullptr;
430  }
431  Av = sl * CQ;
432  Ac = 0.5 * sl * Av;
433  DVx = gV[1] * gB[2] - gV[2] * gB[1];
434  DVy = gV[2] * gB[0] - gV[0] * gB[2];
435  DVz = gV[0] * gB[1] - gV[1] * gB[0];
436 
437  P[0] = gP[0] + gV[0] * sl + Ac * DVx;
438  P[1] = gP[1] + gV[1] * sl + Ac * DVy;
439  P[2] = gP[2] + gV[2] * sl + Ac * DVz;
440  V[0] = gV[0] + Av * DVx;
441  V[1] = gV[1] + Av * DVy;
442  V[2] = gV[2] + Av * DVz;
443 
444  getMagneticField(P, gB, fieldCache);
445 
446  for (i = 0; i < 3; i++) gBf[i] = gB[i];
447  for (i = 0; i < 3; i++) {
448  dBds[i] = (gBf[i] - gBi[i]) / sl;
449  gB[i] = gBi[i];
450  }
451  nStep = nStepMax;
452  while (nStep > 0) {
453  c = D[0] * gP[0] + D[1] * gP[1] + D[2] * gP[2] + D[3];
454  b = D[0] * gV[0] + D[1] * gV[1] + D[2] * gV[2];
455  a = 0.5 * CQ * (gB[0] * (D[1] * gV[2] - D[2] * gV[1]) +
456  gB[1] * (D[2] * gV[0] - D[0] * gV[2]) +
457  gB[2] * (D[0] * gV[1] - D[1] * gV[0]));
458 
459  ratio = 4 * a * c / (b * b);
460  useExpansion = fabs(ratio) <= 0.1;
461 
462  if (useExpansion) {
463  sl = -c / b;
464  sl = sl * (1 - a * sl / b);
465  } else {
466  descr = b * b - 4.0 * a * c;
467  if (descr < 0.0) {
468  // printf("D<0 - extrapolation failed\n");
469  return nullptr;
470  }
471  int signb = (b < 0.0) ? -1 : 1;
472  sl = (-b + signb * sqrt(descr)) / (2 * a);
473  }
474 
475  ds = sl / nStep;
476  Av = ds * CQ;
477  Ac = 0.5 * ds * Av;
478  DVx = gV[1] * gB[2] - gV[2] * gB[1];
479  DVy = gV[2] * gB[0] - gV[0] * gB[2];
480  DVz = gV[0] * gB[1] - gV[1] * gB[0];
481 
482  P[0] = gP[0] + gV[0] * ds + Ac * DVx;
483  P[1] = gP[1] + gV[1] * ds + Ac * DVy;
484  P[2] = gP[2] + gV[2] * ds + Ac * DVz;
485  V[0] = gV[0] + Av * DVx;
486  V[1] = gV[1] + Av * DVy;
487  V[2] = gV[2] + Av * DVz;
488  for (i = 0; i < 3; i++) {
489  gV[i] = V[i];
490  gP[i] = P[i];
491  }
492  for (i = 0; i < 3; i++) gB[i] += dBds[i] * ds;
493  nStep--;
494  }
495  pSE->transformPointToLocal(gP, lP);
496  Rf[0] = lP[0];
497  Rf[1] = lP[1];
498  Rf[2] = atan2(V[1], V[0]);
499 
500  if (fabs(V[2]) > 1.0) {
501  return nullptr;
502  }
503 
504  Rf[3] = acos(V[2]);
505  Rf[4] = pTS->getTrackState(4);
506 
507  gV[0] = sint * cosf;
508  gV[1] = sint * sinf;
509  gV[2] = cost;
510 
511  for (i = 0; i < 4; i++) D[i] = pSE->getPar(i);
512  for (i = 0; i < 3; i++) gP[i] = gPi[i];
513 
514  for (i = 0; i < 3; i++) {
515  gB[i] = 0.5 * (gBi[i] + gBf[i]);
516  }
517 
518  c = D[0] * gP[0] + D[1] * gP[1] + D[2] * gP[2] + D[3];
519  b = D[0] * gV[0] + D[1] * gV[1] + D[2] * gV[2];
520  a = CQ * (gB[0] * (D[1] * gV[2] - D[2] * gV[1]) +
521  gB[1] * (D[2] * gV[0] - D[0] * gV[2]) +
522  gB[2] * (D[0] * gV[1] - D[1] * gV[0]));
523 
524  ratio = 4 * a * c / (b * b);
525  useExpansion = fabs(ratio) <= 0.1;
526 
527  if (useExpansion) {
528  s = -c / b;
529  s = s * (1 - a * s / b);
530  } else {
531  descr = b * b - 4.0 * a * c;
532  if (descr < 0.0) {
533  // printf("D<0 - extrapolation failed\n");
534  return nullptr;
535  }
536  int signb = (b < 0.0) ? -1 : 1;
537  s = (-b + signb * sqrt(descr)) / (2 * a);
538  }
539 
540  Av = s * CQ;
541  Ac = 0.5 * s * Av;
542  Cc = 0.5 * s * s * C;
543 
544  DVx = gV[1] * gB[2] - gV[2] * gB[1];
545  DVy = gV[2] * gB[0] - gV[0] * gB[2];
546  DVz = gV[0] * gB[1] - gV[1] * gB[0];
547 
548  P[0] = gP[0] + gV[0] * s + Ac * DVx;
549  P[1] = gP[1] + gV[1] * s + Ac * DVy;
550  P[2] = gP[2] + gV[2] * s + Ac * DVz;
551 
552  V[0] = gV[0] + Av * DVx;
553  V[1] = gV[1] + Av * DVy;
554  V[2] = gV[2] + Av * DVz;
555 
556  pSE->transformPointToLocal(P, lP);
557 
558  memset(&Jm[0][0], 0, sizeof(Jm));
559 
560  for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) M[i][j] = pSE->getRotMatrix(i, j);
561 
562  double coeff[3], dadVx, dadVy, dadVz, dadQ, dsdx, dsdy, dsdz, dsdVx, dsdVy, dsdVz, dsdQ;
563  coeff[0] = -c * c / (b * b * b);
564  coeff[1] = c * (1.0 + 3.0 * c * a / (b * b)) / (b * b);
565  coeff[2] = -(1.0 + 2.0 * c * a / (b * b)) / b;
566 
567  dadVx = 0.5 * CQ * (-D[1] * gB[2] + D[2] * gB[1]);
568  dadVy = 0.5 * CQ * (D[0] * gB[2] - D[2] * gB[0]);
569  dadVz = 0.5 * CQ * (-D[0] * gB[1] + D[1] * gB[0]);
570  dadQ = 0.5 * C * (D[0] * DVx + D[1] * DVy + D[2] * DVz);
571 
572  dsdx = coeff[2] * D[0];
573  dsdy = coeff[2] * D[1];
574  dsdz = coeff[2] * D[2];
575  dsdVx = coeff[0] * dadVx + coeff[1] * D[0];
576  dsdVy = coeff[0] * dadVy + coeff[1] * D[1];
577  dsdVz = coeff[0] * dadVz + coeff[1] * D[2];
578  dsdQ = coeff[0] * dadQ;
579 
580  Jm[0][0] = 1.0 + V[0] * dsdx;
581  Jm[0][1] = V[0] * dsdy;
582  Jm[0][2] = V[0] * dsdz;
583 
584  Jm[0][3] = s + V[0] * dsdVx;
585  Jm[0][4] = V[0] * dsdVy + Ac * gB[2];
586  Jm[0][5] = V[0] * dsdVz - Ac * gB[1];
587  Jm[0][6] = V[0] * dsdQ + Cc * DVx;
588 
589  Jm[1][0] = V[1] * dsdx;
590  Jm[1][1] = 1.0 + V[1] * dsdy;
591  Jm[1][2] = V[1] * dsdz;
592 
593  Jm[1][3] = V[1] * dsdVx - Ac * gB[2];
594  Jm[1][4] = s + V[1] * dsdVy;
595  Jm[1][5] = V[1] * dsdVz + Ac * gB[0];
596  Jm[1][6] = V[1] * dsdQ + Cc * DVy;
597 
598  Jm[2][0] = V[2] * dsdx;
599  Jm[2][1] = V[2] * dsdy;
600  Jm[2][2] = 1.0 + V[2] * dsdz;
601  Jm[2][3] = V[2] * dsdVx + Ac * gB[1];
602  Jm[2][4] = V[2] * dsdVy - Ac * gB[0];
603  Jm[2][5] = s + V[2] * dsdVz;
604  Jm[2][6] = V[2] * dsdQ + Cc * DVz;
605 
606  Jm[3][0] = dsdx * CQ * DVx;
607  Jm[3][1] = dsdy * CQ * DVx;
608  Jm[3][2] = dsdz * CQ * DVx;
609 
610  Jm[3][3] = 1.0 + dsdVx * CQ * DVx;
611  Jm[3][4] = CQ * (dsdVy * DVx + s * gB[2]);
612  Jm[3][5] = CQ * (dsdVz * DVx - s * gB[1]);
613 
614  Jm[3][6] = (CQ * dsdQ + C * s) * DVx;
615 
616  Jm[4][0] = dsdx * CQ * DVy;
617  Jm[4][1] = dsdy * CQ * DVy;
618  Jm[4][2] = dsdz * CQ * DVy;
619 
620  Jm[4][3] = CQ * (dsdVx * DVy - s * gB[2]);
621  Jm[4][4] = 1.0 + dsdVy * CQ * DVy;
622  Jm[4][5] = CQ * (dsdVz * DVy + s * gB[0]);
623 
624  Jm[4][6] = (CQ * dsdQ + C * s) * DVy;
625 
626  Jm[5][0] = dsdx * CQ * DVz;
627  Jm[5][1] = dsdy * CQ * DVz;
628  Jm[5][2] = dsdz * CQ * DVz;
629  Jm[5][3] = CQ * (dsdVx * DVz + s * gB[1]);
630  Jm[5][4] = CQ * (dsdVy * DVz - s * gB[0]);
631  Jm[5][5] = 1.0 + dsdVz * CQ * DVz;
632  Jm[5][6] = (CQ * dsdQ + C * s) * DVz;
633 
634  Jm[6][6] = 1.0;
635 
636  memset(&J1[0][0], 0, sizeof(J1));
637 
638  J1[0][0] = M[0][0];
639  J1[0][1] = M[0][1];
640  J1[0][2] = M[0][2];
641  J1[1][0] = M[1][0];
642  J1[1][1] = M[1][1];
643  J1[1][2] = M[1][2];
644  J1[2][3] = -V[1] / (V[0] * V[0] + V[1] * V[1]);
645  J1[2][4] = V[0] / (V[0] * V[0] + V[1] * V[1]);
646  J1[3][5] = -1.0 / sqrt(1 - V[2] * V[2]);
647  J1[4][6] = 1.0;
648 
649  for (i = 0; i < 7; i++) {
650  for (j = 0; j < 2; j++)
651  Buf[j][i] = J1[j][0] * Jm[0][i] + J1[j][1] * Jm[1][i] + J1[j][2] * Jm[2][i];
652  Buf[2][i] = J1[2][3] * Jm[3][i] + J1[2][4] * Jm[4][i];
653  Buf[3][i] = J1[3][5] * Jm[5][i];
654  Buf[4][i] = Jm[6][i];
655  }
656 
657  if (pSB != nullptr) {
658  for (i = 0; i < 5; i++) {
659  J[i][0] = Buf[i][0] * J0[0][0] + Buf[i][1] * J0[1][0] + Buf[i][2] * J0[2][0];
660  J[i][1] = Buf[i][0] * J0[0][1] + Buf[i][1] * J0[1][1] + Buf[i][2] * J0[2][1];
661  J[i][2] = Buf[i][3] * J0[3][2] + Buf[i][4] * J0[4][2];
662  J[i][3] = Buf[i][3] * J0[3][3] + Buf[i][4] * J0[4][3] + Buf[i][5] * J0[5][3];
663  J[i][4] = Buf[i][6];
664  }
665  } else {
666  for (i = 0; i < 5; i++) {
667  J[i][0] = Buf[i][0] * J0[0][0] + Buf[i][1] * J0[1][0];
668  J[i][1] = Buf[i][2];
669  J[i][2] = Buf[i][0] * J0[0][2] + Buf[i][1] * J0[1][2] + Buf[i][3] * J0[3][2] + Buf[i][4] * J0[4][2];
670  J[i][3] = Buf[i][3] * J0[3][3] + Buf[i][4] * J0[4][3] + Buf[i][5] * J0[5][3];
671  J[i][4] = Buf[i][6];
672  }
673  }
674  } else {
675  Rf[0] = pTS->getTrackState(0);
676  Rf[1] = pTS->getTrackState(1);
677  Rf[2] = pTS->getTrackState(2);
678  Rf[3] = pTS->getTrackState(3);
679  Rf[4] = pTS->getTrackState(4);
680  memset(&J[0][0], 0, sizeof(J));
681  for (i = 0; i < 5; i++) J[i][i] = 1.0;
682  }
683 
684  for (i = 0; i < 5; i++) for (j = 0; j < 5; j++) {
685  AG[i][j] = 0.0;
686  for (m = 0; m < 5; m++) AG[i][j] += J[i][m] * pTS->getTrackCovariance(m, j);
687  }
688  for (i = 0; i < 5; i++) for (j = i; j < 5; j++) {
689  Gf[i][j] = 0.0;
690  for (m = 0; m < 5; m++) Gf[i][j] += AG[i][m] * J[j][m];
691  Gf[j][i] = Gf[i][j];
692  }
693 
694  auto pTE = std::make_unique<Trk::TrkTrackState>(pTS);
695 
696  pTE->setTrackState(Rf);
697  pTE->setTrackCovariance(Gf);
698  pTE->attachToSurface(pSE);
699  pTE->applyMaterialEffects();
700 
701  for (i = 0; i < 5; i++) for (j = 0; j < 5; j++) {
702  Gi[i][j] = pTE->getTrackCovariance(i, j);
703  }
704 
705  matrixInversion5x5(Gi);
706 
707  for (i = 0; i < 5; i++) for (j = 0; j < 5; j++) {
708  A[i][j] = 0.0;
709  for (m = 0; m < 5; m++) A[i][j] += AG[m][i] * Gi[m][j];
710  }
711  pTE->setPreviousState(pTS);
712  pTE->setSmootherGain(A);
713 
714  return pTE;
715 }
716 
717 bool Trk::DistributedKalmanFilter::runForwardKalmanFilter(PVPNodes& pvpNodes,
718  PVPTrackStates& pvpTrackStates,
719  TrkTrackState* pInitState,
720  MagField::AtlasFieldCache& fieldCache) const {
721  bool OK = true;
722  Trk::TrkPlanarSurface* pSB = nullptr, * pSE;
723 
724  pvpTrackStates.push_back(std::make_unique<Trk::TrkTrackState>(pInitState));
725  Trk::TrkTrackState* pTS = pvpTrackStates.back().get();
726  for (std::unique_ptr<TrkBaseNode>& node : pvpNodes) {
727  pSE = node->getSurface();
728  std::unique_ptr<Trk::TrkTrackState> pNS = extrapolate(pTS, pSB, pSE, fieldCache);
729  pSB = pSE;
730  if (pNS) {
731  Trk::TrkTrackState* state = pNS.get();
732  pvpTrackStates.push_back(std::move(pNS));
733 
734  bool isNaN = std::isnan(pTS->getTrackState(0)) ||
735  std::isnan(pTS->getTrackState(1)) ||
736  std::isnan(pTS->getTrackState(2)) ||
737  std::isnan(pTS->getTrackState(3)) ||
738  std::isnan(pTS->getTrackState(4));
739  if (isNaN) {
740  OK = false;
741  break;
742  }
743 
744  node->validateMeasurement(state);
745  node->updateTrackState(state);
746  pTS = state;
747  } else {
748  OK = false;
749  break;
750  }
751  }
752  return OK;
753 }
754 
755 void Trk::DistributedKalmanFilter::runSmoother(PVPTrackStates& pvpTrackStates) {
756  for (std::unique_ptr<TrkTrackState>& state : pvpTrackStates) {
757  state->runSmoother();
758  }
759 }
760 
761 void Trk::DistributedKalmanFilter::calculateLRsolution(PVPNodes& pvpNodes) {
762  for (std::unique_ptr<TrkBaseNode>& node : pvpNodes) {
763  node->updateInternal();
764  }
765 }
766 
767 int Trk::DistributedKalmanFilter::findOutliers(PVPNodes& pvpNodes, double cut) {
768  double dchi2;
769  int nOutl = 0;
770 
771  for (std::unique_ptr<TrkBaseNode>& node : pvpNodes) {
772  dchi2 = (node->getChi2Distance(node->getTrackState())) / node->getNdof();
773  if (dchi2 > cut) {
774  if (node->isValidated()) nOutl++;
775  node->setNodeState(0);
776  } else node->setNodeState(1);
777  }
778  return nOutl;
779 }
780 
781 Trk::TrackStateOnSurface* Trk::DistributedKalmanFilter::createTrackStateOnSurface(Trk::TrkBaseNode* pN) const {
782  TrackStateOnSurface* pTSS = nullptr;
783  char type = pN->getNodeType();
784 
785  std::unique_ptr<Trk::TrackParameters> pTP {};
786  if (type == 0) return pTSS;
787 
788  TrkTrackState* pTS = pN->getTrackState();
789  const Trk::PrepRawData* pPRD = pN->getPrepRawData();
790 
791  if ((type == 1) || (type == 2)) {
792  const Trk::Surface& rS = pPRD->detectorElement()->surface();
793  const Trk::PlaneSurface* pPS = dynamic_cast<const Trk::PlaneSurface*>(&rS);
794  if (pPS == nullptr) return pTSS;
795 
796  AmgSymMatrix(5) pM;
797  ;
798  for (int i = 0; i < 5; i++) {
799  for (int j = 0; j < 5; j++) {
800  pM(i, j) = pTS->getTrackCovariance(i, j);
801  }
802  }
803  pTP = std::make_unique<Trk::AtaPlane>(pTS->getTrackState(0),
804  pTS->getTrackState(1),
805  pTS->getTrackState(2),
806  pTS->getTrackState(3),
807  pTS->getTrackState(4), *pPS,
808  std::move(pM));
809  } else if (type == 3) {
810  const Trk::Surface& rS = pPRD->detectorElement()->surface(pPRD->identify());
811  const Trk::StraightLineSurface* pLS = dynamic_cast<const Trk::StraightLineSurface*>(&rS);
812  if (pLS == nullptr) return pTSS;
813 
814  AmgSymMatrix(5) pM;
815  for (int i = 0; i < 5; i++) {
816  for (int j = 0; j < 5; j++) {
817  (pM) (i, j) = pTS->getTrackCovariance(i, j);
818  }
819  }
820  if ((pTS->getTrackState(2) < -M_PI) || (pTS->getTrackState(2) > M_PI)) printf("Phi is beyond the range\n");
821  pTP = std::make_unique<Trk::AtaStraightLine>(pTS->getTrackState(0),
822  pTS->getTrackState(1),
823  pTS->getTrackState(2),
824  pTS->getTrackState(3),
825  pTS->getTrackState(4),
826  *pLS,
827  std::move(pM));
828  }
829  if (pTP == nullptr) return nullptr;
830 
831  auto pRIO = std::unique_ptr<Trk::RIO_OnTrack>(m_ROTcreator->correct(*pPRD, *pTP));
832  if (pRIO == nullptr) {
833  return nullptr;
834  }
835  auto pFQ = Trk::FitQualityOnSurface(pN->getChi2(), pN->getNdof());
836  pTSS = new Trk::TrackStateOnSurface(pFQ, std::move(pRIO), std::move(pTP));
837  return pTSS;
838 }
839 
840 // fit a set of PrepRawData objects
841 std::unique_ptr<Trk::Track>
842 Trk::DistributedKalmanFilter::fit(
843  const EventContext& ctx,
844  const Trk::PrepRawDataSet& prepRDColl,
845  const Trk::TrackParameters& estimatedParametersNearOrigine,
846  const RunOutlierRemoval runOutlier,
847  const Trk::ParticleHypothesis matEffects) const {
848  const double outlCut = 12.0;
849  const double rlSili = 0.022;
850  const double rlTrt = 0.005;
851 
852  TrkPlanarSurface* pS;
853  TrkTrackState* pTS, * pInitState;
854  double Rk[5], C[3], N[3], M[3][3];
855  int i, nAmbHits = 0;
856  const Perigee* pP;
857 
858  Trk::Track* fittedTrack = nullptr;
859  Eigen::Vector3d mx, my, mz;
860  const Trk::PerigeeSurface perSurf;
861 
862  bool badTrack = false;
863 
864 
865  // protection, if empty PrepRawDataSet
866  if (prepRDColl.empty()) {
867  ATH_MSG_ERROR("try to fit empty vec<PrepRawData>, reject fit");
868  return nullptr;
869  }
870  if (runOutlier) ATH_MSG_VERBOSE("-> Outlier removal switched on");
871  if (matEffects != nonInteracting) ATH_MSG_VERBOSE("-> Material Effects switched on");
872 
873  // 0. Sort PRD set
874 
875  Trk::PrepRawDataComparisonFunction* compareHits = new Trk::PrepRawDataComparisonFunction
876  (estimatedParametersNearOrigine.position(),
877  estimatedParametersNearOrigine.momentum());
878  PrepRawDataSet inputPRDColl = PrepRawDataSet(prepRDColl);
879  if (!is_sorted(inputPRDColl.begin(), inputPRDColl.end(), *compareHits)) std::sort(
880  inputPRDColl.begin(), inputPRDColl.end(), *compareHits);
881  delete compareHits;
882  ATH_MSG_VERBOSE(" List of sorted PRDs: ");
883  for (const auto *prdIt : inputPRDColl) {
884  if (!prdIt->detectorElement()) continue;
885  ATH_MSG_VERBOSE(m_idHelper->print_to_string(prdIt->identify())
886  << " radius= "
887  << prdIt
888  ->detectorElement()
889  ->surface(prdIt->identify())
890  .center()
891  .mag());
892  }
893 
894  // 1. Create initial track state:
895 
896  pP = dynamic_cast<const Perigee*>(&estimatedParametersNearOrigine);
897  if (pP == nullptr) {
898  ATH_MSG_DEBUG(" fit needs perigee parameters - creating it");
899  const Trk::TrackParameters* pTP =
900  m_extrapolator->extrapolateDirectly(ctx,
901  estimatedParametersNearOrigine,
902  perSurf,
903  Trk::anyDirection,
904  false,
905  Trk::nonInteracting).release();
906 
907  pP = dynamic_cast<const Perigee*>(pTP);
908  if (pP == nullptr) {
909  ATH_MSG_ERROR("failed to create perigee parameters, reject fit");
910  return nullptr;
911  }
912  }
913  Rk[0] = pP->parameters()[Trk::d0];
914  Rk[1] = pP->parameters()[Trk::z0];
915  Rk[2] = pP->parameters()[Trk::phi0];
916  Rk[3] = pP->parameters()[Trk::theta];
917  Rk[4] = pP->parameters()[Trk::qOverP];
918  pTS = new TrkTrackState(Rk);
919 
920  if (matEffects == Trk::nonInteracting) {
921  pTS->setScatteringMode(0);
922  ATH_MSG_VERBOSE("-> Multiple Scattering treatment switched off");
923  } else {
924  if ((matEffects == Trk::pion) || (matEffects == Trk::muon)) {
925  pTS->setScatteringMode(1);
926  ATH_MSG_VERBOSE("-> Multiple Scattering treatment switched on");
927  } else {
928  if (matEffects == Trk::electron) {
929  pTS->setScatteringMode(2);
930  ATH_MSG_VERBOSE(
931  "-> Multiple Scattering and Bremm treatments switched on");
932  } else ATH_MSG_WARNING("Material setting " << matEffects << "not supported !");
933  }
934  }
935 
936  pInitState = pTS;
937 
938  // 2. Create filtering nodes:
939 
940  PVPNodes pvpNodes;
941  PVPSurfaces pvpSurfaces;
942  PVPTrackStates pvpTrackStates;
943 
944  for (const auto *prdIt : inputPRDColl) {
945  if (prdIt->detectorElement() == nullptr) {
946  ATH_MSG_WARNING(" PrepRawData has no detector element - drop it");
947  continue;
948  }
949  Identifier ID = prdIt->identify();
950  if (m_idHelper->is_indet(ID) || m_idHelper->is_sct(ID) ||
951  m_idHelper->is_pixel(ID) || m_idHelper->is_trt(ID)) {
952  if (m_idHelper->is_pixel(ID) || m_idHelper->is_sct(ID)) {
953  const Trk::Surface& rSurf = prdIt->detectorElement()->surface();
954  N[0] = rSurf.normal().x();
955  N[1] = rSurf.normal().y();
956  N[2] = rSurf.normal().z();
957  C[0] = rSurf.center().x();
958  C[1] = rSurf.center().y();
959  C[2] = rSurf.center().z();
960  mx = rSurf.transform().rotation().block(0, 0, 3, 1);
961  my = rSurf.transform().rotation().block(0, 1, 3, 1);
962  mz = rSurf.transform().rotation().block(0, 2, 3, 1);
963  for (i = 0; i < 3; i++) {
964  M[i][0] = mx[i];
965  M[i][1] = my[i];
966  M[i][2] = mz[i];
967  }
968  pvpSurfaces.push_back(std::make_unique<TrkPlanarSurface>(C, N, M, rlSili));
969  pS = pvpSurfaces.back().get();
970  if (m_idHelper->is_pixel(ID)) {
971  pvpNodes.push_back(std::make_unique<TrkPixelNode>(pS, 200.0, prdIt));
972  } else {
973  if (fabs(N[2]) < 0.1) {
974  pvpNodes.push_back(std::make_unique<TrkClusterNode>(pS, 200.0, prdIt));
975  } else {
976  pvpNodes.push_back(std::make_unique<TrkEndCapClusterNode>(pS, 200.0, prdIt));
977  }
978  }
979  continue;
980  }
981  if (m_idHelper->is_trt(ID)) {
982  nAmbHits++;
983 
984  const Trk::Surface& rSurf = prdIt->detectorElement()->surface(ID);
985  double len = 500.0;
986  const Trk::SurfaceBounds& rBounds =
987  prdIt->detectorElement()->surface().bounds();
988 
989  C[0] = rSurf.center().x();
990  C[1] = rSurf.center().y();
991  C[2] = rSurf.center().z();
992  double cosFs = C[0] / sqrt(C[0] * C[0] + C[1] * C[1]);
993  double sinFs = C[1] / sqrt(C[0] * C[0] + C[1] * C[1]);
994 
995  const Amg::Vector3D& rNorm =
996  prdIt->detectorElement()->surface().normal();
997  bool isBarrel = true;
998  if (fabs(rNorm.z()) > 0.2) isBarrel = false;
999 
1000  if (isBarrel) {
1001  const Trk::RectangleBounds& bBounds =
1002  dynamic_cast<const Trk::RectangleBounds&>(rBounds);
1003  len = bBounds.halflengthY();
1004  N[0] = cosFs;
1005  N[1] = sinFs;
1006  N[2] = 0.0;
1007 
1008  if (C[2] < 0.0) {
1009  M[0][0] = -sinFs;
1010  M[0][1] = 0.0;
1011  M[0][2] = cosFs;
1012  M[1][0] = cosFs;
1013  M[1][1] = 0.0;
1014  M[1][2] = sinFs;
1015  M[2][0] = 0.0;
1016  M[2][1] = 1.0;
1017  M[2][2] = 0.0;
1018  } else {
1019  M[0][0] = sinFs;
1020  M[0][1] = 0.0;
1021  M[0][2] = cosFs;
1022  M[1][0] = -cosFs;
1023  M[1][1] = 0.0;
1024  M[1][2] = sinFs;
1025  M[2][0] = 0.0;
1026  M[2][1] = -1.0;
1027  M[2][2] = 0.0;
1028  }
1029  } else {
1030  const Trk::DiscBounds& ecBounds =
1031  dynamic_cast<const Trk::DiscBounds&>(rBounds);
1032  len = 0.5 * (ecBounds.rMax() - ecBounds.rMin());
1033  N[0] = 0.0;
1034  N[1] = 0.0;
1035  N[2] = 1.0;
1036  M[0][0] = -sinFs;
1037  M[0][1] = -cosFs;
1038  M[0][2] = 0.0;
1039  M[1][0] = cosFs;
1040  M[1][1] = -sinFs;
1041  M[1][2] = 0.0;
1042  M[2][0] = 0.0;
1043  M[2][1] = 0.0;
1044  M[2][2] = 1.0;
1045  }
1046  pvpSurfaces.push_back(std::make_unique<TrkPlanarSurface>(C, N, M, rlTrt));
1047  pS = pvpSurfaces.back().get();
1048  pvpNodes.push_back(std::make_unique<TrkTrtNode>(pS, 200.0, -len, len, prdIt));
1049 
1050  continue;
1051  }
1052  }
1053  ATH_MSG_WARNING(" PrepRawData is neither identified nor supported !");
1054  ATH_MSG_WARNING("RIO ID prints as " << m_idHelper->print_to_string(ID));
1055  }
1056  ATH_MSG_VERBOSE(" Number of nodes/surfaces created: "
1057  << pvpNodes.size() << "/" << pvpSurfaces.size());
1058 
1059  SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle {
1060  m_fieldCacheCondObjInputKey, ctx
1061  };
1062 
1063  if (!readHandle.isValid()) {
1064  std::stringstream msg;
1065  msg << "Failed to retrieve magmnetic field conditions data "
1066  << m_fieldCacheCondObjInputKey.key() << ".";
1067  throw std::runtime_error(msg.str());
1068  }
1069  const AtlasFieldCacheCondObj* fieldCondObj {
1070  *readHandle
1071  };
1072 
1073  MagField::AtlasFieldCache fieldCache;
1074  fieldCondObj->getInitializedCache(fieldCache);
1075 
1076  // 3. Sort vector of filtering nodes according to their distances from the
1077  // perigee point
1078 
1079  // m_sortFilteringNodes(pInitState);
1080 
1081  // 4. Main algorithm: filter and smoother (Rauch-Tung-Striebel)
1082 
1083  if (runForwardKalmanFilter(pvpNodes, pvpTrackStates, pInitState, fieldCache)) {
1084  runSmoother(pvpTrackStates);
1085  if (runOutlier) {
1086  int nOutl = findOutliers(pvpNodes, outlCut);
1087  ATH_MSG_VERBOSE(nOutl << " outliers removed ");
1088 
1089  if ((nOutl * 1.0 /pvpNodes.size()) > 0.3) {
1090  ATH_MSG_VERBOSE(" two many outliers - bad track ");
1091  badTrack = true;
1092  }
1093  if ((nOutl != 0) && (!badTrack)) {
1094  pvpTrackStates.clear();
1095  runForwardKalmanFilter(pvpNodes, pvpTrackStates, pInitState, fieldCache);
1096  runSmoother(pvpTrackStates);
1097  nOutl = findOutliers(pvpNodes, outlCut);
1098  ATH_MSG_VERBOSE(nOutl << " new outliers found ");
1099  }
1100  }
1101  if ((nAmbHits != 0) && (!badTrack)) {
1102  calculateLRsolution(pvpNodes);
1103  pvpTrackStates.clear();
1104  runForwardKalmanFilter(pvpNodes, pvpTrackStates, pInitState, fieldCache);
1105  runSmoother(pvpTrackStates);
1106  }
1107 
1108  // 5. Create and store back all the stuff
1109  if (!badTrack) {
1110  pTS = pvpTrackStates.begin()->get();
1111  auto pMP {
1112  createMeasuredPerigee(pTS)
1113  };
1114  ATH_MSG_DEBUG("Fitted perigee: d0=" << pMP->parameters()[Trk::d0]
1115  << " z0=" << pMP->parameters()[Trk::z0]
1116  << " phi=" << pMP->parameters()[Trk::phi0]
1117  << " theta=" << pMP->parameters()[Trk::theta]
1118  << " qOverP=" << pMP->parameters()[Trk::qOverP]
1119  << " pT="
1120  << std::sin(pMP->parameters()[Trk::theta]) /
1121  pMP->parameters()[Trk::qOverP]);
1122 
1123  auto pvTS = std::make_unique<Trk::TrackStates>();
1124 
1125  std::bitset<Trk::TrackStateOnSurface::NumberOfTrackStateOnSurfaceTypes>
1126  typePattern;
1127  typePattern.set(Trk::TrackStateOnSurface::Perigee);
1128  const TrackStateOnSurface* pTSOS =
1129  new TrackStateOnSurface(nullptr, std::move(pMP), nullptr, typePattern);
1130 
1131  pvTS->push_back(pTSOS);
1132 
1133  double chi2 = 0.0;
1134  int ndof = -5;
1135 
1136  for (std::unique_ptr<TrkBaseNode>& node : pvpNodes) {
1137  if (node->isValidated()) {
1138  TrackStateOnSurface* pTSS = createTrackStateOnSurface(node.get());
1139  if (pTSS != nullptr) {
1140  pvTS->push_back(pTSS);
1141  chi2 += node->getChi2();
1142  ndof += node->getNdof();
1143  }
1144  }
1145  }
1146  ATH_MSG_DEBUG("Total Chi2: " << chi2 << " DoF=" << ndof);
1147  auto pFQ = std::make_unique<Trk::FitQuality>(chi2, ndof);
1148  ATH_MSG_DEBUG(pvTS->size() << " new RIO_OnTrack(s) created");
1149  TrackInfo info(TrackInfo::DistributedKalmanFilter, matEffects);
1150  fittedTrack = new Track(info, std::move(pvTS), std::move(pFQ));
1151  } else {
1152  fittedTrack = nullptr;
1153  }
1154  } else {
1155  ATH_MSG_WARNING("Extrapolation failed ");
1156  fittedTrack = nullptr;
1157  }
1158  delete pInitState;
1159  return std::unique_ptr<Trk::Track>(fittedTrack);
1160 }
1161 
1162 std::unique_ptr<Trk::Track>
1163 Trk::DistributedKalmanFilter::fit(
1164  const EventContext& ctx,
1165  const Trk::MeasurementSet& inputRotColl,
1166  const Trk::TrackParameters& estimatedParametersNearOrigine,
1167  const Trk::RunOutlierRemoval runOutlier,
1168  const Trk::ParticleHypothesis matEffects) const {
1169  if (inputRotColl.empty()) {
1170  ATH_MSG_ERROR("try to fit empty vec<MeasurementBase>, reject fit");
1171  return nullptr;
1172  }
1173 
1174  PrepRawDataSet prepRDColl;
1175  DataVector<const MeasurementBase>::const_iterator it = inputRotColl.begin();
1176  DataVector<const MeasurementBase>::const_iterator itEnd = inputRotColl.end();
1177  for (; it != itEnd; ++it) {
1178  if (!(*it)) {
1179  ATH_MSG_WARNING("This track contains empty MeasurementBase "
1180  << "objects! Skipped this MB..");
1181  } else {
1182  const RIO_OnTrack* pROT = dynamic_cast<const RIO_OnTrack*>(*it);
1183  const PrepRawData* prepRD = (pROT) ? (pROT->prepRawData()) : nullptr;
1184  if (!prepRD) {
1185  ATH_MSG_WARNING("RIO_OnTrack->prepRawRata() "
1186  << "returns no object, ignore... ");
1187  } else {
1188  prepRDColl.push_back(prepRD);
1189  }
1190  }
1191  }
1192  return fit(ctx,
1193  prepRDColl, estimatedParametersNearOrigine, runOutlier, matEffects);
1194 }
1195 
1196 std::unique_ptr<Trk::Track>
1197 Trk::DistributedKalmanFilter::fit(
1198  const EventContext& ctx,
1199  const Trk::Track& inputTrack,
1200  const Trk::MeasurementSet& inputRotColl,
1201  const Trk::RunOutlierRemoval runOutlier,
1202  const Trk::ParticleHypothesis matEffects) const {
1203  ATH_MSG_VERBOSE("--> enter DistributedKalmanFilter::fit(Track,easurementSet,)");
1204 
1205  // protection against track not having any parameters
1206  if (inputTrack.trackParameters()->empty()) {
1207  ATH_MSG_ERROR("need estimated track parameters near "
1208  << "origin, reject fit");
1209  return nullptr;
1210  }
1211  // protection against not having RIO_OnTracks
1212  if (inputTrack.measurementsOnTrack()->empty()) {
1213  ATH_MSG_ERROR("try to refit track with empty "
1214  << "vec<RIO_OnTrack>, reject fit");
1215  return nullptr;
1216  }
1217 
1218  // create PrepRawData subset
1219  PrepRawDataSet prepRDColl;
1220 
1221  // collect PrepRawData pointers from input track ROTs
1222  msg(MSG::VERBOSE) << "extract PrepRawDataSet from Track" << endmsg;
1223  DataVector<const MeasurementBase>::const_iterator it = inputTrack.measurementsOnTrack()->begin();
1224  DataVector<const MeasurementBase>::const_iterator itEnd = inputTrack.measurementsOnTrack()->end();
1225  for (; it != itEnd; ++it) {
1226  if (!(*it)) {
1227  ATH_MSG_WARNING("This track contains empty MeasurementBase "
1228  << "objects! Skipped this MB..");
1229  } else {
1230  const RIO_OnTrack* testROT = dynamic_cast<const RIO_OnTrack*>(*it);
1231  const PrepRawData* prepRD = (testROT) ? (testROT)->prepRawData() : nullptr;
1232  if (!prepRD) {
1233  ATH_MSG_WARNING("RIO_OnTrack->prepRawRata() "
1234  << "returns no object, ignore... ");
1235  } else {
1236  prepRDColl.push_back(prepRD);
1237  }
1238  }
1239  }
1240  const Trk::TrackParameters* minPar =
1241  *(std::min_element(inputTrack.trackParameters()->begin(),
1242  inputTrack.trackParameters()->end(),
1243  minTrkPar()));
1244 
1245  it = inputRotColl.begin();
1246  itEnd = inputRotColl.end();
1247  for (; it != itEnd; ++it) {
1248  if (!(*it)) {
1249  ATH_MSG_WARNING("This track contains empty MeasurementBase "
1250  << "objects! Skipped this MB..");
1251  } else {
1252  const RIO_OnTrack* pROT = dynamic_cast<const RIO_OnTrack*>(*it);
1253  const PrepRawData* prepRD = (pROT) ? (pROT->prepRawData()) : nullptr;
1254  if (!prepRD) {
1255  ATH_MSG_WARNING("RIO_OnTrack->prepRawRata() "
1256  << "returns no object, ignore... ");
1257  } else {
1258  prepRDColl.push_back(prepRD);
1259  }
1260  }
1261  }
1262  return fit(ctx, prepRDColl, *minPar, runOutlier, matEffects);
1263 }
1264 
1265 // fit a set of RIO_OnTrack objects
1266 std::unique_ptr<Trk::Track>
1267 Trk::DistributedKalmanFilter::fit(
1268  const EventContext& ctx,
1269  const Trk::RIO_OnTrackSet& inputRotColl,
1270  const Trk::TrackParameters& estimatedParametersNearOrigine,
1271  const Trk::RunOutlierRemoval runOutlier,
1272  const Trk::ParticleHypothesis matEffects) const {
1273  ATH_MSG_VERBOSE("--> entering DistributedKalmanFilter::fit"
1274  << "(RIO_OnTrackSet,TrackParameters,,)");
1275 
1276  // protection, if empty RIO_OnTrackSet
1277  if (inputRotColl.empty()) {
1278  ATH_MSG_ERROR("try to fit track+vec<ROT> with empty "
1279  << "extended vec<RIO_OnTrack>, reject fit");
1280  return nullptr;
1281  }
1282 
1283  if (runOutlier) ATH_MSG_VERBOSE("-> Outlier removal switched on");
1284  if (matEffects != Trk::nonInteracting) ATH_MSG_VERBOSE("-> Material Effects switched on");
1285 
1286  // for the time being, disintegrate ROTSet back to PrepRawData set.
1287  ATH_MSG_VERBOSE("copy & downgrade ROTSet into PRDset - "
1288  << "this is improvised.");
1289  PrepRawDataSet prepRDColl;
1290 
1291  RIO_OnTrackSet::const_iterator it = inputRotColl.begin();
1292  RIO_OnTrackSet::const_iterator itEnd = inputRotColl.end();
1293 
1294  for (; it != itEnd; ++it) {
1295  const PrepRawData* prepRD = ((*it)->prepRawData());
1296  if (!prepRD) {
1297  ATH_MSG_WARNING("RIO_OnTrack->prepRawRata() "
1298  << "returns no object, ignore... ");
1299  } else {
1300  prepRDColl.push_back(prepRD);
1301  }
1302  }
1303 
1304  // fit set of PrepRawData using main method, start with first Track Parameter in inputTrack
1305  ATH_MSG_VERBOSE("call fit(PRDSet,TP,,)");
1306  return fit(
1307  ctx, prepRDColl, estimatedParametersNearOrigine, runOutlier, matEffects);
1308 }
1309 
grepfile.info
info
Definition: grepfile.py:38
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
Trk::anyDirection
@ anyDirection
Definition: PropDirection.h:22
IDTPM::ndof
float ndof(const U &p)
Definition: TrackParametersHelper.h:142
Trk::RectangleBounds
Definition: RectangleBounds.h:38
cost
int cost(std::vector< std::string > &files, node &n, const std::string &directory="", bool deleteref=false, bool relocate=false)
Definition: hcg.cxx:921
TrkDetElementBase.h
Trk::RIO_OnTrackSet
std::vector< const RIO_OnTrack * > RIO_OnTrackSet
vector of detector hits on a track
Definition: FitterTypes.h:34
Trk::TrackInfo
Contains information about the 'fitter' of this track.
Definition: Tracking/TrkEvent/TrkTrack/TrkTrack/TrackInfo.h:32
Trk::TrackStateOnSurface::Perigee
@ Perigee
This represents a perigee, and so will contain a Perigee object only.
Definition: TrackStateOnSurface.h:117
checkxAOD.ds
ds
Definition: Tools/PyUtils/bin/checkxAOD.py:257
DataModel_detail::const_iterator
Const iterator class for DataVector/DataList.
Definition: DVLIterator.h:82
Trk::TrkPlanarSurface::transformPointToLocal
void transformPointToLocal(const double *, double *)
Definition: TrkPlanarSurface.cxx:113
Trk::PrepRawDataSet
std::vector< const PrepRawData * > PrepRawDataSet
vector of clusters and drift circles
Definition: FitterTypes.h:26
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
dqm_algorithms::tools::findOutliers
void findOutliers(std::vector< binContainer > &input, double &mean, double &scale, int &nIn, int nIterations, double exponent, double threshold, double SBCF=1., double nStop=8.)
Definition: AlgorithmHelper.cxx:888
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
athena.path
path
python interpreter configuration --------------------------------------—
Definition: athena.py:126
DiscBounds.h
TrackParameters.h
MeasurementBase.h
fitman.my
my
Definition: fitman.py:523
SG::ReadCondHandle
Definition: ReadCondHandle.h:44
ID
std::vector< Identifier > ID
Definition: CalibHitIDCheck.h:24
RectangleBounds.h
Surface.h
Trk::TrkPlanarSurface::transformPointToGlobal
void transformPointToGlobal(const double *, double *)
Definition: TrkPlanarSurface.cxx:122
Trk::Track
The ATLAS Track class.
Definition: Tracking/TrkEvent/TrkTrack/TrkTrack/Track.h:73
AtlasFieldCacheCondObj
Definition: AtlasFieldCacheCondObj.h:19
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
DistributedKalmanFilter.h
Trk::DistributedKalmanFilter::createTrackStateOnSurface
TrackStateOnSurface * createTrackStateOnSurface(TrkBaseNode *) const
Definition: DistributedKalmanFilter.cxx:781
Trk::PerigeeSurface
Definition: PerigeeSurface.h:43
Trk::ParametersBase::position
const Amg::Vector3D & position() const
Access method for the position.
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
Trk::ParametersT
Dummy class used to allow special convertors to be called for surfaces owned by a detector element.
Definition: EMErrorDetail.h:25
Trk::SurfaceBounds
Definition: SurfaceBounds.h:47
Trk::TrkTrackState
Definition: Tracking/TrkFitter/TrkDistributedKalmanFilter/TrkDistributedKalmanFilter/TrkTrackState.h:24
xAOD::JetInput::Track
@ Track
Definition: JetContainerInfo.h:61
Trk::TrkPlanarSurface::getRotMatrix
double getRotMatrix(int, int)
Definition: TrkPlanarSurface.cxx:84
Trk::TrkPlanarSurface
Definition: Tracking/TrkFitter/TrkDistributedKalmanFilter/TrkDistributedKalmanFilter/TrkPlanarSurface.h:25
skel.it
it
Definition: skel.GENtoEVGEN.py:423
Trk::DistributedKalmanFilter::calculateLRsolution
static void calculateLRsolution(PVPNodes &pvpNodes)
Definition: DistributedKalmanFilter.cxx:761
M_PI
#define M_PI
Definition: ActiveFraction.h:11
minTrkPar
Definition: DistributedKalmanFilter.cxx:168
Trk::one
@ one
Definition: TrkDetDescr/TrkSurfaces/TrkSurfaces/RealQuadraticEquation.h:22
Trk::DistributedKalmanFilter::initialize
virtual StatusCode initialize() override
Definition: DistributedKalmanFilter.cxx:146
mc.diff
diff
Definition: mc.SFGenPy8_MuMu_DD.py:14
Trk::z0
@ z0
Definition: ParamDefs.h:70
PropDirection.h
fitman.mx
mx
Definition: fitman.py:520
FitQualityOnSurface.h
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:158
Trk::RIO_OnTrack
Definition: RIO_OnTrack.h:70
JetTiledMap::N
@ N
Definition: TiledEtaPhiMap.h:44
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
Trk::TrkBaseNode::getNodeType
virtual char getNodeType()
Definition: TrkBaseNode.cxx:61
Trk::Surface::center
const Amg::Vector3D & center() const
Returns the center position of the Surface.
Trk::DistributedKalmanFilter::m_option_sortingRefPoint
std::vector< double > m_option_sortingRefPoint
Definition: DistributedKalmanFilter.h:180
Trk::DiscBounds::rMax
double rMax() const
This method returns outer radius.
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition: EventPrimitives.h:52
PrepRawData.h
Trk::DistributedKalmanFilter::DistributedKalmanFilter
DistributedKalmanFilter(const std::string &, const std::string &, const IInterface *)
Definition: DistributedKalmanFilter.cxx:131
Trk::RunOutlierRemoval
bool RunOutlierRemoval
switch to toggle quality processing after fit
Definition: FitterTypes.h:22
Trk::TrkBaseNode::getNdof
int getNdof() const
Definition: TrkBaseNode.cxx:65
Track.h
Trk::DistributedKalmanFilter::integrate
double integrate(double Rk[5], TrkPlanarSurface *pSB, TrkPlanarSurface *pSE, double *Rf, MagField::AtlasFieldCache &fieldCache) const
Definition: DistributedKalmanFilter.cxx:225
Trk::AmgSymMatrix
AmgSymMatrix(5) &GXFTrackState
Definition: GXFTrackState.h:156
Trk::FitQualityOnSurface
Definition: FitQualityOnSurface.h:19
Trk::ParticleHypothesis
ParticleHypothesis
Definition: ParticleHypothesis.h:25
AtlasDetectorID.h
This class provides an interface to generate or decode an identifier for the upper levels of the dete...
Trk::DiscBounds::rMin
double rMin() const
This method returns inner radius.
TrkTrackState.h
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
Trk::two
@ two
Definition: TrkDetDescr/TrkSurfaces/TrkSurfaces/RealQuadraticEquation.h:23
Trk::DistributedKalmanFilter::finalize
virtual StatusCode finalize() override
Definition: DistributedKalmanFilter.cxx:160
lumiFormat.i
int i
Definition: lumiFormat.py:92
Trk::DistributedKalmanFilter::runForwardKalmanFilter
bool runForwardKalmanFilter(PVPNodes &pvpNodes, PVPTrackStates &pvpTrackStates, TrkTrackState *, MagField::AtlasFieldCache &fieldCache) const
Definition: DistributedKalmanFilter.cxx:717
Identifier
Definition: DetectorDescription/Identifier/Identifier/Identifier.h:32
beamspotman.n
n
Definition: beamspotman.py:731
TrkBaseNode.h
Trk::theta
@ theta
Definition: ParamDefs.h:72
Trk::TrkBaseNode
Definition: Tracking/TrkFitter/TrkDistributedKalmanFilter/TrkDistributedKalmanFilter/TrkBaseNode.h:23
Trk::TrkTrackState::setScatteringMode
void setScatteringMode(int)
Definition: TrkTrackState.cxx:130
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
Trk::electron
@ electron
Definition: ParticleHypothesis.h:27
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
fitman.mz
mz
Definition: fitman.py:526
Trk::pion
@ pion
Definition: ParticleHypothesis.h:29
PrepRawDataComparisonFunction.h
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:522
BindingsTest.cut
cut
This script demonstrates how to call a C++ class from Python Also how to use PyROOT is shown.
Definition: BindingsTest.py:13
Trk::DistributedKalmanFilter::findOutliers
static int findOutliers(PVPNodes &pvpNodes, double)
Definition: DistributedKalmanFilter.cxx:767
Trk::TrkDetElementBase::surface
virtual const Surface & surface() const =0
Return surface associated with this detector element.
Trk::Surface::normal
virtual const Amg::Vector3D & normal() const
Returns the normal vector of the Surface (i.e.
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Trk::TrkPlanarSurface::getPar
double getPar(int)
Definition: TrkPlanarSurface.cxx:80
python.TransformConfig.descr
descr
print "%s.properties()" % self.__name__
Definition: TransformConfig.py:360
Trk::ParametersBase
Definition: ParametersBase.h:55
Trk::TrkPlanarSurface::getInvRotMatrix
double getInvRotMatrix(int, int)
Definition: TrkPlanarSurface.cxx:88
Trk::muon
@ muon
Definition: ParticleHypothesis.h:28
Trk::TrkBaseNode::getTrackState
TrkTrackState * getTrackState()
Definition: TrkBaseNode.cxx:43
ParticleHypothesis.h
MuonValidation_CreateResolutionProfiles.fit
def fit(h, emin, emax)
Definition: MuonValidation_CreateResolutionProfiles.py:69
Trk::MeasurementSet
std::vector< const MeasurementBase * > MeasurementSet
vector of fittable measurements
Definition: FitterTypes.h:30
Trk::PrepRawData
Definition: PrepRawData.h:62
Trk::PrepRawData::identify
Identifier identify() const
return the identifier
Trk::Track::trackParameters
const DataVector< const TrackParameters > * trackParameters() const
Return a pointer to a vector of TrackParameters.
Definition: Tracking/TrkEvent/TrkTrack/src/Track.cxx:97
python.PyKernel.detStore
detStore
Definition: PyKernel.py:41
Trk::TrackStateOnSurface
represents the track state (measurement, material, fit parameters and quality) at a surface.
Definition: TrackStateOnSurface.h:71
Trk::d0
@ d0
Definition: ParamDefs.h:69
TrackInfo.h
RIO_OnTrack.h
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
Trk::nonInteracting
@ nonInteracting
Definition: ParticleHypothesis.h:25
Trk::RectangleBounds::halflengthY
double halflengthY() const
for consitant naming
isNaN
bool isNaN(const T &t)
Definition: FlexErrArrayGroup.h:36
TrkFilteringNodes.h
Trk::DistributedKalmanFilter::~DistributedKalmanFilter
virtual ~DistributedKalmanFilter()
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Trk::RIO_OnTrack::prepRawData
virtual const Trk::PrepRawData * prepRawData() const =0
returns the PrepRawData (also known as RIO) object to which this RIO_OnTrack is associated.
Trk::DistributedKalmanFilter::runSmoother
static void runSmoother(PVPTrackStates &pvpTrackStates)
Definition: DistributedKalmanFilter.cxx:755
Trk::Track::measurementsOnTrack
const DataVector< const MeasurementBase > * measurementsOnTrack() const
return a pointer to a vector of MeasurementBase (NOT including any that come from outliers).
Definition: Tracking/TrkEvent/TrkTrack/src/Track.cxx:178
Trk::ParametersBase::momentum
const Amg::Vector3D & momentum() const
Access method for the momentum.
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
a
TList * a
Definition: liststreamerinfos.cxx:10
python.compareTCTs.ratio
ratio
Definition: compareTCTs.py:295
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Trk::PlaneSurface
Definition: PlaneSurface.h:64
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
Trk::TrackInfo::DistributedKalmanFilter
@ DistributedKalmanFilter
Fast Kalman filter from HLT with simplified material effects.
Definition: Tracking/TrkEvent/TrkTrack/TrkTrack/TrackInfo.h:59
MagField::AtlasFieldCache
Local cache for magnetic field (based on MagFieldServices/AtlasFieldSvcTLS.h)
Definition: AtlasFieldCache.h:43
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:73
python.LArCondContChannels.isBarrel
isBarrel
Definition: LArCondContChannels.py:659
Trk::DistributedKalmanFilter::fit
virtual std::unique_ptr< Track > fit(const EventContext &ctx, const Track &, const RunOutlierRemoval runOutlier=false, const ParticleHypothesis matEffects=Trk::nonInteracting) const override
RE-FIT A TRACK.
Definition: DistributedKalmanFilter.cxx:176
MuonHough::extrapolate
float extrapolate(const MuonLayerHough::Maximum &ref, const MuonLayerHough::Maximum &ex, bool doparabolic=false)
Definition: MuonLayerHough.cxx:519
Trk::TrkBaseNode::getChi2
double getChi2() const
Definition: TrkBaseNode.cxx:69
Trk::PrepRawDataComparisonFunction
Class providing comparison function, or relational definition, for PrepRawData.
Definition: PrepRawDataComparisonFunction.h:33
Trk::TrkTrackState::getTrackState
double getTrackState(int i)
Definition: Tracking/TrkFitter/TrkDistributedKalmanFilter/TrkDistributedKalmanFilter/TrkTrackState.h:47
Trk::DistributedKalmanFilter::PVPTrackStates
std::vector< std::unique_ptr< TrkTrackState > > PVPTrackStates
Definition: DistributedKalmanFilter.h:122
Trk::TrkBaseNode::getPrepRawData
virtual const PrepRawData * getPrepRawData()
Definition: TrkBaseNode.cxx:39
Trk::DistributedKalmanFilter::PVPNodes
std::vector< std::unique_ptr< TrkBaseNode > > PVPNodes
Definition: DistributedKalmanFilter.h:120
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:14
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
AthAlgTool
Definition: AthAlgTool.h:26
FitQuality.h
Trk::Surface
Definition: Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:75
Trk::Surface::transform
const Amg::Transform3D & transform() const
Returns HepGeom::Transform3D by reference.
python.compressB64.c
def c
Definition: compressB64.py:93
TrkPlanarSurface.h
python.AutoConfigFlags.msg
msg
Definition: AutoConfigFlags.py:7
Trk::DistributedKalmanFilter::extrapolate
std::unique_ptr< Trk::TrkTrackState > extrapolate(TrkTrackState *, TrkPlanarSurface *, TrkPlanarSurface *, MagField::AtlasFieldCache &fieldCache) const
Definition: DistributedKalmanFilter.cxx:315
Trk::phi0
@ phi0
Definition: ParamDefs.h:71
TrackStateOnSurface.h
node
Definition: memory_hooks-stdcmalloc.h:74
Trk::TrkTrackState::getTrackCovariance
double getTrackCovariance(int i, int j)
Definition: Tracking/TrkFitter/TrkDistributedKalmanFilter/TrkDistributedKalmanFilter/TrkTrackState.h:51
Trk::StraightLineSurface
Definition: StraightLineSurface.h:51
Trk::DiscBounds
Definition: DiscBounds.h:44
fitman.k
k
Definition: fitman.py:528
Trk::DistributedKalmanFilter::PVPSurfaces
std::vector< std::unique_ptr< TrkPlanarSurface > > PVPSurfaces
Definition: DistributedKalmanFilter.h:121
Trk::PrepRawData::detectorElement
virtual const TrkDetElementBase * detectorElement() const =0
return the detector element corresponding to this PRD The pointer will be zero if the det el is not d...
Generated on Thu Jul 4 2024 21:14:09 for ATLAS Offline Software by   doxygen 1.8.18