ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
StepperTests.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file StepperTests.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2018-2020 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 #include <boost/test/unit_test.hpp>
10 
11 #include <fstream>
12 
35 
36 namespace tt = boost::test_tools;
37 using namespace Acts::UnitLiterals;
38 
39 namespace Acts {
40 namespace Test {
41 
43 
44 // Create a test context
47 
49 struct PropState {
51  PropState(EigenStepper<ConstantBField>::State sState) : stepping(sState) {}
55  struct {
56  double mass = 42.;
57  double tolerance = 1e-4;
58  double stepSizeCutOff = 0.;
59  unsigned int maxRungeKuttaStepTrials = 10000;
60  } options;
61 };
62 
66 struct EndOfWorld {
68  double maxX = 1_m;
69 
71  EndOfWorld() = default;
72 
80  template <typename propagator_state_t, typename stepper_t>
81  bool operator()(propagator_state_t& state, const stepper_t& stepper) const {
82  const double tolerance = state.options.targetTolerance;
83  if (maxX - std::abs(stepper.position(state.stepping).x()) <= tolerance ||
84  std::abs(stepper.position(state.stepping).y()) >= 0.5_m ||
85  std::abs(stepper.position(state.stepping).z()) >= 0.5_m)
86  return true;
87  return false;
88  }
89 };
90 
94 struct StepCollector {
98  struct this_result {
99  // Position of the propagator after each step
100  std::vector<Vector3D> position;
101  // Momentum of the propagator after each step
102  std::vector<Vector3D> momentum;
103  };
104 
106 
115  template <typename propagator_state_t, typename stepper_t>
116  void operator()(propagator_state_t& state, const stepper_t& stepper,
117  result_type& result) const {
118  result.position.push_back(stepper.position(state.stepping));
119  result.momentum.push_back(stepper.momentum(state.stepping) *
120  stepper.direction(state.stepping));
121  }
122 };
123 
125 BOOST_AUTO_TEST_CASE(eigen_stepper_state_test) {
126  // Set up some variables
128  double stepSize = 123.;
129  double tolerance = 234.;
130  ConstantBField bField(Vector3D(1., 2.5, 33.33));
131 
132  Vector3D pos(1., 2., 3.);
133  Vector3D mom(4., 5., 6.);
134  double time = 7.;
135  double charge = -1.;
136 
137  // Test charged parameters without covariance matrix
138  CurvilinearParameters cp(std::nullopt, pos, mom, charge, time);
140  stepSize, tolerance);
141 
142  // Test the result & compare with the input/test for reasonable members
143  BOOST_TEST(esState.jacToGlobal == BoundToFreeMatrix::Zero());
144  BOOST_TEST(esState.jacTransport == FreeMatrix::Identity());
145  BOOST_TEST(esState.derivative == FreeVector::Zero());
146  BOOST_TEST(!esState.covTransport);
147  BOOST_TEST(esState.cov == Covariance::Zero());
148  BOOST_TEST(esState.pos == pos);
149  BOOST_TEST(esState.dir == mom.normalized());
150  BOOST_TEST(esState.p == mom.norm());
151  BOOST_TEST(esState.q == charge);
152  BOOST_TEST(esState.t == time);
153  BOOST_TEST(esState.navDir == ndir);
154  BOOST_TEST(esState.pathAccumulated == 0.);
155  BOOST_TEST(esState.stepSize == ndir * stepSize);
156  BOOST_TEST(esState.previousStepSize == 0.);
157  BOOST_TEST(esState.tolerance == tolerance);
158 
159  // Test without charge and covariance matrix
160  NeutralCurvilinearParameters ncp(std::nullopt, pos, mom, time);
162  stepSize, tolerance);
163  BOOST_TEST(esState.q == 0.);
164 
165  // Test with covariance matrix
166  Covariance cov = 8. * Covariance::Identity();
167  ncp = NeutralCurvilinearParameters(cov, pos, mom, time);
169  stepSize, tolerance);
170  BOOST_TEST(esState.jacToGlobal != BoundToFreeMatrix::Zero());
171  BOOST_TEST(esState.covTransport);
172  BOOST_TEST(esState.cov == cov);
173 }
174 
177 BOOST_AUTO_TEST_CASE(eigen_stepper_test) {
178  // Set up some variables for the state
180  double stepSize = 123.;
181  double tolerance = 234.;
182  ConstantBField bField(Vector3D(1., 2.5, 33.33));
183 
184  // Construct the parameters
185  Vector3D pos(1., 2., 3.);
186  Vector3D mom(4., 5., 6.);
187  double time = 7.;
188  double charge = -1.;
189  Covariance cov = 8. * Covariance::Identity();
190  CurvilinearParameters cp(cov, pos, mom, charge, time);
191 
192  // Build the state and the stepper
194  stepSize, tolerance);
195  EigenStepper<ConstantBField> es(bField);
196 
197  // Test the getters
198  BOOST_TEST(es.position(esState) == esState.pos);
199  BOOST_TEST(es.direction(esState) == esState.dir);
200  BOOST_TEST(es.momentum(esState) == esState.p);
201  BOOST_TEST(es.charge(esState) == esState.q);
202  BOOST_TEST(es.time(esState) == esState.t);
203  //~ BOOST_TEST(es.overstepLimit(esState) == tolerance);
204  BOOST_TEST(es.getField(esState, pos) == bField.getField(pos));
205 
206  // Step size modifies
207  const std::string originalStepSize = esState.stepSize.toString();
208 
209  es.setStepSize(esState, 1337.);
210  BOOST_TEST(esState.previousStepSize == ndir * stepSize);
211  BOOST_TEST(esState.stepSize == 1337.);
212 
213  es.releaseStepSize(esState);
214  BOOST_TEST(esState.stepSize == -123.);
215  BOOST_TEST(es.outputStepSize(esState) == originalStepSize);
216 
217  // Test the curvilinear state construction
218  auto curvState = es.curvilinearState(esState);
219  auto curvPars = std::get<0>(curvState);
220  CHECK_CLOSE_ABS(curvPars.position(), cp.position(), 1e-6);
221  CHECK_CLOSE_ABS(curvPars.momentum(), cp.momentum(), 1e-6);
222  CHECK_CLOSE_ABS(curvPars.charge(), cp.charge(), 1e-6);
223  CHECK_CLOSE_ABS(curvPars.time(), cp.time(), 1e-6);
224  BOOST_TEST(curvPars.covariance().has_value());
225  BOOST_TEST(*curvPars.covariance() != cov);
226  CHECK_CLOSE_COVARIANCE(std::get<1>(curvState),
227  BoundMatrix(BoundMatrix::Identity()), 1e-6);
228  CHECK_CLOSE_ABS(std::get<2>(curvState), 0., 1e-6);
229 
230  // Test the update method
231  Vector3D newPos(2., 4., 8.);
232  Vector3D newMom(3., 9., 27.);
233  double newTime(321.);
234  es.update(esState, newPos, newMom.normalized(), newMom.norm(), newTime);
235  BOOST_TEST(esState.pos == newPos);
236  BOOST_TEST(esState.dir == newMom.normalized());
237  BOOST_TEST(esState.p == newMom.norm());
238  BOOST_TEST(esState.q == charge);
239  BOOST_TEST(esState.t == newTime);
240 
241  // The covariance transport
242  esState.cov = cov;
243  es.covarianceTransport(esState);
244  BOOST_TEST(esState.cov != cov);
245  BOOST_TEST(esState.jacToGlobal != BoundToFreeMatrix::Zero());
246  BOOST_TEST(esState.jacTransport == FreeMatrix::Identity());
247  BOOST_TEST(esState.derivative == FreeVector::Zero());
248 
249  // Perform a step without and with covariance transport
250  esState.cov = cov;
251  PropState ps(esState);
252 
253  ps.stepping.covTransport = false;
254  double h = es.step(ps).value();
255  BOOST_TEST(ps.stepping.stepSize == h);
256  CHECK_CLOSE_COVARIANCE(ps.stepping.cov, cov, 1e-6);
257  BOOST_TEST(ps.stepping.pos.norm() != newPos.norm());
258  BOOST_TEST(ps.stepping.dir != newMom.normalized());
259  BOOST_TEST(ps.stepping.q == charge);
260  BOOST_TEST(ps.stepping.t < newTime);
261  BOOST_TEST(ps.stepping.derivative == FreeVector::Zero());
262  BOOST_TEST(ps.stepping.jacTransport == FreeMatrix::Identity());
263 
264  ps.stepping.covTransport = true;
265  double h2 = es.step(ps).value();
266  BOOST_TEST(h2 == h);
267  CHECK_CLOSE_COVARIANCE(ps.stepping.cov, cov, 1e-6);
268  BOOST_TEST(ps.stepping.pos.norm() != newPos.norm());
269  BOOST_TEST(ps.stepping.dir != newMom.normalized());
270  BOOST_TEST(ps.stepping.q == charge);
271  BOOST_TEST(ps.stepping.t < newTime);
272  BOOST_TEST(ps.stepping.derivative != FreeVector::Zero());
273  BOOST_TEST(ps.stepping.jacTransport != FreeMatrix::Identity());
274 
276  auto plane = Surface::makeShared<PlaneSurface>(pos, mom.normalized());
277  BoundParameters bp(tgContext, cov, pos, mom, charge, time, plane);
279  stepSize, tolerance);
280 
281  // Test the intersection in the context of a surface
282  auto targetSurface = Surface::makeShared<PlaneSurface>(
283  pos + ndir * 2. * mom.normalized(), mom.normalized());
284  es.updateSurfaceStatus(esState, *targetSurface, BoundaryCheck(false));
285  BOOST_TEST(esState.stepSize.value(ConstrainedStep::actor), ndir * 2.);
286 
287  // Test the step size modification in the context of a surface
288  es.updateStepSize(
289  esState,
290  targetSurface->intersect(esState.geoContext, esState.pos,
291  esState.navDir * esState.dir, false),
292  false);
293  BOOST_TEST(esState.stepSize == 2.);
294  esState.stepSize = ndir * stepSize;
295  es.updateStepSize(
296  esState,
297  targetSurface->intersect(esState.geoContext, esState.pos,
298  esState.navDir * esState.dir, false),
299  true);
300  BOOST_TEST(esState.stepSize == 2.);
301 
302  // Test the bound state construction
303  auto boundState = es.boundState(esState, *plane);
304  auto boundPars = std::get<0>(boundState);
305  CHECK_CLOSE_ABS(boundPars.position(), bp.position(), 1e-6);
306  CHECK_CLOSE_ABS(boundPars.momentum(), bp.momentum(), 1e-6);
307  CHECK_CLOSE_ABS(boundPars.charge(), bp.charge(), 1e-6);
308  CHECK_CLOSE_ABS(boundPars.time(), bp.time(), 1e-6);
309  BOOST_TEST(boundPars.covariance().has_value());
310  BOOST_TEST(*boundPars.covariance() != cov);
311  CHECK_CLOSE_COVARIANCE(std::get<1>(boundState),
312  BoundMatrix(BoundMatrix::Identity()), 1e-6);
313  CHECK_CLOSE_ABS(std::get<2>(boundState), 0., 1e-6);
314 
315  // Update in context of a surface
316  BoundParameters bpTarget(tgContext, 2. * cov, 2. * pos, 2. * mom,
317  -1. * charge, 2. * time, targetSurface);
318  es.update(esState, bpTarget);
319  BOOST_TEST(esState.pos == 2. * pos);
320  BOOST_TEST(esState.dir == mom.normalized());
321  BOOST_TEST(esState.p == 2. * mom.norm());
322  BOOST_TEST(esState.q == 1. * charge);
323  BOOST_TEST(esState.t == 2. * time);
324  CHECK_CLOSE_COVARIANCE(esState.cov, Covariance(2. * cov), 1e-6);
325 
326  // Transport the covariance in the context of a surface
327  es.covarianceTransport(esState, *plane);
328  BOOST_TEST(esState.cov != cov);
329  BOOST_TEST(esState.jacToGlobal != BoundToFreeMatrix::Zero());
330  BOOST_TEST(esState.jacTransport == FreeMatrix::Identity());
331  BOOST_TEST(esState.derivative == FreeVector::Zero());
332 }
333 
342 
343 // Test case a). The DenseEnvironmentExtension should state that it is not
344 // valid in this case.
345 BOOST_AUTO_TEST_CASE(step_extension_vacuum_test) {
348  vConf.position = {0.5_m, 0., 0.};
349  vConf.length = {1_m, 1_m, 1_m};
351  conf.volumeCfg.push_back(vConf);
352  conf.position = {0.5_m, 0., 0.};
353  conf.length = {1_m, 1_m, 1_m};
354 
355  // Build detector
356  cvb.setConfig(conf);
358  tgbCfg.trackingVolumeBuilders.push_back(
359  [=](const auto& context, const auto& inner, const auto& vb) {
360  return cvb.trackingVolume(context, inner, vb);
361  });
362  TrackingGeometryBuilder tgb(tgbCfg);
363  std::shared_ptr<const TrackingGeometry> vacuum =
365 
366  // Build navigator
367  Navigator naviVac(vacuum);
368  naviVac.resolvePassive = true;
369  naviVac.resolveMaterial = true;
370  naviVac.resolveSensitive = true;
371 
372  // Set initial parameters for the particle track
373  Covariance cov = Covariance::Identity();
374  Vector3D startParams(0., 0., 0.), startMom(1_GeV, 0., 0.);
376  startMom, 1., 0.);
377 
378  // Create action list for surface collection
380  AbortList<EndOfWorld> abortList;
381 
382  // Set options for propagator
385  propOpts(tgContext, mfContext);
386  propOpts.actionList = aList;
387  propOpts.abortList = abortList;
388  propOpts.maxSteps = 100;
389  propOpts.maxStepSize = 1.5_m;
390 
391  // Build stepper and propagator
392  ConstantBField bField(Vector3D(0., 0., 0.));
393  EigenStepper<
397  es(bField);
398  Propagator<EigenStepper<ConstantBField,
402  Navigator>
403  prop(es, naviVac);
404 
405  // Launch and collect results
406  const auto& result = prop.propagate(sbtp, propOpts).value();
407  const StepCollector::this_result& stepResult =
408  result.get<typename StepCollector::result_type>();
409 
410  // Check that the propagation happend without interactions
411  for (const auto& pos : stepResult.position) {
412  CHECK_SMALL(pos.y(), 1_um);
413  CHECK_SMALL(pos.z(), 1_um);
414  if (pos == stepResult.position.back())
415  CHECK_CLOSE_ABS(pos.x(), 1_m, 1_um);
416  }
417  for (const auto& mom : stepResult.momentum) {
418  CHECK_CLOSE_ABS(mom, startMom, 1_keV);
419  }
420 
421  // Rebuild and check the choice of extension
422  ActionList<StepCollector> aListDef;
423 
424  // Set options for propagator
426  propOptsDef(tgContext, mfContext);
427  propOptsDef.actionList = aListDef;
428  propOptsDef.abortList = abortList;
429  propOptsDef.maxSteps = 100;
430  propOptsDef.maxStepSize = 1.5_m;
431 
433  bField);
434  Propagator<
436  Navigator>
437  propDef(esDef, naviVac);
438 
439  // Launch and collect results
440  const auto& resultDef = propDef.propagate(sbtp, propOptsDef).value();
441  const StepCollector::this_result& stepResultDef =
442  resultDef.get<typename StepCollector::result_type>();
443 
444  // Check that the right extension was chosen
445  // If chosen correctly, the number of elements should be identical
446  BOOST_TEST(stepResult.position.size() == stepResultDef.position.size());
447  for (unsigned int i = 0; i < stepResult.position.size(); i++) {
448  CHECK_CLOSE_ABS(stepResult.position[i], stepResultDef.position[i], 1_um);
449  }
450  BOOST_TEST(stepResult.momentum.size() == stepResultDef.momentum.size());
451  for (unsigned int i = 0; i < stepResult.momentum.size(); i++) {
452  CHECK_CLOSE_ABS(stepResult.momentum[i], stepResultDef.momentum[i], 1_keV);
453  }
454 }
455 // Test case b). The DefaultExtension should state that it is invalid here.
456 BOOST_AUTO_TEST_CASE(step_extension_material_test) {
459  vConf.position = {0.5_m, 0., 0.};
460  vConf.length = {1_m, 1_m, 1_m};
461  vConf.volumeMaterial = std::make_shared<const HomogeneousVolumeMaterial>(
462  Material(352.8, 394.133, 9.012, 4., 1.848e-3));
464  conf.volumeCfg.push_back(vConf);
465  conf.position = {0.5_m, 0., 0.};
466  conf.length = {1_m, 1_m, 1_m};
467 
468  // Build detector
469  cvb.setConfig(conf);
471  tgbCfg.trackingVolumeBuilders.push_back(
472  [=](const auto& context, const auto& inner, const auto& vb) {
473  return cvb.trackingVolume(context, inner, vb);
474  });
475  TrackingGeometryBuilder tgb(tgbCfg);
476  std::shared_ptr<const TrackingGeometry> material =
478 
479  // Build navigator
480  Navigator naviMat(material);
481  naviMat.resolvePassive = true;
482  naviMat.resolveMaterial = true;
483  naviMat.resolveSensitive = true;
484 
485  // Set initial parameters for the particle track
486  Covariance cov = Covariance::Identity();
487  Vector3D startParams(0., 0., 0.), startMom(5_GeV, 0., 0.);
489  startMom, 1., 0.);
490 
491  // Create action list for surface collection
493  AbortList<EndOfWorld> abortList;
494 
495  // Set options for propagator
498  propOpts(tgContext, mfContext);
499  propOpts.actionList = aList;
500  propOpts.abortList = abortList;
501  propOpts.maxSteps = 100;
502  propOpts.maxStepSize = 1.5_m;
503  propOpts.debug = true;
504 
505  // Build stepper and propagator
506  ConstantBField bField(Vector3D(0., 0., 0.));
507  EigenStepper<
511  es(bField);
512  Propagator<EigenStepper<ConstantBField,
516  Navigator>
517  prop(es, naviMat);
518 
519  // Launch and collect results
520  const auto& result = prop.propagate(sbtp, propOpts).value();
521  const StepCollector::this_result& stepResult =
522  result.get<typename StepCollector::result_type>();
523 
524  // Check that there occured interaction
525  for (const auto& pos : stepResult.position) {
526  CHECK_SMALL(pos.y(), 1_um);
527  CHECK_SMALL(pos.z(), 1_um);
528  if (pos == stepResult.position.front()) {
529  CHECK_SMALL(pos.x(), 1_um);
530  } else {
531  BOOST_CHECK_GT(std::abs(pos.x()), 1_um);
532  }
533  }
534  for (const auto& mom : stepResult.momentum) {
535  CHECK_SMALL(mom.y(), 1_keV);
536  CHECK_SMALL(mom.z(), 1_keV);
537  if (mom == stepResult.momentum.front()) {
538  CHECK_CLOSE_ABS(mom.x(), 5_GeV, 1_keV);
539  } else {
540  BOOST_CHECK_LT(mom.x(), 5_GeV);
541  }
542  }
543 
544  // Rebuild and check the choice of extension
545  // Set options for propagator
548  propOptsDense(tgContext, mfContext);
549  propOptsDense.actionList = aList;
550  propOptsDense.abortList = abortList;
551  propOptsDense.maxSteps = 100;
552  propOptsDense.maxStepSize = 1.5_m;
553  propOptsDense.debug = true;
554 
555  // Build stepper and propagator
557  esDense(bField);
558  Propagator<EigenStepper<ConstantBField,
560  Navigator>
561  propDense(esDense, naviMat);
562 
563  // Launch and collect results
564  const auto& resultDense = propDense.propagate(sbtp, propOptsDense).value();
565  const StepCollector::this_result& stepResultDense =
566  resultDense.get<typename StepCollector::result_type>();
567 
568  // Check that the right extension was chosen
569  // If chosen correctly, the number of elements should be identical
570  BOOST_TEST(stepResult.position.size() == stepResultDense.position.size());
571  for (unsigned int i = 0; i < stepResult.position.size(); i++) {
572  CHECK_CLOSE_ABS(stepResult.position[i], stepResultDense.position[i], 1_um);
573  }
574  BOOST_TEST(stepResult.momentum.size() == stepResultDense.momentum.size());
575  for (unsigned int i = 0; i < stepResult.momentum.size(); i++) {
576  CHECK_CLOSE_ABS(stepResult.momentum[i], stepResultDense.momentum[i], 1_keV);
577  }
578 
580 
581  // Re-launch the configuration with magnetic field
582  bField.setField(0., 1_T, 0.);
583  EigenStepper<
584  ConstantBField,
585  StepperExtensionList<DefaultExtension, DenseEnvironmentExtension>,
586  detail::HighestValidAuctioneer>
587  esB(bField);
588  Propagator<EigenStepper<ConstantBField,
589  StepperExtensionList<DefaultExtension,
590  DenseEnvironmentExtension>,
591  detail::HighestValidAuctioneer>,
592  Navigator>
593  propB(esB, naviMat);
594 
595  const auto& resultB = propB.propagate(sbtp, propOptsDense).value();
596  const StepCollector::this_result& stepResultB =
597  resultB.get<typename StepCollector::result_type>();
598 
599  // Check that there occured interaction
600  for (const auto& pos : stepResultB.position) {
601  if (pos == stepResultB.position.front()) {
602  CHECK_SMALL(pos, 1_um);
603  } else {
604  BOOST_CHECK_GT(std::abs(pos.x()), 1_um);
605  CHECK_SMALL(pos.y(), 1_um);
606  BOOST_CHECK_GT(std::abs(pos.z()), 1_um);
607  }
608  }
609  for (const auto& mom : stepResultB.momentum) {
610  if (mom == stepResultB.momentum.front()) {
611  CHECK_CLOSE_ABS(mom, startMom, 1_keV);
612  } else {
613  BOOST_CHECK_NE(mom.x(), 5_GeV);
614  CHECK_SMALL(mom.y(), 1_keV);
615  BOOST_CHECK_NE(mom.z(), 0.);
616  }
617  }
618 }
619 // Test case c). Both should be involved in their part of the detector
620 BOOST_AUTO_TEST_CASE(step_extension_vacmatvac_test) {
623  vConfVac1.position = {0.5_m, 0., 0.};
624  vConfVac1.length = {1_m, 1_m, 1_m};
625  vConfVac1.name = "First vacuum volume";
627  vConfMat.position = {1.5_m, 0., 0.};
628  vConfMat.length = {1_m, 1_m, 1_m};
629  vConfMat.volumeMaterial = std::make_shared<const HomogeneousVolumeMaterial>(
630  Material(352.8, 394.133, 9.012, 4., 1.848e-3));
631  vConfMat.name = "Material volume";
633  vConfVac2.position = {2.5_m, 0., 0.};
634  vConfVac2.length = {1_m, 1_m, 1_m};
635  vConfVac2.name = "Second vacuum volume";
637  conf.volumeCfg = {vConfVac1, vConfMat, vConfVac2};
638  conf.position = {1.5_m, 0., 0.};
639  conf.length = {3_m, 1_m, 1_m};
640 
641  // Build detector
642  cvb.setConfig(conf);
644  tgbCfg.trackingVolumeBuilders.push_back(
645  [=](const auto& context, const auto& inner, const auto& vb) {
646  return cvb.trackingVolume(context, inner, vb);
647  });
648  TrackingGeometryBuilder tgb(tgbCfg);
649  std::shared_ptr<const TrackingGeometry> det = tgb.trackingGeometry(tgContext);
650 
651  // Build navigator
652  Navigator naviDet(det);
653  naviDet.resolvePassive = true;
654  naviDet.resolveMaterial = true;
655  naviDet.resolveSensitive = true;
656 
657  // Set initial parameters for the particle track
658  Covariance cov = Covariance::Identity();
659  Vector3D startParams(0., 0., 0.), startMom(5_GeV, 0., 0.);
661  startMom, 1., 0.);
662 
663  // Create action list for surface collection
664  AbortList<EndOfWorld> abortList;
665  abortList.get<EndOfWorld>().maxX = 3_m;
666 
668 
669  // Set options for propagator
672  propOpts(tgContext, mfContext);
673  propOpts.abortList = abortList;
674  propOpts.maxSteps = 100;
675  propOpts.maxStepSize = 1.5_m;
676 
677  // Build stepper and propagator
678  ConstantBField bField(Vector3D(0., 1_T, 0.));
679  EigenStepper<
683  es(bField);
684  Propagator<EigenStepper<ConstantBField,
688  Navigator>
689  prop(es, naviDet);
690 
691  // Launch and collect results
692  const auto& result = prop.propagate(sbtp, propOpts).value();
693  const StepCollector::this_result& stepResult =
694  result.get<typename StepCollector::result_type>();
695 
696  // Manually set the extensions for each step and propagate through each
697  // volume by propagation to the boundaries
698  // Collect boundaries
699  std::vector<Surface const*> surs;
700  std::vector<std::shared_ptr<const BoundarySurfaceT<TrackingVolume>>>
701  boundaries = det->lowestTrackingVolume(tgContext, {0.5_m, 0., 0.})
702  ->boundarySurfaces();
703  for (auto& b : boundaries) {
704  if (b->surfaceRepresentation().center(tgContext).x() == 1_m) {
705  surs.push_back(&(b->surfaceRepresentation()));
706  break;
707  }
708  }
709  boundaries =
710  det->lowestTrackingVolume(tgContext, {1.5_m, 0., 0.})->boundarySurfaces();
711  for (auto& b : boundaries) {
712  if (b->surfaceRepresentation().center(tgContext).x() == 2_m) {
713  surs.push_back(&(b->surfaceRepresentation()));
714  break;
715  }
716  }
717  boundaries =
718  det->lowestTrackingVolume(tgContext, {2.5_m, 0., 0.})->boundarySurfaces();
719  for (auto& b : boundaries) {
720  if (b->surfaceRepresentation().center(tgContext).x() == 3_m) {
721  surs.push_back(&(b->surfaceRepresentation()));
722  break;
723  }
724  }
725 
726  // Build launcher through vacuum
727  // Set options for propagator
728 
731  propOptsDef(tgContext, mfContext);
732  abortList.get<EndOfWorld>().maxX = 1_m;
733  propOptsDef.abortList = abortList;
734  propOptsDef.maxSteps = 100;
735  propOptsDef.maxStepSize = 1.5_m;
736  propOptsDef.debug = false;
737 
738  // Build stepper and propagator
740  bField);
741  Propagator<
743  Navigator>
744  propDef(esDef, naviDet);
745 
746  // Launch and collect results
747  const auto& resultDef =
748  propDef.propagate(sbtp, *(surs[0]), propOptsDef).value();
749  const StepCollector::this_result& stepResultDef =
750  resultDef.get<typename StepCollector::result_type>();
751 
752  // Check the exit situation of the first volume
753  std::pair<Vector3D, Vector3D> endParams, endParamsControl;
754  for (unsigned int i = 0; i < stepResultDef.position.size(); i++) {
755  if (1_m - stepResultDef.position[i].x() < 1e-4) {
756  endParams =
757  std::make_pair(stepResultDef.position[i], stepResultDef.momentum[i]);
758  break;
759  }
760  }
761  for (unsigned int i = 0; i < stepResult.position.size(); i++) {
762  if (1_m - stepResult.position[i].x() < 1e-4) {
763  endParamsControl =
764  std::make_pair(stepResult.position[i], stepResult.momentum[i]);
765  break;
766  }
767  }
768 
769  if (propOptsDef.debug) {
770  const auto debugString =
771  resultDef.template get<DebugOutput::result_type>().debugString;
772  std::cout << debugString << std::endl;
773  }
774 
775  CHECK_CLOSE_ABS(endParams.first, endParamsControl.first, 1_um);
776  CHECK_CLOSE_ABS(endParams.second, endParamsControl.second, 1_um);
777 
778  BOOST_TEST(endParams.first.x() == endParamsControl.first.x(),
779  tt::tolerance(1e-5));
780  BOOST_TEST(endParams.first.y() == endParamsControl.first.y(),
781  tt::tolerance(1e-5));
782  BOOST_TEST(endParams.first.z() == endParamsControl.first.z(),
783  tt::tolerance(1e-5));
784  BOOST_TEST(endParams.second.x() == endParamsControl.second.x(),
785  tt::tolerance(1e-5));
786  BOOST_TEST(endParams.second.y() == endParamsControl.second.y(),
787  tt::tolerance(1e-5));
788  BOOST_TEST(endParams.second.z() == endParamsControl.second.z(),
789  tt::tolerance(1e-5));
790 
791  // Build launcher through material
792  // Set initial parameters for the particle track by using the result of the
793  // first volume
794  startParams = endParams.first;
795  startMom = endParams.second;
797  cov, startParams, startMom, 1., 0.);
798 
799  // Set options for propagator
802  propOptsDense(tgContext, mfContext);
803  abortList.get<EndOfWorld>().maxX = 2_m;
804  propOptsDense.abortList = abortList;
805  propOptsDense.maxSteps = 1000;
806  propOptsDense.maxStepSize = 1.5_m;
807  propOptsDense.tolerance = 1e-8;
808 
809  // Build stepper and propagator
811  esDense(bField);
812  Propagator<EigenStepper<ConstantBField,
814  Navigator>
815  propDense(esDense, naviDet);
816 
817  // Launch and collect results
818  const auto& resultDense =
819  propDense.propagate(sbtpPiecewise, *(surs[1]), propOptsDense).value();
820  const StepCollector::this_result& stepResultDense =
821  resultDense.get<typename StepCollector::result_type>();
822 
823  // Check the exit situation of the second volume
824  for (unsigned int i = 0; i < stepResultDense.position.size(); i++) {
825  if (2_m - stepResultDense.position[i].x() < 1e-4) {
826  endParams = std::make_pair(stepResultDense.position[i],
827  stepResultDense.momentum[i]);
828  break;
829  }
830  }
831  for (unsigned int i = 0; i < stepResult.position.size(); i++) {
832  if (2_m - stepResult.position[i].x() < 1e-4) {
833  endParamsControl =
834  std::make_pair(stepResult.position[i], stepResult.momentum[i]);
835  break;
836  }
837  }
838 
839  CHECK_CLOSE_ABS(endParams.first, endParamsControl.first, 1_um);
840  CHECK_CLOSE_ABS(endParams.second, endParamsControl.second, 1_um);
841 }
842 
843 // Test case a). The DenseEnvironmentExtension should state that it is not
844 // valid in this case.
845 BOOST_AUTO_TEST_CASE(step_extension_trackercalomdt_test) {
846  double rotationAngle = M_PI * 0.5;
847  Vector3D xPos(cos(rotationAngle), 0., sin(rotationAngle));
848  Vector3D yPos(0., 1., 0.);
849  Vector3D zPos(-sin(rotationAngle), 0., cos(rotationAngle));
850  MaterialProperties matProp(352.8, 407., 9.012, 4., 1.848e-3, 0.5_mm);
851 
854  sConf1.position = Vector3D(0.3_m, 0., 0.);
855  sConf1.rotation.col(0) = xPos;
856  sConf1.rotation.col(1) = yPos;
857  sConf1.rotation.col(2) = zPos;
858  sConf1.rBounds =
859  std::make_shared<const RectangleBounds>(RectangleBounds(0.5_m, 0.5_m));
860  sConf1.surMat = std::shared_ptr<const ISurfaceMaterial>(
861  new HomogeneousSurfaceMaterial(matProp));
862  sConf1.thickness = 1._mm;
864  lConf1.surfaceCfg = sConf1;
865 
867  sConf2.position = Vector3D(0.6_m, 0., 0.);
868  sConf2.rotation.col(0) = xPos;
869  sConf2.rotation.col(1) = yPos;
870  sConf2.rotation.col(2) = zPos;
871  sConf2.rBounds =
872  std::make_shared<const RectangleBounds>(RectangleBounds(0.5_m, 0.5_m));
873  sConf2.surMat = std::shared_ptr<const ISurfaceMaterial>(
874  new HomogeneousSurfaceMaterial(matProp));
875  sConf2.thickness = 1._mm;
877  lConf2.surfaceCfg = sConf2;
878 
880  muConf1.position = {2.3_m, 0., 0.};
881  muConf1.length = {20._cm, 20._cm, 20._cm};
882  muConf1.volumeMaterial =
883  std::shared_ptr<const IVolumeMaterial>(new HomogeneousVolumeMaterial(
884  Material(352.8, 407., 9.012, 4., 1.848e-3)));
885  muConf1.name = "MDT1";
887  muConf2.position = {2.7_m, 0., 0.};
888  muConf2.length = {20._cm, 20._cm, 20._cm};
889  muConf2.volumeMaterial =
890  std::shared_ptr<const IVolumeMaterial>(new HomogeneousVolumeMaterial(
891  Material(352.8, 407., 9.012, 4., 1.848e-3)));
892  muConf2.name = "MDT2";
893 
895  vConf1.position = {0.5_m, 0., 0.};
896  vConf1.length = {1._m, 1._m, 1._m};
897  vConf1.layerCfg = {lConf1, lConf2};
898  vConf1.name = "Tracker";
900  vConf2.position = {1.5_m, 0., 0.};
901  vConf2.length = {1._m, 1._m, 1._m};
902  vConf2.volumeMaterial =
903  std::shared_ptr<const IVolumeMaterial>(new HomogeneousVolumeMaterial(
904  Material(352.8, 407., 9.012, 4., 1.848e-3)));
905  vConf2.name = "Calorimeter";
907  vConf3.position = {2.5_m, 0., 0.};
908  vConf3.length = {1._m, 1._m, 1._m};
909  vConf3.volumeCfg = {muConf1, muConf2};
910  vConf3.name = "Muon system";
912  conf.volumeCfg = {vConf1, vConf2, vConf3};
913  conf.position = {1.5_m, 0., 0.};
914  conf.length = {3._m, 1._m, 1._m};
915 
916  // Build detector
917  cvb.setConfig(conf);
919  tgbCfg.trackingVolumeBuilders.push_back(
920  [=](const auto& context, const auto& inner, const auto& vb) {
921  return cvb.trackingVolume(context, inner, vb);
922  });
923  TrackingGeometryBuilder tgb(tgbCfg);
924  std::shared_ptr<const TrackingGeometry> detector =
926 
927  // Build navigator
928  Navigator naviVac(detector);
929  naviVac.resolvePassive = true;
930  naviVac.resolveMaterial = true;
931  naviVac.resolveSensitive = true;
932 
933  // Set initial parameters for the particle track
934  Covariance cov = Covariance::Identity();
935  Vector3D startParams(0., 0., 0.), startMom(1._GeV, 0., 0.);
937  startMom, 1., 0.);
938 
939  // Set options for propagator
942  propOpts(tgContext, mfContext);
943  propOpts.abortList.get<EndOfWorld>().maxX = 3._m;
944 
945  // Build stepper and propagator
946  ConstantBField bField(Vector3D(0., 0., 0.));
947  EigenStepper<
951  es(bField);
952  Propagator<EigenStepper<ConstantBField,
956  Navigator>
957  prop(es, naviVac);
958 
959  // Launch and collect results
960  const auto& result = prop.propagate(sbtp, propOpts).value();
961  const StepCollector::this_result& stepResult =
962  result.get<typename StepCollector::result_type>();
963 
964  // Test that momentum changes only occured at the right detector parts
965  double lastMomentum = stepResult.momentum[0].x();
966  for (unsigned int i = 0; i < stepResult.position.size(); i++) {
967  // Test for changes
968  if ((stepResult.position[i].x() > 0.3_m &&
969  stepResult.position[i].x() < 0.6_m) ||
970  (stepResult.position[i].x() > 0.6_m &&
971  stepResult.position[i].x() <= 1._m) ||
972  (stepResult.position[i].x() > 1._m &&
973  stepResult.position[i].x() <= 2._m) ||
974  (stepResult.position[i].x() > 2.2_m &&
975  stepResult.position[i].x() <= 2.4_m) ||
976  (stepResult.position[i].x() > 2.6_m &&
977  stepResult.position[i].x() <= 2.8_m)) {
978  BOOST_TEST(stepResult.momentum[i].x() <= lastMomentum);
979  lastMomentum = stepResult.momentum[i].x();
980  } else
981  // Test the absence of momentum loss
982  {
983  if (stepResult.position[i].x() < 0.3_m ||
984  (stepResult.position[i].x() > 2._m &&
985  stepResult.position[i].x() <= 2.2_m) ||
986  (stepResult.position[i].x() > 2.4_m &&
987  stepResult.position[i].x() <= 2.6_m) ||
988  (stepResult.position[i].x() > 2.8_m &&
989  stepResult.position[i].x() <= 3._m)) {
990  BOOST_TEST(stepResult.momentum[i].x() == lastMomentum);
991  }
992  }
993  }
994 }
995 } // namespace Test
996 } // namespace Acts