9 #include <boost/test/unit_test.hpp>
36 namespace tt = boost::test_tools;
37 using namespace Acts::UnitLiterals;
57 double tolerance = 1
e-4;
58 double stepSizeCutOff = 0.;
59 unsigned int maxRungeKuttaStepTrials = 10000;
80 template <
typename propagator_state_t,
typename stepper_t>
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)
115 template <
typename propagator_state_t,
typename stepper_t>
116 void operator()(propagator_state_t& state,
const stepper_t&
stepper,
118 result.
position.push_back(stepper.position(state.stepping));
119 result.
momentum.push_back(stepper.momentum(state.stepping) *
120 stepper.direction(state.stepping));
128 double stepSize = 123.;
129 double tolerance = 234.;
140 stepSize, tolerance);
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);
162 stepSize, tolerance);
163 BOOST_TEST(esState.q == 0.);
169 stepSize, tolerance);
170 BOOST_TEST(esState.jacToGlobal != BoundToFreeMatrix::Zero());
171 BOOST_TEST(esState.covTransport);
172 BOOST_TEST(esState.cov == cov);
180 double stepSize = 123.;
181 double tolerance = 234.;
194 stepSize, tolerance);
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);
207 const std::string originalStepSize = esState.stepSize.toString();
210 BOOST_TEST(esState.previousStepSize == ndir * stepSize);
211 BOOST_TEST(esState.stepSize == 1337.);
214 BOOST_TEST(esState.stepSize == -123.);
219 auto curvPars = std::get<0>(curvState);
224 BOOST_TEST(curvPars.covariance().has_value());
225 BOOST_TEST(*curvPars.covariance() != cov);
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);
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());
253 ps.stepping.covTransport =
false;
254 double h = es.
step(ps).value();
255 BOOST_TEST(ps.stepping.stepSize == h);
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());
264 ps.stepping.covTransport =
true;
265 double h2 = es.
step(ps).value();
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());
276 auto plane = Surface::makeShared<PlaneSurface>(
pos, mom.normalized());
279 stepSize, tolerance);
282 auto targetSurface = Surface::makeShared<PlaneSurface>(
283 pos + ndir * 2. * mom.normalized(), mom.normalized());
290 targetSurface->intersect(esState.geoContext, esState.pos,
291 esState.navDir * esState.dir,
false),
293 BOOST_TEST(esState.stepSize == 2.);
294 esState.stepSize = ndir * stepSize;
297 targetSurface->intersect(esState.geoContext, esState.pos,
298 esState.navDir * esState.dir,
false),
300 BOOST_TEST(esState.stepSize == 2.);
309 BOOST_TEST(boundPars.covariance().has_value());
310 BOOST_TEST(*boundPars.covariance() != cov);
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);
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());
359 [=](
const auto&
context,
const auto& inner,
const auto& vb) {
363 std::shared_ptr<const TrackingGeometry> vacuum =
374 Vector3D startParams(0., 0., 0.), startMom(1
_GeV, 0., 0.);
406 const auto& result = prop.propagate(sbtp, propOpts).value();
437 propDef(esDef, naviVac);
440 const auto& resultDef = propDef.propagate(sbtp, propOptsDef).value();
447 for (
unsigned int i = 0; i < stepResult.
position.size(); i++) {
451 for (
unsigned int i = 0; i < stepResult.
momentum.size(); i++) {
461 vConf.
volumeMaterial = std::make_shared<const HomogeneousVolumeMaterial>(
462 Material(352.8, 394.133, 9.012, 4., 1.848
e-3));
465 conf.position = {0.5_m, 0., 0.};
466 conf.length = {1
_m, 1
_m, 1_m};
472 [=](
const auto&
context,
const auto& inner,
const auto& vb) {
476 std::shared_ptr<const TrackingGeometry>
material =
487 Vector3D startParams(0., 0., 0.), startMom(5
_GeV, 0., 0.);
503 propOpts.
debug =
true;
520 const auto& result = prop.propagate(sbtp, propOpts).value();
540 BOOST_CHECK_LT(
mom.x(), 5
_GeV);
553 propOptsDense.
debug =
true;
561 propDense(esDense, naviMat);
564 const auto& resultDense = propDense.propagate(sbtp, propOptsDense).value();
571 for (
unsigned int i = 0; i < stepResult.
position.size(); i++) {
575 for (
unsigned int i = 0; i < stepResult.
momentum.size(); i++) {
585 StepperExtensionList<DefaultExtension, DenseEnvironmentExtension>,
586 detail::HighestValidAuctioneer>
590 DenseEnvironmentExtension>,
591 detail::HighestValidAuctioneer>,
595 const auto& resultB = propB.propagate(sbtp, propOptsDense).value();
613 BOOST_CHECK_NE(
mom.x(), 5
_GeV);
615 BOOST_CHECK_NE(
mom.z(), 0.);
623 vConfVac1.
position = {0.5_m, 0., 0.};
625 vConfVac1.
name =
"First vacuum volume";
627 vConfMat.
position = {1.5_m, 0., 0.};
629 vConfMat.
volumeMaterial = std::make_shared<const HomogeneousVolumeMaterial>(
630 Material(352.8, 394.133, 9.012, 4., 1.848
e-3));
631 vConfMat.
name =
"Material volume";
633 vConfVac2.
position = {2.5_m, 0., 0.};
635 vConfVac2.
name =
"Second vacuum volume";
637 conf.
volumeCfg = {vConfVac1, vConfMat, vConfVac2};
645 [=](
const auto&
context,
const auto& inner,
const auto& vb) {
659 Vector3D startParams(0., 0., 0.), startMom(5
_GeV, 0., 0.);
692 const auto& result = prop.propagate(sbtp, propOpts).value();
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()));
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()));
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()));
736 propOptsDef.
debug =
false;
744 propDef(esDef, naviDet);
747 const auto& resultDef =
748 propDef.propagate(sbtp, *(surs[0]), propOptsDef).value();
753 std::pair<Vector3D, Vector3D> endParams, endParamsControl;
754 for (
unsigned int i = 0; i < stepResultDef.
position.size(); i++) {
761 for (
unsigned int i = 0; i < stepResult.
position.size(); i++) {
769 if (propOptsDef.
debug) {
770 const auto debugString =
771 resultDef.template get<DebugOutput::result_type>().debugString;
772 std::cout << debugString << std::endl;
778 BOOST_TEST(endParams.first.x() == endParamsControl.first.x(),
779 tt::tolerance(1
e-5));
780 BOOST_TEST(endParams.first.y() == endParamsControl.first.y(),
781 tt::tolerance(1
e-5));
782 BOOST_TEST(endParams.first.z() == endParamsControl.first.z(),
783 tt::tolerance(1
e-5));
784 BOOST_TEST(endParams.second.x() == endParamsControl.second.x(),
785 tt::tolerance(1
e-5));
786 BOOST_TEST(endParams.second.y() == endParamsControl.second.y(),
787 tt::tolerance(1
e-5));
788 BOOST_TEST(endParams.second.z() == endParamsControl.second.z(),
789 tt::tolerance(1
e-5));
794 startParams = endParams.first;
795 startMom = endParams.second;
797 cov, startParams, startMom, 1., 0.);
815 propDense(esDense, naviDet);
818 const auto& resultDense =
819 propDense.propagate(sbtpPiecewise, *(surs[1]), propOptsDense).value();
824 for (
unsigned int i = 0; i < stepResultDense.
position.size(); i++) {
825 if (2
_m - stepResultDense.
position[i].x() < 1
e-4) {
826 endParams = std::make_pair(stepResultDense.
position[i],
831 for (
unsigned int i = 0; i < stepResult.
position.size(); i++) {
846 double rotationAngle =
M_PI * 0.5;
847 Vector3D xPos(cos(rotationAngle), 0., sin(rotationAngle));
849 Vector3D zPos(-sin(rotationAngle), 0., cos(rotationAngle));
860 sConf1.
surMat = std::shared_ptr<const ISurfaceMaterial>(
873 sConf2.
surMat = std::shared_ptr<const ISurfaceMaterial>(
881 muConf1.
length = {20._cm, 20._cm, 20._cm};
884 Material(352.8, 407., 9.012, 4., 1.848
e-3)));
885 muConf1.
name =
"MDT1";
888 muConf2.
length = {20._cm, 20._cm, 20._cm};
891 Material(352.8, 407., 9.012, 4., 1.848
e-3)));
892 muConf2.
name =
"MDT2";
896 vConf1.
length = {1._m, 1._m, 1._m};
898 vConf1.
name =
"Tracker";
901 vConf2.
length = {1._m, 1._m, 1._m};
904 Material(352.8, 407., 9.012, 4., 1.848
e-3)));
905 vConf2.
name =
"Calorimeter";
908 vConf3.
length = {1._m, 1._m, 1._m};
910 vConf3.
name =
"Muon system";
912 conf.
volumeCfg = {vConf1, vConf2, vConf3};
914 conf.
length = {3._m, 1._m, 1._m};
920 [=](
const auto&
context,
const auto& inner,
const auto& vb) {
924 std::shared_ptr<const TrackingGeometry>
detector =
935 Vector3D startParams(0., 0., 0.), startMom(1.
_GeV, 0., 0.);
960 const auto& result = prop.propagate(sbtp, propOpts).value();
965 double lastMomentum = stepResult.
momentum[0].x();
966 for (
unsigned int i = 0; i < stepResult.
position.size(); i++) {
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();
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);