Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FitterTestsCommon.hpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FitterTestsCommon.hpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2021 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 
9 #pragma once
10 
11 #include <boost/test/unit_test.hpp>
12 
29 
30 #include <iterator>
31 
32 using namespace Acts::UnitLiterals;
33 using namespace Acts::Test;
34 
42  double distanceMax = std::numeric_limits<double>::max();
43 
50  template <typename traj_t>
51  bool operator()(typename traj_t::ConstTrackStateProxy state) const {
52  // can't determine an outlier w/o a measurement or predicted parameters
53  if (not state.hasCalibrated() or not state.hasPredicted()) {
54  return false;
55  }
56  auto residuals = (state.effectiveCalibrated() -
57  state.effectiveProjector() * state.predicted())
58  .eval();
59  auto distance = residuals.norm();
60  return (distanceMax <= distance);
61  }
62 };
63 
67  double momentumMax = std::numeric_limits<double>::max();
68 
74  template <typename traj_t>
75  bool operator()(typename traj_t::ConstTrackStateProxy state) const {
76  // can't determine an outlier w/o a measurement or predicted parameters
77  auto momentum = fabs(1 / state.filtered()[Acts::eBoundQOverP]);
78  std::cout << "momentum : " << momentum << std::endl;
79  return (momentum <= momentumMax);
80  }
81 };
82 
83 // Construct a straight-line propagator.
84 auto makeStraightPropagator(std::shared_ptr<const Acts::TrackingGeometry> geo) {
86  cfg.resolvePassive = false;
87  cfg.resolveMaterial = true;
88  cfg.resolveSensitive = true;
92  stepper, std::move(navigator));
93 }
94 
95 // Construct a propagator using a constant magnetic field along z.
96 template <typename stepper_t>
98  std::shared_ptr<const Acts::TrackingGeometry> geo, double bz) {
100  cfg.resolvePassive = false;
101  cfg.resolveMaterial = true;
102  cfg.resolveSensitive = true;
104  auto field =
105  std::make_shared<Acts::ConstantBField>(Acts::Vector3(0.0, 0.0, bz));
106  stepper_t stepper(std::move(field));
108  std::move(stepper), std::move(navigator));
109 }
110 
111 // Put all this in a struct to avoid that all these objects are exposed as
112 // global objects in the header
113 struct FitterTester {
114  using Rng = std::default_random_engine;
115 
116  // Context objects
120 
121  // detector geometry
123  std::shared_ptr<const Acts::TrackingGeometry> geometry = geometryStore();
124 
126 
127  // expected number of measurements for the given detector
128  constexpr static size_t nMeasurements = 6u;
129 
130  // detector resolutions
131  MeasurementResolution resPixel = {MeasurementType::eLoc01, {25_um, 50_um}};
132  MeasurementResolution resStrip0 = {MeasurementType::eLoc0, {100_um}};
133  MeasurementResolution resStrip1 = {MeasurementType::eLoc1, {150_um}};
136  {Acts::GeometryIdentifier().setVolume(3).setLayer(2), resStrip0},
137  {Acts::GeometryIdentifier().setVolume(3).setLayer(4), resStrip1},
138  {Acts::GeometryIdentifier().setVolume(3).setLayer(6), resStrip0},
139  {Acts::GeometryIdentifier().setVolume(3).setLayer(8), resStrip1},
140  };
141 
142  // simulation propagator
145 
146  static std::vector<Acts::SourceLink> prepareSourceLinks(
147  const std::vector<TestSourceLink>& sourceLinks) {
148  std::vector<Acts::SourceLink> result;
149  std::transform(sourceLinks.begin(), sourceLinks.end(),
150  std::back_inserter(result),
151  [](const auto& sl) { return Acts::SourceLink{sl}; });
152  return result;
153  }
154 
156  // The testing functions
158 
159  template <typename fitter_t, typename fitter_options_t, typename parameters_t>
160  void test_ZeroFieldNoSurfaceForward(const fitter_t& fitter,
161  fitter_options_t options,
162  const parameters_t& start, Rng& rng,
163  const bool expected_reversed,
164  const bool expected_smoothed,
165  const bool doDiag) const {
166  auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
167  resolutions, rng);
168 
169  auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
170  BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
171 
172  // this is the default option. set anyway for consistency
173  options.referenceSurface = nullptr;
174 
176  Acts::ConstTrackAccessor<bool> smoothed{"smoothed"};
177 
178  auto doTest = [&](bool diag) {
181  if (diag) {
182  tracks.addColumn<bool>("reversed");
183  tracks.addColumn<bool>("smoothed");
184 
185  BOOST_CHECK(tracks.hasColumn("reversed"));
186  BOOST_CHECK(tracks.hasColumn("smoothed"));
187  }
188 
189  auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), start,
190  options, tracks);
191  BOOST_REQUIRE(res.ok());
192 
193  const auto track = res.value();
194  BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
195  BOOST_CHECK(!track.hasReferenceSurface());
196  BOOST_CHECK_EQUAL(track.nMeasurements(), sourceLinks.size());
197  BOOST_CHECK_EQUAL(track.nHoles(), 0u);
198 
199  if (diag) {
200  // check the output status flags
201  BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
202  BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
203  }
204  };
205 
206  if (doDiag) {
207  doTest(true);
208  } // with reversed & smoothed columns
209  doTest(false); // without the extra columns
210  }
211 
212  template <typename fitter_t, typename fitter_options_t, typename parameters_t>
213  void test_ZeroFieldWithSurfaceForward(const fitter_t& fitter,
214  fitter_options_t options,
215  const parameters_t& start, Rng& rng,
216  const bool expected_reversed,
217  const bool expected_smoothed,
218  const bool doDiag) const {
219  auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
220  resolutions, rng);
221  auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
222  BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
223 
224  // initial fitter options configured for backward filtering mode
225  // backward filtering requires a reference surface
226  options.referenceSurface = &start.referenceSurface();
227  // this is the default option. set anyway for consistency
228  options.propagatorPlainOptions.direction = Acts::Direction::Forward;
229 
232  tracks.addColumn<bool>("reversed");
233  tracks.addColumn<bool>("smoothed");
234 
235  auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), start,
236  options, tracks);
237  BOOST_REQUIRE(res.ok());
238 
239  const auto& track = res.value();
240  BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
241  BOOST_CHECK(track.hasReferenceSurface());
242  BOOST_CHECK_EQUAL(track.nMeasurements(), sourceLinks.size());
243  BOOST_CHECK_EQUAL(track.nHoles(), 0u);
244 
245  BOOST_CHECK(tracks.hasColumn("reversed"));
246  BOOST_CHECK(tracks.hasColumn("smoothed"));
247 
249  Acts::ConstTrackAccessor<bool> smoothed{"smoothed"};
250 
251  // check the output status flags
252  if (doDiag) {
253  BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
254  BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
255  }
256 
257  // count the number of `smoothed` states
258  if (expected_reversed && expected_smoothed) {
259  size_t nSmoothed = 0;
260  for (const auto ts : track.trackStatesReversed()) {
261  nSmoothed += ts.hasSmoothed();
262  }
263  BOOST_CHECK_EQUAL(nSmoothed, sourceLinks.size());
264  }
265  }
266 
267  template <typename fitter_t, typename fitter_options_t, typename parameters_t>
268  void test_ZeroFieldWithSurfaceBackward(const fitter_t& fitter,
269  fitter_options_t options,
270  const parameters_t& start, Rng& rng,
271  const bool expected_reversed,
272  const bool expected_smoothed,
273  const bool doDiag) const {
274  auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
275  resolutions, rng);
276  auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
277  BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
278 
279  // create a track near the tracker exit for outward->inward filtering
280  Acts::Vector4 posOuter = start.fourPosition(geoCtx);
281  posOuter[Acts::ePos0] = 3_m;
283  posOuter, start.direction(), start.qOverP(), start.covariance(),
285 
286  options.referenceSurface = &startOuter.referenceSurface();
287  options.propagatorPlainOptions.direction = Acts::Direction::Backward;
288 
291  tracks.addColumn<bool>("reversed");
292  tracks.addColumn<bool>("smoothed");
293 
294  auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), startOuter,
295  options, tracks);
296  BOOST_CHECK(res.ok());
297 
298  const auto& track = res.value();
299  BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
300  BOOST_CHECK(track.hasReferenceSurface());
301  BOOST_CHECK_EQUAL(track.nMeasurements(), sourceLinks.size());
302  BOOST_CHECK_EQUAL(track.nHoles(), 0u);
303 
305  Acts::ConstTrackAccessor<bool> smoothed{"smoothed"};
306  // check the output status flags
307  if (doDiag) {
308  BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
309  BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
310  }
311 
312  // count the number of `smoothed` states
313  if (expected_reversed && expected_smoothed) {
314  size_t nSmoothed = 0;
315  for (const auto ts : track.trackStatesReversed()) {
316  nSmoothed += ts.hasSmoothed();
317  }
318  BOOST_CHECK_EQUAL(nSmoothed, sourceLinks.size());
319  }
320  }
321 
322  template <typename fitter_t, typename fitter_options_t, typename parameters_t>
323  void test_ZeroFieldWithSurfaceAtExit(const fitter_t& fitter,
324  fitter_options_t options,
325  const parameters_t& start, Rng& rng,
326  const bool expected_reversed,
327  const bool expected_smoothed,
328  const bool doDiag) const {
329  auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
330  resolutions, rng);
331  auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
332  BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
333 
334  // create a boundless target surface near the tracker exit
335  Acts::Vector3 center(3._m, 0., 0.);
336  Acts::Vector3 normal(1., 0., 0.);
337  auto targetSurface =
338  Acts::Surface::makeShared<Acts::PlaneSurface>(center, normal);
339 
340  options.referenceSurface = targetSurface.get();
341 
344  tracks.addColumn<bool>("reversed");
345  tracks.addColumn<bool>("smoothed");
346 
347  auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), start,
348  options, tracks);
349  BOOST_REQUIRE(res.ok());
350 
351  const auto& track = res.value();
352  BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
353  BOOST_CHECK(track.hasReferenceSurface());
354  BOOST_CHECK_EQUAL(track.nMeasurements(), sourceLinks.size());
355  BOOST_CHECK_EQUAL(track.nHoles(), 0u);
356 
358  Acts::ConstTrackAccessor<bool> smoothed{"smoothed"};
359 
360  // check the output status flags
361  if (doDiag) {
362  BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
363  BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
364  }
365  }
366 
367  template <typename fitter_t, typename fitter_options_t, typename parameters_t>
368  void test_ZeroFieldShuffled(const fitter_t& fitter, fitter_options_t options,
369  const parameters_t& start, Rng& rng,
370  const bool expected_reversed,
371  const bool expected_smoothed,
372  const bool doDiag) const {
373  auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
374  resolutions, rng);
375  auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
376  BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
377 
378  options.referenceSurface = &start.referenceSurface();
379 
380  Acts::BoundVector parameters = Acts::BoundVector::Zero();
381 
384  tracks.addColumn<bool>("reversed");
385  tracks.addColumn<bool>("smoothed");
386 
388  Acts::ConstTrackAccessor<bool> smoothed{"smoothed"};
389 
390  // fit w/ all hits in order
391  {
392  auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), start,
393  options, tracks);
394  BOOST_REQUIRE(res.ok());
395 
396  const auto& track = res.value();
397  BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
398  BOOST_CHECK_EQUAL(track.nMeasurements(), sourceLinks.size());
399  BOOST_REQUIRE(track.hasReferenceSurface());
400  parameters = track.parameters();
401  BOOST_CHECK_EQUAL(track.nHoles(), 0u);
402 
403  // check the output status flags
404  if (doDiag) {
405  BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
406  BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
407  }
408  }
409  // fit w/ all hits in random order
410  {
411  decltype(sourceLinks) shuffledSourceLinks = sourceLinks;
412  std::shuffle(shuffledSourceLinks.begin(), shuffledSourceLinks.end(), rng);
413  auto res = fitter.fit(shuffledSourceLinks.begin(),
414  shuffledSourceLinks.end(), start, options, tracks);
415  BOOST_REQUIRE(res.ok());
416 
417  const auto& track = res.value();
418  BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
419  BOOST_REQUIRE(track.hasReferenceSurface());
420  // check consistency w/ un-shuffled measurements
421  CHECK_CLOSE_ABS(track.parameters(), parameters, 1e-5);
422  BOOST_CHECK_EQUAL(track.nMeasurements(), sourceLinks.size());
423  // check the output status flags
424  if (doDiag) {
425  BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
426  BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
427  }
428  }
429  }
430 
431  template <typename fitter_t, typename fitter_options_t, typename parameters_t>
432  void test_ZeroFieldWithHole(const fitter_t& fitter, fitter_options_t options,
433  const parameters_t& start, Rng& rng,
434  const bool expected_reversed,
435  const bool expected_smoothed,
436  const bool doDiag) const {
437  auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
438  resolutions, rng);
439  auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
440  BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
441 
444  tracks.addColumn<bool>("reversed");
445  tracks.addColumn<bool>("smoothed");
446 
448  Acts::ConstTrackAccessor<bool> smoothed{"smoothed"};
449 
450  // always keep the first and last measurement. leaving those in seems to not
451  // count the respective surfaces as holes.
452  for (size_t i = 1u; (i + 1u) < sourceLinks.size(); ++i) {
453  // remove the i-th measurement
454  auto withHole = sourceLinks;
455  withHole.erase(std::next(withHole.begin(), i));
456  BOOST_REQUIRE_EQUAL(withHole.size() + 1u, sourceLinks.size());
457  BOOST_TEST_INFO("Removed measurement " << i);
458 
459  auto res =
460  fitter.fit(withHole.begin(), withHole.end(), start, options, tracks);
461  BOOST_REQUIRE(res.ok());
462 
463  const auto& track = res.value();
464  BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
465  BOOST_REQUIRE(!track.hasReferenceSurface());
466  BOOST_CHECK_EQUAL(track.nMeasurements(), withHole.size());
467  // check the output status flags
468  if (doDiag) {
469  BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
470  BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
471  }
472  BOOST_CHECK_EQUAL(track.nHoles(), 1u);
473  }
474  BOOST_CHECK_EQUAL(tracks.size(), sourceLinks.size() - 2);
475  }
476 
477  template <typename fitter_t, typename fitter_options_t, typename parameters_t>
478  void test_ZeroFieldWithOutliers(const fitter_t& fitter,
479  fitter_options_t options,
480  const parameters_t& start, Rng& rng,
481  const bool expected_reversed,
482  const bool expected_smoothed,
483  const bool doDiag) const {
484  auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
485  resolutions, rng);
486  auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
487  auto outlierSourceLinks =
488  prepareSourceLinks(measurements.outlierSourceLinks);
489  BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
490  BOOST_REQUIRE_EQUAL(outlierSourceLinks.size(), nMeasurements);
491 
494  tracks.addColumn<bool>("reversed");
495  tracks.addColumn<bool>("smoothed");
496 
498  Acts::ConstTrackAccessor<bool> smoothed{"smoothed"};
499 
500  for (size_t i = 0; i < sourceLinks.size(); ++i) {
501  // replace the i-th measurement with an outlier
502  auto withOutlier = sourceLinks;
503  withOutlier[i] = outlierSourceLinks[i];
504  BOOST_REQUIRE_EQUAL(withOutlier.size(), sourceLinks.size());
505  BOOST_TEST_INFO("Replaced measurement " << i << " with outlier");
506 
507  auto res = fitter.fit(withOutlier.begin(), withOutlier.end(), start,
508  options, tracks);
509  BOOST_REQUIRE(res.ok());
510 
511  const auto& track = res.value();
512  BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
513  // count the number of outliers
514  size_t nOutliers = 0;
515  for (const auto state : track.trackStatesReversed()) {
516  nOutliers += state.typeFlags().test(Acts::TrackStateFlag::OutlierFlag);
517  }
518  BOOST_CHECK_EQUAL(nOutliers, 1u);
519  BOOST_REQUIRE(!track.hasReferenceSurface());
520  BOOST_CHECK_EQUAL(track.nMeasurements(), withOutlier.size() - 1u);
521  // check the output status flags
522  if (doDiag) {
523  BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
524  BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
525  }
526  BOOST_CHECK_EQUAL(track.nHoles(), 0u);
527  }
528  BOOST_CHECK_EQUAL(tracks.size(), sourceLinks.size());
529  }
530 
531  template <typename fitter_t, typename fitter_options_t, typename parameters_t>
532  void test_ZeroFieldWithReverseFiltering(const fitter_t& fitter,
533  fitter_options_t options,
534  const parameters_t& start, Rng& rng,
535  const bool expected_reversed,
536  const bool expected_smoothed,
537  const bool doDiag) const {
538  auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
539  resolutions, rng);
540 
543  tracks.addColumn<bool>("reversed");
544  tracks.addColumn<bool>("smoothed");
545 
547  Acts::ConstTrackAccessor<bool> smoothed{"smoothed"};
548 
549  auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
550 
551  const auto& outlierSourceLinks = measurements.outlierSourceLinks;
552  BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
553  BOOST_REQUIRE_EQUAL(outlierSourceLinks.size(), nMeasurements);
554 
555  // create a boundless target surface near the tracker entry
556  Acts::Vector3 center(-3._m, 0., 0.);
557  Acts::Vector3 normal(1., 0., 0.);
558  auto targetSurface =
559  Acts::Surface::makeShared<Acts::PlaneSurface>(center, normal);
560 
561  options.referenceSurface = targetSurface.get();
562 
563  auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), start,
564  options, tracks);
565  BOOST_REQUIRE(res.ok());
566  const auto& track = res.value();
567 
568  // Track of 1 GeV with a threshold set at 0.1 GeV, reversed filtering should
569  // not be used
570  if (doDiag) {
571  BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
572  BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
573  }
574  }
575 
576  // TODO this is not really Kalman fitter specific. is probably better tested
577  // with a synthetic trajectory.
578  template <typename fitter_t, typename fitter_options_t, typename parameters_t>
579  void test_GlobalCovariance(const fitter_t& fitter, fitter_options_t options,
580  const parameters_t& start, Rng& rng) const {
581  auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
582  resolutions, rng);
583  auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
584  BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
585 
588 
589  auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), start,
590  options, tracks);
591  BOOST_REQUIRE(res.ok());
592 
593  // Calculate global track parameters covariance matrix
594  const auto& track = res.value();
595  auto [trackParamsCov, stateRowIndices] =
597  tracks.trackStateContainer(), track.tipIndex());
598  BOOST_CHECK_EQUAL(trackParamsCov.rows(),
599  sourceLinks.size() * Acts::eBoundSize);
600  BOOST_CHECK_EQUAL(stateRowIndices.size(), sourceLinks.size());
601  // Each smoothed track state will have eBoundSize rows/cols in the global
602  // covariance. stateRowIndices is a map of the starting row/index with the
603  // state tip as the key. Thus, the last track state (i.e. the state
604  // corresponding track.tipIndex()) has a starting row/index =
605  // eBoundSize * (nMeasurements - 1), i.e. 6*(6-1) = 30.
606  BOOST_CHECK_EQUAL(stateRowIndices.at(track.tipIndex()),
607  Acts::eBoundSize * (nMeasurements - 1));
608  }
609 };