ATLAS Offline Software
TrkVKalVrtFitterTestAlg.cxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration.
3  */
11 #include "EventPrimitives/EventPrimitivesHelpers.h"
12 #include "xAODTracking/Vertex.h"
13 #include "xAODTracking/TrackParticle.h"
14 #include "TrkTrack/Track.h"
15 #include "TrkParticleBase/TrackParticleBase.h"
16 //undef NDEBUG after EDM
17 #undef NDEBUG
18 #include "TrkVKalVrtFitter/VxCascadeInfo.h"
19 #include "TrkVKalVrtFitterTestAlg.h"
20 #include "TrkVKalVrtFitter/IVertexCascadeFitter.h"
21 #include "TestTools/FLOATassert.h"
22 #include "GaudiKernel/SystemOfUnits.h"
23 #include "TrkVKalVrtFitter/IVKalState.h" //unique_ptr needs the deleter
24 #include "AthContainers/ConstAccessor.h"
25 
26 #include <cassert>
27 #include <cmath>
28 
29 #include "CLHEP/Vector/LorentzVector.h"
30 
31 using Gaudi::Units::mm;
32 using Gaudi::Units::MeV;
33 using Gaudi::Units::GeV;
34 
35 
36 namespace {
37 
38 
39 template <class T>
40 std::vector<const T*> asVec (const std::vector<std::unique_ptr<T> >& v)
41 {
42  std::vector<const T*> ret;
43  ret.reserve(v.size());
44 
45  for (const std::unique_ptr<T>& p : v) {
46  ret.push_back (p.get());
47  }
48  return ret;
49 }
50 
51 
52 std::vector<const Trk::NeutralParameters*>
53 asVec (const std::vector<std::unique_ptr<Trk::NeutralPerigee> >& v)
54 {
55  std::vector<const Trk::NeutralParameters*> ret;
56  ret.reserve(v.size());
57 
58  for (const std::unique_ptr<Trk::NeutralPerigee>& p : v) {
59  ret.push_back (p.get());
60  }
61  return ret;
62 }
63 
64 
65 AmgSymMatrix(5) cov5()
66 {
67  AmgSymMatrix(5) m;
68  m.setIdentity();
69  return m;
70 }
71 
72 
73 using PerigeeUVec_t = std::vector<std::unique_ptr<Trk::Perigee> >;
74 
75 using NeutralUVec_t = std::vector<std::unique_ptr<Trk::NeutralPerigee> >;
76 NeutralUVec_t makeNeutrals1()
77 {
78  Amg::Vector3D pos0 { 0, 0, 0 };
79 
80  Amg::Vector3D pos1a { 3*mm, 0.5*mm, -7*mm };
81  Amg::Vector3D mom1a { 1000*MeV, 900*MeV, 2000*MeV };
82  Amg::Vector3D pos1b { -1*mm, 2.5*mm, 1*mm };
83  Amg::Vector3D mom1b { 800*MeV, 1000*MeV, 300*MeV };
84  Amg::Vector3D pos1c { -1.5*mm, 1*mm, -3*mm };
85  Amg::Vector3D mom1c { 500*MeV, 4000*MeV, 800*MeV };
86 
87  NeutralUVec_t ret;
88 
89  ret.emplace_back (std::make_unique<Trk::NeutralPerigee>(pos1a, mom1a, 1, pos1a, cov5()));
90  ret.emplace_back (std::make_unique<Trk::NeutralPerigee>(pos1b, mom1b, 1, pos1b, cov5()));
91  ret.emplace_back (std::make_unique<Trk::NeutralPerigee>(pos1c, mom1c, 1, pos1c, cov5()));
92 
93  return ret;
94 }
95 
96 
97 
98 
99 // Copied from TrackParticleCreatorTool.
100 void setDefiningParameters( xAOD::TrackParticle& tp,
101  const Trk::Perigee& perigee )
102 {
103  tp.setDefiningParameters(perigee.parameters()[Trk::d0],
104  perigee.parameters()[Trk::z0],
105  perigee.parameters()[Trk::phi0],
106  perigee.parameters()[Trk::theta],
107  perigee.parameters()[Trk::qOverP]);
108  const AmgSymMatrix(5)* covMatrix = perigee.covariance();
109  // see https://its.cern.ch/jira/browse/ATLASRECTS-645 for justification to comment out the following line
110  // assert(covMatrix && covMatrix->rows()==5&& covMatrix->cols()==5);
111  std::vector<float> covMatrixVec;
112  if( !covMatrix ) ;//ATH_MSG_WARNING("Setting Defining parameters without error matrix");
113  else Amg::compress(*covMatrix,covMatrixVec);
114  tp.setDefiningParametersCovMatrixVec(covMatrixVec);
115  const Amg::Vector3D& surfaceCenter = perigee.associatedSurface().center();
116  tp.setParametersOrigin(surfaceCenter.x(), surfaceCenter.y(), surfaceCenter.z() );
117 }
118 void setDefiningParameters( xAOD::NeutralParticle& tp,
119  const Trk::NeutralPerigee& perigee )
120 {
121  tp.setDefiningParameters(perigee.parameters()[Trk::d0],
122  perigee.parameters()[Trk::z0],
123  perigee.parameters()[Trk::phi0],
124  perigee.parameters()[Trk::theta],
125  perigee.parameters()[Trk::qOverP]);
126  const AmgSymMatrix(5)* covMatrix = perigee.covariance();
127  // see https://its.cern.ch/jira/browse/ATLASRECTS-645 for justification to comment out the following line
128  // assert(covMatrix && covMatrix->rows()==5&& covMatrix->cols()==5);
129  std::vector<float> covMatrixVec;
130  if( !covMatrix ) ;//ATH_MSG_WARNING("Setting Defining parameters without error matrix");
131  else Amg::compress(*covMatrix,covMatrixVec);
132  tp.setDefiningParametersCovMatrixVec(covMatrixVec);
133  const Amg::Vector3D& surfaceCenter = perigee.associatedSurface().center();
134  tp.setParametersOrigin(surfaceCenter.x(), surfaceCenter.y(), surfaceCenter.z() );
135 }
136 
137 
138 using xAODTPUVec_t = std::vector<std::unique_ptr<xAOD::TrackParticle> >;
139 xAODTPUVec_t makexAODTP (PerigeeUVec_t&& perigees)
140 {
141  xAODTPUVec_t tracks;
142 
143  for (std::unique_ptr<Trk::Perigee>& p : perigees) {
144  auto tp = std::make_unique<xAOD::TrackParticle>();
145  tp->makePrivateStore();
146  setDefiningParameters (*tp, *p);
147  tracks.push_back (std::move (tp));
148  }
149 
150  return tracks;
151 }
152 
153 
154 using xAODNPUVec_t = std::vector<std::unique_ptr<xAOD::NeutralParticle> >;
155 xAODNPUVec_t makexAODNP (NeutralUVec_t&& perigees)
156 {
157  xAODNPUVec_t tracks;
158 
159  for (std::unique_ptr<Trk::NeutralPerigee>& p : perigees) {
160  auto tp = std::make_unique<xAOD::NeutralParticle>();
161  tp->makePrivateStore();
162  setDefiningParameters (*tp, *p);
163  tracks.push_back (std::move (tp));
164  }
165 
166  return tracks;
167 }
168 
169 
170 void dumpCovariance (const AmgSymMatrix(5)& m)
171 {
172  for (int i=0; i < 5; i++) {
173  for (int j=0; j < 5; j++) {
174  std::cout << m(i,j) << ", ";
175  }
176  }
177 }
178 
179 
180 [[maybe_unused]]
181 void dumpVertex (const xAOD::Vertex& v)
182 {
183  std::cout << "vvv\n";
184  std::cout << v.x() << " " << v.y() << " " << v.z() << "\n";
185  static const SG::ConstAccessor<short> vertexTypeAcc ("vertexType");
186  if (vertexTypeAcc.isAvailable(v)) {
187  std::cout << "vertexType " << v.vertexType() << "\n";
188  }
189  std::cout << "chi2/ndof " << v.chiSquared() << " " << v.numberDoF() << "\n";
190 
191  std::cout << "cov ";
192  for (float f : v.covariance()) {
193  std::cout << f << " ";
194  }
195  std::cout << "\n";
196 
197  static const SG::ConstAccessor<std::vector<ElementLink<xAOD::TrackParticleContainer> > > trackParticleLinksAcc ("trackParticleLinks");
198  if (trackParticleLinksAcc.isAvailable(v)) {
199  std::cout << "tplinks ";
200  for (const ElementLink< xAOD::TrackParticleContainer >& l : v.trackParticleLinks()) {
201  std::cout << l.dataID() << "/" << l.index() << " ";
202  }
203  std::cout << "\n";
204  }
205 
206  static const SG::ConstAccessor<float> trackWeightsAcc ("trackWeights");
207  if (trackWeightsAcc.isAvailable(v)) {
208  std::cout << "wt ";
209  for (float f : v.trackWeights()) {
210  std::cout << f << " ";
211  }
212  std::cout << "\n";
213  }
214 
215  static const SG::ConstAccessor<std::vector<ElementLink<xAOD::NeutralParticleContainer> > > neutralParticleLinksAcc ("neutralParticleLinks");
216  if (neutralParticleLinksAcc.isAvailable(v)) {
217  std::cout << "nplinks ";
218  for (const ElementLink< xAOD::NeutralParticleContainer >& l : v.neutralParticleLinks()) {
219  std::cout << l.dataID() << "/" << l.index() << " ";
220  }
221  std::cout << "\n";
222  }
223 
224  static const SG::ConstAccessor<float> neutralWeightsAcc ("neutralWeights");
225  if (neutralWeightsAcc.isAvailable(v)) {
226  std::cout << "wt ";
227  for (float f : v.neutralWeights()) {
228  std::cout << f << " ";
229  }
230  std::cout << "\n";
231  }
232 
233  std::cout << v.vxTrackAtVertexAvailable() << "\n";
234  for (const Trk::VxTrackAtVertex& vv : v.vxTrackAtVertex()) {
235  vv.dump (std::cout);
236  std::cout << "cov ";
237  if (vv.perigeeAtVertex()) {
238  dumpCovariance (*vv.perigeeAtVertex()->covariance());
239  }
240  else {
241  std::cout << "(null)";
242  }
243  std::cout << "\n";
244  }
245 }
246 
247 
248 void assertVec3D (const char* which,
249  const Amg::Vector3D& a,
250  const Amg::Vector3D& b,
251  double thresh = 2e-5)
252 {
253  thresh = std::max (thresh, 2e-5);
254  if ( ! Athena_test::isEqual (a.x(), b.x(), thresh) ||
255  ! Athena_test::isEqual (a.y(), b.y(), thresh) ||
256  ! Athena_test::isEqual (a.z(), b.z(), thresh) )
257  {
258  std::cerr << "TrkVKalVrtFitterTestAlg::assertVec3D mismatch " << which
259  << ": ["
260  << a.x() << ", "
261  << a.y() << ", "
262  << a.z() << "] / ["
263  << b.x() << ", "
264  << b.y() << ", "
265  << b.z() << "]\n";
266  std::abort();
267  }
268 }
269 
270 
271 void comparePerigee (const Trk::TrackParameters* a,
272  const Trk::TrackParameters* b,
273  double thresh = 2e-5)
274 {
275  if (!a && !b) return;
276  if (!a || !b) std::abort();
277  assertVec3D ("perigee pos", a->position(), b->position(), thresh);
278  assertVec3D ("perigee mom", a->momentum(), b->momentum(), thresh);
279  assert (a->charge() == b->charge());
280  assertVec3D ("perigee surf",
281  a->associatedSurface().center(),
282  b->associatedSurface().center(), thresh);
283 }
284 
285 
286 void compareVertex (const xAOD::Vertex& a, const xAOD::Vertex& b,
287  double thresh = 2e-5)
288 {
289  assertVec3D ("vertex pos", a.position(), b.position(), thresh);
290  assert (Athena_test::isEqual (a.chiSquared(), b.chiSquared(),
291  std::max (5e-5, thresh)) );
292  assert (Athena_test::isEqual (a.numberDoF(), b.numberDoF(), 1e-5) );
293  assert (a.covariance().size() == b.covariance().size());
294  for (unsigned int i = 0; i < a.covariance().size(); i++) {
295  if (std::isinf(a.covariance()[i]) && std::isinf(b.covariance()[i])) continue;
296  assert (Athena_test::isEqual (a.covariance()[i], b.covariance()[i], 2e-2) );
297  }
298 
299  assert (a.vxTrackAtVertexAvailable() == b.vxTrackAtVertexAvailable());
300  if (a.vxTrackAtVertexAvailable()) {
301  const std::vector< Trk::VxTrackAtVertex >& avec = a.vxTrackAtVertex();
302  const std::vector< Trk::VxTrackAtVertex >& bvec = b.vxTrackAtVertex();
303  assert (avec.size() == bvec.size());
304  for (unsigned int i = 0; i < avec.size(); i++) {
305  comparePerigee (avec[i].initialPerigee(), bvec[i].initialPerigee(), thresh);
306  comparePerigee (avec[i].perigeeAtVertex(), bvec[i].perigeeAtVertex(), thresh);
307  assert (Athena_test::isEqual (avec[i].trackQuality().chiSquared(),
308  bvec[i].trackQuality().chiSquared(),
309  3e-2));
310  assert (avec[i].trackQuality().numberDoF() ==
311  bvec[i].trackQuality().numberDoF());
312  }
313  }
314 }
315 
316 
317 void setRefittedPerigee (xAOD::Vertex& v, unsigned i,
318  float charge,
319  const Amg::Vector3D& mom,
320  const std::vector<float>& c)
321 {
322  std::vector< Trk::VxTrackAtVertex >& vec = v.vxTrackAtVertex();
323  if (vec.size() <= i) vec.resize(i+1);
324 
325  AmgSymMatrix(5) cov = cov5();
326  for (int i=0; i < 5; i++) {
327  for (int j=0; j < 5; j++) {
328  unsigned ipos = i*5 + j;
329  (cov)(i,j) = ipos < c.size() ? c[ipos] : 0;
330  }
331  }
332 
333  const Amg::Vector3D& pos = v.position();
334  auto p = std::make_unique<Trk::Perigee> (pos, mom, charge, pos,
335  cov);
336  vec[i].setPerigeeAtVertex (p.release());
337 }
338 
339 
340 void setFitQuality (xAOD::Vertex& v, unsigned i, float chi2, int ndof)
341 {
342  std::vector< Trk::VxTrackAtVertex >& vec = v.vxTrackAtVertex();
343  if (vec.size() <= i) vec.resize(i+1);
344  vec[i].setTrackQuality (Trk::FitQuality (chi2, ndof));
345 }
346 
347 
348 } // anonymous namespace
349 
350 
351 namespace Trk {
352 
353 
357 StatusCode TrkVKalVrtFitterTestAlg::initialize()
358 {
359  ATH_CHECK( m_fitter.retrieve() );
360  return StatusCode::SUCCESS;
361 }
362 
363 
367 StatusCode TrkVKalVrtFitterTestAlg::execute()
368 {
369  ATH_MSG_VERBOSE ("execute");
370 
371  ATH_CHECK( test1() );
372  ATH_CHECK( test2() );
373  ATH_CHECK( test3() );
374 
375  return StatusCode::SUCCESS;
376 }
377 
378 
379 
380 // Neutral, no constraint.
381 // (Mixed charged+neutral seems not to work.)
382 StatusCode TrkVKalVrtFitterTestAlg::test1()
383 {
384  xAOD::Vertex exp_v0;
385  exp_v0.makePrivateStore();
386  exp_v0.setPosition ({-5.58116, -8.50767, -4.76804});
387  exp_v0.setFitQuality (0.541406, 3);
388  exp_v0.setCovariance(std::vector<float>{
389  259.128, 394.824, 764.131, 144.775, 266.345, 155.231,
390  -30152, -6184.95, -2642.5, 4.60831e+13, 10912.8, -4760.06,
391  6318.67, 1.0494e+14, 2.4226e+14, -12079, -23269.1, -40859.3,
392  2.00894e+13, 4.51754e+13, 9.12659e+12});
393  setFitQuality(exp_v0, 0, 0.311, 2);
394  setFitQuality (exp_v0, 1, 0.138, 2);
395  setFitQuality (exp_v0, 2, 0.082, 2);
396 
397  NeutralUVec_t neutrals = makeNeutrals1();
398  std::unique_ptr<xAOD::Vertex> v1 (m_fitter->fit (std::vector<const Trk::TrackParameters*>(),
399  asVec (neutrals)));
400  compareVertex (*v1, exp_v0);
401 
402  return StatusCode::SUCCESS;
403 }
404 
405 
406 
407 // Neutral + Vector3D constraint
408 StatusCode TrkVKalVrtFitterTestAlg::test2()
409 {
410  xAOD::Vertex exp_v0;
411  exp_v0.makePrivateStore();
412  exp_v0.setPosition ({ 24.131, 38.8186, 13.2978});
413  exp_v0.setFitQuality (0.564942, 3);
414  exp_v0.setCovariance(std::vector<float>{
415  21964.5, 32891.5, 51244.9, 13612.1, 20963.9,
416  9573.41, 816688, 1.09315e+06, 393680, 1.92433e+14,
417  1.10804e+06, 1.87392e+06, 616429, 3.20119e+14, 5.41674e+14,
418  401502, 640572, 338398, 1.97024e+14, 3.42617e+14,
419  2.25906e+14});
420  setFitQuality(exp_v0, 0, 0.306, 2);
421  setFitQuality(exp_v0, 1, 0.00660789, 2);
422  setFitQuality(exp_v0, 2, 0.266975, 2);
423 
424  Amg::Vector3D pnt1(5, 6, -3);
425 
426  NeutralUVec_t neutrals = makeNeutrals1();
427  std::unique_ptr<xAOD::Vertex> v1 (m_fitter->fit (std::vector<const Trk::TrackParameters*>(),
428  asVec (neutrals),
429  pnt1));
430  compareVertex (*v1, exp_v0);
431 
432  // ??? This gives a different result due to different reference frame
433  // handling in CvtNeutralParameters vs CvtNeutralParticle
434  xAOD::Vertex exp_v1;
435  exp_v1.makePrivateStore();
436  exp_v1.setPosition ({27.7411, 46.5526, 14.3721});
437  exp_v1.setFitQuality (0.537727, 3);
438  exp_v1.setCovariance(std::vector<float>{
439  32598.3, 50372.7, 80773.1, 18537.3, 29396.9,
440  12058.1, 850485, 1.1479e+06, 380180, 1.54295e+14,
441  1.21881e+06, 2.11418e+06, 643860, 2.55981e+14, 4.31589e+14,
442  417309, 683027, 356607, 1.51101e+14, 2.64837e+14,
443  1.77033e+14});
444  setFitQuality(exp_v1, 0, 0.281, 2);
445  setFitQuality(exp_v1, 1, 0.010, 2);
446  setFitQuality(exp_v1, 2, 0.248, 2);
447 
448  xAODNPUVec_t xaodnp = makexAODNP(makeNeutrals1());
449  std::unique_ptr<xAOD::Vertex> v2 (m_fitter->fit (std::vector<const xAOD::TrackParticle*>(),
450  asVec (xaodnp),
451  pnt1));
452  compareVertex (*v2, exp_v1);
453 
454  return StatusCode::SUCCESS;
455 }
456 
457 
458 // Neutral + Vertex constraint
459 StatusCode TrkVKalVrtFitterTestAlg::test3()
460 {
461  xAOD::Vertex exp_v0;
462  exp_v0.makePrivateStore();
463  exp_v0.setPosition ({5.05085, 6.19647, -2.81445});
464  exp_v0.setFitQuality (2.20220, 6);
465  exp_v0.setCovariance (std::vector<float>
466  { 0.955506, 0.0266418, 0.95407, -0.0110041,
467  -0.0101241, 0.953395, 530.063, -411.644,
468  -325.151, 3.03883e+14, -711.821, 1110.92,
469  -551.227, 2.09286e+14, 2.49946e+14, 123.813,
470  386.78, 892.698, -3.14264e+12, 1.04832e+14,
471  1.10053e+14 });
472  setFitQuality (exp_v0, 0, 0.396, 2);
473  setFitQuality (exp_v0, 1, 0.752, 2);
474  setFitQuality (exp_v0, 2, 1.055, 2);
475 
476  Amg::Vector3D pnt1 (5, 6, -3);
477  xAOD::Vertex pnt2;
478  pnt2.makePrivateStore();
479  pnt2.setPosition (pnt1);
480  AmgSymMatrix(3) pnt2covar;
481  pnt2covar.setIdentity();
482  pnt2.setCovariancePosition (pnt2covar);
483 
484  NeutralUVec_t neutrals = makeNeutrals1();
485  std::unique_ptr<xAOD::Vertex> v1 (m_fitter->fit (std::vector<const Trk::TrackParameters*>(),
486  asVec (neutrals),
487  pnt2));
488  compareVertex (*v1, exp_v0);
489 
490  // ??? This gives a different result due to different reference frame
491  // handling in CvtNeutralParameters vs CvtNeutralParticle
492  xAOD::Vertex exp_v1;
493  exp_v1.makePrivateStore();
494  exp_v1.setPosition ({5.00912, 6.23591, -2.82707});
495  exp_v1.setFitQuality (2.25019, 6);
496  exp_v1.setCovariance (std::vector<float>
497  { 0.968686, 0.0228402, 0.951705, -0.00910123,
498  -0.0098848, 0.950267, 435.6, -418.438,
499  -276.022, 3.17041e+14, -401.428, 987.01,
500  -336.323, 1.592e+14, 1.35435e+14, 183.46,
501  272.759, 988.474, -1.22294e+13, 5.57451e+13,
502  7.21632e+13 });
503  setFitQuality (exp_v1, 0, 0.313, 2);
504  setFitQuality (exp_v1, 1, 0.692, 2);
505  setFitQuality (exp_v1, 2, 1.246, 2);
506 
507  xAODNPUVec_t xaodnp = makexAODNP (makeNeutrals1());
508  std::unique_ptr<xAOD::Vertex> v2 (m_fitter->fit (std::vector<const xAOD::TrackParticle*>(),
509  asVec (xaodnp),
510  pnt2));
511  compareVertex (*v2, exp_v1);
512 
513  return StatusCode::SUCCESS;
514 }
515 
516 
517 } // namespace Trk
518 
519 
520 
521 namespace {
522 
523 
524 AmgSymMatrix(5) cov5a()
525 {
526  AmgSymMatrix(5) m;
527  m.setZero();
528  for (int i=0; i < 5; i++) {
529  (m)(i,i) = 1e-2;
530  }
531  (m)(1,1)=1;
532  return m;
533 }
534 
535 
536 using PerigeeUVec_t = std::vector<std::unique_ptr<Trk::Perigee> >;
537 PerigeeUVec_t makePerigees2()
538 {
539  Amg::Vector3D pos1a { 10*mm, 0*mm, -5*mm };
540  Amg::Vector3D mom1a { 1000*MeV, 0*MeV, 0*MeV };
541  Amg::Vector3D pos1b { 10.5*mm, 0.5*mm, -5.5*mm };
542  Amg::Vector3D mom1b { 800*MeV, 200*MeV, 200*MeV };
543  Amg::Vector3D pos1c { 9.5*mm, -0.5*mm, -4.5*mm };
544  Amg::Vector3D mom1c { 700*MeV, -300*MeV, -200*MeV };
545 
546  PerigeeUVec_t ret;
547 
548  ret.emplace_back (std::make_unique<Trk::Perigee>(pos1a, mom1a, 1, pos1a, cov5a()));
549  ret.emplace_back (std::make_unique<Trk::Perigee>(pos1b, mom1b, -1, pos1a, cov5a()));
550  ret.emplace_back (std::make_unique<Trk::Perigee>(pos1c, mom1c, -1, pos1a, cov5a()));
551 
552  return ret;
553 }
554 
555 
556 PerigeeUVec_t makePerigees3()
557 {
558  Amg::Vector3D pos1a { 5*mm, 0*mm, -2.5*mm };
559  Amg::Vector3D mom1a { 600*MeV, 400*MeV, 200*MeV };
560  Amg::Vector3D pos1b { 5.1*mm, 0.3*mm, -2.6*mm };
561  Amg::Vector3D mom1b { 700*MeV, -300*MeV, -150*MeV };
562 
563  PerigeeUVec_t ret;
564 
565  ret.emplace_back (std::make_unique<Trk::Perigee>(pos1a, mom1a, 1, pos1a, cov5a()));
566  ret.emplace_back (std::make_unique<Trk::Perigee>(pos1b, mom1b, -1, pos1a, cov5a()));
567 
568  return ret;
569 }
570 
571 
572 Amg::Vector3D exp_mom (const double moms[][4], int imom)
573 {
574  return {moms[imom][0], moms[imom][1], moms[imom][2]};
575 }
576 
577 
578 } // anonymous namespace
579 
580 
581 namespace Trk {
582 
583 
584 // Simple cascade fitter test.
585 StatusCode TrkVKalVrtFitterTestAlg::test4()
586 {
587  Trk::IVertexCascadeFitter* fitter =
588  dynamic_cast<Trk::IVertexCascadeFitter*> (m_fitter.get());
589 
590  xAODTPUVec_t tracks1 = makexAODTP (makePerigees2());
591 
592  std::unique_ptr<IVKalState> state (fitter->makeState());
593  Trk::VertexID v1 = fitter->startVertex (asVec (tracks1),
594  std::vector<double> {100*MeV, 150*MeV, 200*MeV},
595  *state,
596  930*MeV);
597 
598  xAODTPUVec_t tracks2 = makexAODTP (makePerigees3());
599 
600  fitter->nextVertex (asVec (tracks2),
601  std::vector<double> {130*MeV, 160*MeV},
602  std::vector<Trk::VertexID> {v1},
603  *state,
604  2000*MeV);
605 
606  std::unique_ptr<Trk::VxCascadeInfo> info1 (fitter->fitCascade(*state));
607  info1->setSVOwnership (true);
608 
609 #if 0
610  std::cout << info1->fitChi2() << " " << info1->nDoF() << "\n";
611  for (const std::vector<TLorentzVector>& vv : info1->getParticleMoms()) {
612  std::cout << "===\n";
613  for (const TLorentzVector& vvv : vv) {
614  std::cout << vvv.X() << " " << vvv.Y() << " " << vvv.Z() << " " << vvv.E() << "\n";
615  }
616  }
617  std::cout << "=== vertices\n";
618  for (const xAOD::Vertex* v : info1->vertices()) {
619  dumpVertex (*v);
620  }
621 #endif
622 
623  // Some of the comparison thresholds here have to be quite large,
624  // because the fit is unstable against small changes in the rounding
625  // of, eg, trig functions. The results we get can vary depending
626  // on the libc version used as well as on the exact hardware used.
627 
628  assert (Athena_test::isEqual (info1->fitChi2(), 8.69008, 1e-5));
629  assert (info1->nDoF() == 8);
630 
631 
632  const double exp_moms0[][4] =
633  {
634  {65.8357, -2.03326, -1.46405, 119.752},
635  {755.228, 239.515, 134.648, 817.537},
636  {900.997, -348.825, -292.857, 1029.19},
637  {1013.68, 681.319, 331.188, 1272.13},
638  {522.571, -222.398, -113.273, 600.81},
639  {1719.34, -112.854, -155.908, 1964.27},
640  };
641  const size_t nmoms0 = std::distance (std::begin(exp_moms0), std::end(exp_moms0));
642  size_t imoms0 = 0;
643  for (const std::vector<TLorentzVector>& vv : info1->getParticleMoms()) {
644  for (const TLorentzVector& vvv : vv) {
645  assert (imoms0 < nmoms0);
646  assert( Athena_test::isEqual (vvv.X(), exp_moms0[imoms0][0], 0.1) );
647  assert( Athena_test::isEqual (vvv.Y(), exp_moms0[imoms0][1], 0.1) );
648  assert( Athena_test::isEqual (vvv.Z(), exp_moms0[imoms0][2], 0.1) );
649  assert( Athena_test::isEqual (vvv.E(), exp_moms0[imoms0][3], 0.1) );
650  ++imoms0;
651  }
652  }
653  assert (imoms0 == nmoms0);
654 
655  assert (info1->vertices().size() == 2);
656 
657  xAOD::Vertex exp_v0;
658  exp_v0.makePrivateStore();
659  exp_v0.setPosition({ 7.89827, 0.0514449, -4.04121 });
660  exp_v0.setFitQuality(5.34877, 8);
661  exp_v0.setCovariance(std::vector<float>{
662  0.218298, -0.00769266, 0.0194589, -0.0118989, 0.0107223, 0.208922 });
663  setRefittedPerigee(
664  exp_v0,
665  0,
666  1,
667  exp_mom(exp_moms0, 0),
668  { 0.209404, -4.58753, -0.00519457, 0.00377293, 0.00105397,
669  -4.58753, 154483, 32.603, -0.117209, -35.6007,
670  -0.00519457, 32.603, 0.0113951, -0.000116443, -0.00751125,
671  0.00377293, -0.117209, -0.000116443, 0.009643, 2.37152e-05,
672  0.00105397, -35.6007, -0.00751125, 2.37152e-05, 0.0082042 });
673  setFitQuality(exp_v0, 0, 0.926, 2);
674  setRefittedPerigee(
675  exp_v0,
676  1,
677  -1,
678  exp_mom(exp_moms0, 1),
679  { 0.185428, 0.807536, -0.00324979, -0.000237823, 1.99968e-05,
680  0.807536, 154490, 6.26553, -0.238515, 0.168369,
681  -0.00324979, 6.26553, 0.00856011, -0.00036866, -2.23097e-05,
682  -0.000237823, -0.238515, -0.00036866, 0.00933191, 7.51824e-06,
683  1.99968e-05, 0.168369, -2.23097e-05, 7.51824e-06, 3.74483e-07 });
684  setFitQuality(exp_v0, 1, 4.142, 2);
685  setRefittedPerigee(
686  exp_v0,
687  2,
688  -1,
689  exp_mom(exp_moms0, 2),
690  { 0.191446, -9.49834, -0.00495436, -0.00130953, -4.72358e-05,
691  -9.49834, 154476, 5.42258, -0.231539, 0.236084,
692  -0.00495436, 5.42258, 0.00860893, -0.000425575, 3.61158e-05,
693  -0.00130953, -0.231539, -0.000425575, 0.00920703, -5.93187e-06,
694  -4.72358e-05, 0.236084, 3.61158e-05, -5.93187e-06, 4.8944e-07 });
695  setFitQuality(exp_v0, 2, 0.281, 2);
696 
697  compareVertex(*info1->vertices()[0], exp_v0, 0.1);
698 
699  xAOD::Vertex exp_v1;
700  exp_v1.makePrivateStore();
701  exp_v1.setPosition({ 5.31046, 0.22012, -3.80093 });
702  exp_v1.setFitQuality(3.34131, 8);
703  exp_v1.setCovariance(std::vector<float>{
704  0.0239352, 0.000977903, 0.00708136, 4.16975e-05, 0.00295894, 0.2431 });
705  setRefittedPerigee(
706  exp_v1,
707  0,
708  1,
709  exp_mom(exp_moms0, 3),
710  { 0.166917, -16.9685, -0.000309914, -0.000819924, 1.31935e-05,
711  -16.9685, 6.35682e+08, -139.016, 9.70413, -495.474,
712  -0.000309914, -139.016, 0.0100214, -7.71225e-06, 0.000103046,
713  -0.000819924, 9.70413, -7.71225e-06, 0.00998731, -5.42066e-06,
714  1.31935e-05, -495.474, 0.000103046, -5.42066e-06, 0.000386195 });
715  setFitQuality(exp_v1, 0, 1.986, 2);
716  setRefittedPerigee(
717  exp_v1,
718  1,
719  -1,
720  exp_mom(exp_moms0, 4),
721  { 0.20904, -15.3143, 0.00086181, -0.000463146, -2.83922e-05,
722  -15.3143, 6.35682e+08, 306.365, -3.09595, -2470.81,
723  0.00086181, 306.365, 0.0101467, -2.3178e-06, -0.00116433,
724  -0.000463146, -3.09595, -2.3178e-06, 0.00999678, 8.90989e-07,
725  -2.83922e-05, -2470.81, -0.00116433, 8.90989e-07, 0.00960596 });
726  setFitQuality(exp_v1, 1, 1.356, 2);
727 
728  compareVertex (*info1->vertices()[1], exp_v1, 0.1);
729 
730  return StatusCode::SUCCESS;
731 }
732 
733 
734 } // namespace Trk
Trk::VxSecVertexInfo::setSVOwnership
void setSVOwnership(bool Ownership)
Definition: VxSecVertexInfo.h:118
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
IDTPM::ndof
float ndof(const U &p)
Definition: TrackParametersHelper.h:132
LArSamples::FitterData::fitter
const ShapeFitter * fitter
Definition: ShapeFitter.cxx:23
xAOD::Vertex_v1::setPosition
void setPosition(const Amg::Vector3D &position)
Sets the 3-position.
GeV
#define GeV
Definition: PhysicsAnalysis/TauID/TauAnalysisTools/Root/HelperFunctions.cxx:17
IVertexCascadeFitter.h
xAOD::Vertex_v1::setFitQuality
void setFitQuality(float chiSquared, float numberDoF)
Set the 'Fit Quality' information.
Definition: Vertex_v1.cxx:150
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
Amg::compress
void compress(const AmgSymMatrix(N) &covMatrix, std::vector< float > &vec)
Definition: EventPrimitivesHelpers.h:56
max
#define max(a, b)
Definition: cfImp.cxx:41
Trk::VxTrackAtVertex
The VxTrackAtVertex is a common class for all present TrkVertexFitters The VxTrackAtVertex is designe...
Definition: VxTrackAtVertex.h:77
Trk::VertexID
int VertexID
Definition: IVertexCascadeFitter.h:23
physval_make_web_display.thresh
thresh
Definition: physval_make_web_display.py:36
Trk::TrkVKalVrtFitterTestAlg::execute
virtual StatusCode execute() override
Execute the algorithm.
Definition: TrkVKalVrtFitterTestAlg.cxx:367
Trk::ParametersT
Dummy class used to allow special convertors to be called for surfaces owned by a detector element.
Definition: EMErrorDetail.h:25
Trk::VxSecVertexInfo::vertices
const std::vector< xAOD::Vertex * > & vertices() const
Definition: VxSecVertexInfo.cxx:100
TrackParticleBase.h
EventPrimitivesHelpers.h
python.SystemOfUnits.MeV
int MeV
Definition: SystemOfUnits.py:154
PlotCalibFromCool.begin
begin
Definition: PlotCalibFromCool.py:94
plotBeamSpotVxVal.cov
cov
Definition: plotBeamSpotVxVal.py:201
TrkVKalVrtFitterTestAlg.h
Algorithm for testing TrkVKalVrtFitter.
ParticleTest.tp
tp
Definition: ParticleTest.py:25
Trk::z0
@ z0
Definition: ParamDefs.h:64
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:158
vec
std::vector< size_t > vec
Definition: CombinationsGeneratorTest.cxx:12
SG::ConstAccessor
Helper class to provide constant type-safe access to aux data.
Definition: ConstAccessor.h:55
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
Trk::VxCascadeInfo::nDoF
int nDoF() const
Definition: VxCascadeInfo.h:133
xAOD::Vertex_v1::setCovariance
void setCovariance(const std::vector< float > &value)
Sets the covariance matrix as a simple vector of values.
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition: EventPrimitives.h:50
mergePhysValFiles.end
end
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:93
Track.h
Trk::AmgSymMatrix
AmgSymMatrix(5) &GXFTrackState
Definition: GXFTrackState.h:156
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
Trk::ParametersT::associatedSurface
virtual const S & associatedSurface() const override final
Access to the Surface method.
Trk::IVertexCascadeFitter
Definition: IVertexCascadeFitter.h:30
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
ParticleGun_EoverP_Config.mom
mom
Definition: ParticleGun_EoverP_Config.py:63
lumiFormat.i
int i
Definition: lumiFormat.py:85
Trk::theta
@ theta
Definition: ParamDefs.h:66
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
vector
Definition: MultiHisto.h:13
xAOD::covMatrix
covMatrix
Definition: TrackMeasurement_v1.cxx:19
Trk::NeutralPerigee
ParametersT< NeutralParametersDim, Neutral, PerigeeSurface > NeutralPerigee
Definition: NeutralParameters.h:31
Athena_test::isEqual
bool isEqual(double x1, double x2, double thresh=1e-6)
Definition: FLOATassert.h:26
Trk::VxCascadeInfo::fitChi2
double fitChi2() const
Definition: VxCascadeInfo.h:134
Trk::TrkVKalVrtFitterTestAlg::m_fitter
ToolHandle< Trk::IVertexFitter > m_fitter
Definition: TrkVKalVrtFitterTestAlg.h:50
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:522
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
hist_file_dump.f
f
Definition: hist_file_dump.py:135
Trk::TrkVKalVrtFitterTestAlg::test1
StatusCode test1()
Definition: TrkVKalVrtFitterTestAlg.cxx:382
Trk::FitQuality
Class to represent and store fit qualities from track reconstruction in terms of and number of degre...
Definition: FitQuality.h:97
Trk::ParametersBase
Definition: ParametersBase.h:55
python.Utils.unixtools.which
def which(filename, env=os.environ)
UNIX-style which ---------------------------------------------------------—.
Definition: unixtools.py:39
Trk::TrkVKalVrtFitterTestAlg::test2
StatusCode test2()
Definition: TrkVKalVrtFitterTestAlg.cxx:408
Vertex.h
IVKalState.h
Base class for VKal state object.
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
Trk::d0
@ d0
Definition: ParamDefs.h:63
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
charge
double charge(const T &p)
Definition: AtlasPID.h:538
SG::AuxElement::makePrivateStore
void makePrivateStore()
Create a new (empty) private store for this object.
Definition: AuxElement.cxx:172
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Trk::VxCascadeInfo::getParticleMoms
const std::vector< std::vector< TLorentzVector > > & getParticleMoms() const
Definition: VxCascadeInfo.h:131
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
IDTPM::chiSquared
float chiSquared(const U &p)
Definition: TrackParametersHelper.h:126
VxCascadeInfo.h
Trk::TrkVKalVrtFitterTestAlg::test3
StatusCode test3()
Definition: TrkVKalVrtFitterTestAlg.cxx:459
python.SystemOfUnits.mm
int mm
Definition: SystemOfUnits.py:83
ReadCellNoiseFromCoolCompare.v2
v2
Definition: ReadCellNoiseFromCoolCompare.py:364
TrackParticle.h
python.PyAthena.v
v
Definition: PyAthena.py:154
a
TList * a
Definition: liststreamerinfos.cxx:10
FLOATassert.h
xAOD::Vertex_v1
Class describing a Vertex.
Definition: Vertex_v1.h:42
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:67
Trk::TrkVKalVrtFitterTestAlg::initialize
virtual StatusCode initialize() override
Standard Gaudi initialize method.
Definition: TrkVKalVrtFitterTestAlg.cxx:357
xAOD::TrackParticle_v1
Class describing a TrackParticle.
Definition: TrackParticle_v1.h:43
ConstAccessor.h
Helper class to provide constant type-safe access to aux data.
xAOD::NeutralParticle_v1
Class describing a NeutralParticle.
Definition: NeutralParticle_v1.h:40
PlotCalibFromCool.vv
vv
Definition: PlotCalibFromCool.py:716
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
python.compressB64.c
def c
Definition: compressB64.py:93
xAOD::Vertex_v1::setCovariancePosition
void setCovariancePosition(const AmgSymMatrix(3)&covariancePosition)
Sets the vertex covariance matrix.
Trk::TrkVKalVrtFitterTestAlg::test4
StatusCode test4()
Definition: TrkVKalVrtFitterTestAlg.cxx:585
Trk::phi0
@ phi0
Definition: ParamDefs.h:65
Trk::v
@ v
Definition: ParamDefs.h:78
Generated on Thu Nov 7 2024 21:30:23 for ATLAS Offline Software by   doxygen 1.8.18