Tests/IntegrationTests/NavigatorConsistency.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/data/test_case.hpp>
0010 #include <boost/test/unit_test.hpp>
0011
0012 #include "Acts/Definitions/Algebra.hpp"
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/EventData/TrackParameters.hpp"
0015 #include "Acts/Geometry/GeometryContext.hpp"
0016 #include "Acts/MagneticField/ConstantBField.hpp"
0017 #include "Acts/Propagator/EigenStepper.hpp"
0018 #include "Acts/Propagator/Navigator.hpp"
0019 #include "Acts/Propagator/Propagator.hpp"
0020 #include "Acts/Propagator/StraightLineStepper.hpp"
0021 #include "Acts/Propagator/SurfaceCollector.hpp"
0022 #include "Acts/Propagator/TryAllNavigator.hpp"
0023 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0024 #include "Acts/Utilities/Logger.hpp"
0025 #include "Acts/Utilities/VectorHelpers.hpp"
0026 #include "ActsTests/CommonHelpers/CylindricalTrackingGeometry.hpp"
0027
0028 #include <algorithm>
0029 #include <numbers>
0030
0031 namespace bdata = boost::unit_test::data;
0032
0033 using namespace Acts;
0034 using namespace Acts::UnitLiterals;
0035 using Acts::VectorHelpers::perp;
0036
0037 namespace ActsTests {
0038
0039 // Create a test context
0040 GeometryContext tgContext = GeometryContext();
0041 MagneticFieldContext mfContext = MagneticFieldContext();
0042
0043 ActsTests::CylindricalTrackingGeometry cGeometry(tgContext);
0044 auto tGeometry = cGeometry();
0045
0046 const double Bz = 2_T;
0047 auto bField = std::make_shared<ConstantBField>(Vector3{0, 0, Bz});
0048
0049 using TestSurfaceCollector = SurfaceCollector<SurfaceSelector>;
0050
0051 std::vector<GeometryIdentifier> collectRelevantGeoIds(
0052     const TestSurfaceCollector::result_type& surfaceHits) {
0053   std::vector<GeometryIdentifier> geoIds;
0054   for (const auto& surfaceHit : surfaceHits.collected) {
0055     auto geoId = surfaceHit.surface->geometryId();
0056     auto material = surfaceHit.surface->surfaceMaterial();
0057     if (geoId.sensitive() == 0 && material == nullptr) {
0058       continue;
0059     }
0060     geoIds.push_back(geoId);
0061   }
0062   return geoIds;
0063 }
0064
0065 /// the actual test method that runs the test can be used with several
0066 /// propagator types
0067 ///
0068 /// @tparam propagator_t is the actual propagator type
0069 ///
0070 /// @param prop is the propagator instance
0071 /// @param start start parameters for propagation
0072 /// @param logger A logger instance
0073 template <typename propagator_t>
0074 void runSelfConsistencyTest(const propagator_t& prop,
0075                             const BoundTrackParameters& start,
0076                             const Logger& logger) {
0077   // Actor list
0078   using ActorList = ActorList<TestSurfaceCollector>;
0079   using Options = typename propagator_t::template Options<ActorList>;
0080
0081   // forward surface test
0082   Options fwdOptions(tgContext, mfContext);
0083   fwdOptions.pathLimit = 25_cm;
0084
0085   // get the surface collector and configure it
0086   auto& fwdSurfaceCollector =
0087       fwdOptions.actorList.template get<TestSurfaceCollector>();
0088   fwdSurfaceCollector.selector.selectSensitive = true;
0089   fwdSurfaceCollector.selector.selectMaterial = true;
0090   fwdSurfaceCollector.selector.selectPassive = true;
0091
0092   ACTS_DEBUG(">>> Forward Propagation : start.");
0093   auto fwdResult = prop.propagate(start, fwdOptions).value();
0094   auto fwdSurfaceHits =
0095       fwdResult.template get<TestSurfaceCollector::result_type>().collected;
0096   auto fwdSurfaces = collectRelevantGeoIds(
0097       fwdResult.template get<TestSurfaceCollector::result_type>());
0098
0099   ACTS_DEBUG(">>> Surface hits found on ...");
0100   for (const auto& fwdSteps : fwdSurfaces) {
0101     ACTS_DEBUG("--> Surface with " << fwdSteps);
0102   }
0103   ACTS_DEBUG(">>> Forward Propagation : end.");
0104
0105   // backward surface test
0106   Options bwdOptions(tgContext, mfContext);
0107   bwdOptions.pathLimit = 25_cm;
0108   bwdOptions.direction = Direction::Backward();
0109
0110   // get the surface collector and configure it
0111   auto& bwdMSurfaceCollector =
0112       bwdOptions.actorList.template get<TestSurfaceCollector>();
0113   bwdMSurfaceCollector.selector.selectSensitive = true;
0114   bwdMSurfaceCollector.selector.selectMaterial = true;
0115   bwdMSurfaceCollector.selector.selectPassive = true;
0116
0117   const auto& startSurface = start.referenceSurface();
0118
0119   ACTS_DEBUG(">>> Backward Propagation : start.");
0120   auto bwdResult =
0121       prop.propagate(*fwdResult.endParameters, startSurface, bwdOptions)
0122           .value();
0123   auto bwdSurfaceHits =
0124       bwdResult.template get<TestSurfaceCollector::result_type>().collected;
0125   auto bwdSurfaces = collectRelevantGeoIds(
0126       bwdResult.template get<TestSurfaceCollector::result_type>());
0127
0128   ACTS_DEBUG(">>> Surface hits found on ...");
0129   for (auto& bwdSteps : bwdSurfaces) {
0130     ACTS_DEBUG("--> Surface with " << bwdSteps);
0131   }
0132   ACTS_DEBUG(">>> Backward Propagation : end.");
0133
0134   // forward-backward compatibility test
0135   {
0136     std::ranges::reverse(bwdSurfaces);
0137     BOOST_CHECK_EQUAL_COLLECTIONS(bwdSurfaces.begin(), bwdSurfaces.end(),
0138                                   fwdSurfaces.begin(), fwdSurfaces.end());
0139   }
0140
0141   // stepping from one surface to the next
0142   // now go from surface to surface and check
0143   Options fwdStepOptions(tgContext, mfContext);
0144
0145   // get the surface collector and configure it
0146   auto& fwdStepSurfaceCollector =
0147       fwdOptions.actorList.template get<TestSurfaceCollector>();
0148   fwdStepSurfaceCollector.selector.selectSensitive = true;
0149   fwdStepSurfaceCollector.selector.selectMaterial = true;
0150   fwdStepSurfaceCollector.selector.selectPassive = true;
0151
0152   std::vector<GeometryIdentifier> fwdStepSurfaces;
0153
0154   // move forward step by step through the surfaces
0155   BoundTrackParameters sParameters = start;
0156   std::vector<BoundTrackParameters> stepParameters;
0157   for (auto& fwdSteps : fwdSurfaceHits) {
0158     ACTS_DEBUG(">>> Forward step : "
0159                << sParameters.referenceSurface().geometryId() << " --> "
0160                << fwdSteps.surface->geometryId());
0161
0162     // make a forward step
0163     auto fwdStep =
0164         prop.propagate(sParameters, *fwdSteps.surface, fwdStepOptions).value();
0165
0166     auto fwdStepSurfacesTmp = collectRelevantGeoIds(
0167         fwdStep.template get<TestSurfaceCollector::result_type>());
0168     fwdStepSurfaces.insert(fwdStepSurfaces.end(), fwdStepSurfacesTmp.begin(),
0169                            fwdStepSurfacesTmp.end());
0170
0171     if (fwdStep.endParameters.has_value()) {
0172       // make sure the parameters do not run out of scope
0173       stepParameters.push_back(*fwdStep.endParameters);
0174       sParameters = stepParameters.back();
0175     }
0176   }
0177   // final destination surface
0178   const Surface& dSurface = fwdResult.endParameters->referenceSurface();
0179   ACTS_DEBUG(">>> Forward step : "
0180              << sParameters.referenceSurface().geometryId() << " --> "
0181              << dSurface.geometryId());
0182   auto fwdStepFinal =
0183       prop.propagate(sParameters, dSurface, fwdStepOptions).value();
0184   auto fwdStepSurfacesTmp = collectRelevantGeoIds(
0185       fwdStepFinal.template get<TestSurfaceCollector::result_type>());
0186   fwdStepSurfaces.insert(fwdStepSurfaces.end(), fwdStepSurfacesTmp.begin(),
0187                          fwdStepSurfacesTmp.end());
0188
0189   // TODO forward-forward step compatibility test
0190
0191   // stepping from one surface to the next : backwards
0192   // now go from surface to surface and check
0193   Options bwdStepOptions(tgContext, mfContext);
0194   bwdStepOptions.direction = Direction::Backward();
0195
0196   // get the surface collector and configure it
0197   auto& bwdStepSurfaceCollector =
0198       bwdOptions.actorList.template get<TestSurfaceCollector>();
0199   bwdStepSurfaceCollector.selector.selectSensitive = true;
0200   bwdStepSurfaceCollector.selector.selectMaterial = true;
0201   bwdStepSurfaceCollector.selector.selectPassive = true;
0202
0203   std::vector<GeometryIdentifier> bwdStepSurfaces;
0204
0205   // move forward step by step through the surfaces
0206   sParameters = *fwdResult.endParameters;
0207   for (auto& bwdSteps : bwdSurfaceHits) {
0208     ACTS_DEBUG(">>> Backward step : "
0209                << sParameters.referenceSurface().geometryId() << " --> "
0210                << bwdSteps.surface->geometryId());
0211
0212     // make a forward step
0213     auto bwdStep =
0214         prop.propagate(sParameters, *bwdSteps.surface, bwdStepOptions).value();
0215
0216     auto bwdStepSurfacesTmp = collectRelevantGeoIds(
0217         bwdStep.template get<TestSurfaceCollector::result_type>());
0218     bwdStepSurfaces.insert(bwdStepSurfaces.end(), bwdStepSurfacesTmp.begin(),
0219                            bwdStepSurfacesTmp.end());
0220
0221     if (bwdStep.endParameters.has_value()) {
0222       // make sure the parameters do not run out of scope
0223       stepParameters.push_back(*bwdStep.endParameters);
0224       sParameters = stepParameters.back();
0225     }
0226   }
0227   // final destination surface
0228   const Surface& dbSurface = start.referenceSurface();
0229   ACTS_DEBUG(">>> Backward step : "
0230              << sParameters.referenceSurface().geometryId() << " --> "
0231              << dSurface.geometryId());
0232   auto bwdStepFinal =
0233       prop.propagate(sParameters, dbSurface, bwdStepOptions).value();
0234   auto bwdStepSurfacesTmp = collectRelevantGeoIds(
0235       bwdStepFinal.template get<TestSurfaceCollector::result_type>());
0236   bwdStepSurfaces.insert(bwdStepSurfaces.end(), bwdStepSurfacesTmp.begin(),
0237                          bwdStepSurfacesTmp.end());
0238
0239   // TODO backward-backward step compatibility test
0240
0241   std::ranges::reverse(bwdStepSurfaces);
0242   BOOST_CHECK_EQUAL_COLLECTIONS(bwdStepSurfaces.begin(), bwdStepSurfaces.end(),
0243                                 fwdStepSurfaces.begin(), fwdStepSurfaces.end());
0244 }
0245
0246 /// the actual test method that runs the test can be used with several
0247 /// propagator types
0248 ///
0249 /// @tparam propagator_probe_t is the probe propagator type
0250 /// @tparam propagator_ref_t is the reference propagator type
0251 ///
0252 /// @param propProbe is the probe propagator instance
0253 /// @param propRef is the reference propagator instance
0254 /// @param start start parameters for propagation
0255 /// @param logger A logger instance
0256 template <typename propagator_probe_t, typename propagator_ref_t>
0257 void runConsistencyTest(const propagator_probe_t& propProbe,
0258                         const propagator_ref_t& propRef,
0259                         const BoundTrackParameters& start,
0260                         const Logger& logger) {
0261   // Action list and abort list
0262   using ActorList = ActorList<TestSurfaceCollector>;
0263
0264   auto run = [&](const auto& prop) {
0265     using propagator_t = std::decay_t<decltype(prop)>;
0266     using Options = typename propagator_t::template Options<ActorList>;
0267
0268     // forward surface test
0269     Options fwdOptions(tgContext, mfContext);
0270     fwdOptions.pathLimit = 25_cm;
0271     fwdOptions.stepping.maxStepSize = 1_cm;
0272
0273     // get the surface collector and configure it
0274     auto& fwdSurfaceCollector =
0275         fwdOptions.actorList.template get<TestSurfaceCollector>();
0276     fwdSurfaceCollector.selector.selectSensitive = true;
0277     fwdSurfaceCollector.selector.selectMaterial = true;
0278     fwdSurfaceCollector.selector.selectPassive = true;
0279
0280     auto fwdResult = prop.propagate(start, fwdOptions).value();
0281     auto fwdSurfaces = collectRelevantGeoIds(
0282         fwdResult.template get<TestSurfaceCollector::result_type>());
0283
0284     ACTS_DEBUG(">>> Surface hits found on ...");
0285     for (const auto& fwdSteps : fwdSurfaces) {
0286       ACTS_DEBUG("--> Surface with " << fwdSteps);
0287     }
0288
0289     return fwdSurfaces;
0290   };
0291
0292   ACTS_DEBUG(">>> Probe Propagation : start.");
0293   const auto& probeSurfaces = run(propProbe);
0294   ACTS_DEBUG(">>> Probe Propagation : end.");
0295
0296   ACTS_DEBUG(">>> Reference Propagation : start.");
0297   const auto& refSurfaces = run(propRef);
0298   ACTS_DEBUG(">>> Reference Propagation : end.");
0299
0300   // probe-ref compatibility test
0301   BOOST_CHECK_EQUAL_COLLECTIONS(probeSurfaces.begin(), probeSurfaces.end(),
0302                                 refSurfaces.begin(), refSurfaces.end());
0303 }
0304
0305 Logging::Level logLevel = Logging::INFO;
0306
0307 const int nTestsSelfConsistency = 500;
0308 const int nTestsRefConsistency = 500;
0309
0310 using StraightLinePropagator = Propagator<StraightLineStepper, Navigator>;
0311 using TestEigenStepper = EigenStepper<>;
0312 using EigenPropagator = Propagator<TestEigenStepper, Navigator>;
0313 using Reference1StraightLinePropagator =
0314     Propagator<StraightLineStepper, TryAllNavigator>;
0315 using Reference1EigenPropagator = Propagator<TestEigenStepper, TryAllNavigator>;
0316 using Reference2StraightLinePropagator =
0317     Propagator<StraightLineStepper, TryAllOverstepNavigator>;
0318 using Reference2EigenPropagator =
0319     Propagator<TestEigenStepper, TryAllOverstepNavigator>;
0320
0321 StraightLineStepper slstepper;
0322 TestEigenStepper estepper(bField);
0323
0324 StraightLinePropagator slpropagator(slstepper,
0325                                     Navigator({tGeometry, true, true, false},
0326                                               getDefaultLogger("sl_nav",
0327                                                                Logging::INFO)),
0328                                     getDefaultLogger("sl_prop", Logging::INFO));
0329 EigenPropagator epropagator(estepper,
0330                             Navigator({tGeometry, true, true, false},
0331                                       getDefaultLogger("e_nav", Logging::INFO)),
0332                             getDefaultLogger("e_prop", Logging::INFO));
0333
0334 Reference1StraightLinePropagator refslpropagator1(
0335     slstepper,
0336     TryAllNavigator({tGeometry, true, true, false},
0337                     getDefaultLogger("ref1_sl_nav", Logging::INFO)),
0338     getDefaultLogger("ref1_sl_prop", Logging::INFO));
0339 Reference1EigenPropagator refepropagator1(
0340     estepper,
0341     TryAllNavigator({tGeometry, true, true, false,
0342                      BoundaryTolerance::Infinite()},
0343                     getDefaultLogger("ref1_e_nav", Logging::INFO)),
0344     getDefaultLogger("ref1_e_prop", Logging::INFO));
0345
0346 Reference2EigenPropagator refepropagator2(
0347     estepper,
0348     TryAllOverstepNavigator({tGeometry, true, true, false,
0349                              BoundaryTolerance::Infinite()},
0350                             getDefaultLogger("ref2_e_nav", Logging::INFO)),
0351     getDefaultLogger("ref2_e_prop", Logging::INFO));
0352 Reference2StraightLinePropagator refslpropagator2(
0353     slstepper,
0354     TryAllOverstepNavigator({tGeometry, true, true, false},
0355                             getDefaultLogger("ref2_sl_nav", Logging::INFO)),
0356     getDefaultLogger("ref2_sl_prop", Logging::INFO));
0357
0358 auto eventGen =
0359     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 20,
0360                    bdata::distribution = std::uniform_real_distribution<double>(
0361                        0.5_GeV, 10_GeV))) ^
0362     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 21,
0363                    bdata::distribution = std::uniform_real_distribution<double>(
0364                        -std::numbers::pi, std::numbers::pi))) ^
0365     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 22,
0366                    bdata::distribution = std::uniform_real_distribution<double>(
0367                        1., std::numbers::pi - 1.))) ^
0368     bdata::random(
0369         (bdata::engine = std::mt19937(), bdata::seed = 23,
0370          bdata::distribution = std::uniform_int_distribution<int>(0, 1)));
0371
0372 BoundTrackParameters createStartParameters(double pT, double phi, double theta,
0373                                            int charge) {
0374   double p = pT / std::sin(theta);
0375   double q = -1 + 2 * charge;
0376   return BoundTrackParameters::createCurvilinear(Vector4::Zero(), phi, theta,
0377                                                  q / p, std::nullopt,
0378                                                  ParticleHypothesis::pion());
0379 }
0380
0381 BOOST_DATA_TEST_CASE(NavigatorStraightLineSelfConsistency,
0382                      eventGen ^ bdata::xrange(nTestsSelfConsistency), pT, phi,
0383                      theta, charge, index) {
0384   ACTS_LOCAL_LOGGER(getDefaultLogger("NavigatorTest", logLevel));
0385
0386   BoundTrackParameters start = createStartParameters(pT, phi, theta, charge);
0387
0388   ACTS_DEBUG(">>> Run navigation tests with:\n    pT = "
0389              << pT << "\n    phi = " << phi << "\n    theta = " << theta
0390              << "\n    charge = " << charge << "\n    index = " << index);
0391
0392   ACTS_DEBUG(">>> Test self consistency slpropagator");
0393   runSelfConsistencyTest(slpropagator, start, logger());
0394 }
0395
0396 BOOST_DATA_TEST_CASE(NavigatorEigenSelfConsistency,
0397                      eventGen ^ bdata::xrange(nTestsSelfConsistency), pT, phi,
0398                      theta, charge, index) {
0399   ACTS_LOCAL_LOGGER(getDefaultLogger("NavigatorTest", logLevel));
0400
0401   BoundTrackParameters start = createStartParameters(pT, phi, theta, charge);
0402
0403   ACTS_DEBUG(">>> Run navigation tests with:\n    pT = "
0404              << pT << "\n    phi = " << phi << "\n    theta = " << theta
0405              << "\n    charge = " << charge << "\n    index = " << index);
0406
0407   ACTS_DEBUG(">>> Test self consistency epropagator");
0408   runSelfConsistencyTest(epropagator, start, logger());
0409 }
0410
0411 BOOST_DATA_TEST_CASE(NavigatorRef1StraightLineConsistency,
0412                      eventGen ^ bdata::xrange(nTestsRefConsistency), pT, phi,
0413                      theta, charge, index) {
0414   ACTS_LOCAL_LOGGER(getDefaultLogger("NavigatorTest", logLevel));
0415
0416   BoundTrackParameters start = createStartParameters(pT, phi, theta, charge);
0417
0418   ACTS_DEBUG(">>> Run navigation tests with:\n    pT = "
0419              << pT << "\n    phi = " << phi << "\n    theta = " << theta
0420              << "\n    charge = " << charge << "\n    index = " << index);
0421
0422   ACTS_DEBUG(">>> Test reference 1 consistency slpropagator");
0423   runConsistencyTest(slpropagator, refslpropagator1, start, logger());
0424 }
0425
0426 BOOST_DATA_TEST_CASE(NavigatorRef1EigenConsistency,
0427                      eventGen ^ bdata::xrange(nTestsRefConsistency), pT, phi,
0428                      theta, charge, index) {
0429   ACTS_LOCAL_LOGGER(getDefaultLogger("NavigatorTest", logLevel));
0430
0431   BoundTrackParameters start = createStartParameters(pT, phi, theta, charge);
0432
0433   ACTS_DEBUG(">>> Run navigation tests with:\n    pT = "
0434              << pT << "\n    phi = " << phi << "\n    theta = " << theta
0435              << "\n    charge = " << charge << "\n    index = " << index);
0436
0437   ACTS_DEBUG(">>> Test reference 1 consistency epropagator");
0438   runConsistencyTest(epropagator, refepropagator1, start, logger());
0439 }
0440
0441 BOOST_DATA_TEST_CASE(NavigatorRef2StraightLineConsistency,
0442                      eventGen ^ bdata::xrange(nTestsRefConsistency), pT, phi,
0443                      theta, charge, index) {
0444   ACTS_LOCAL_LOGGER(getDefaultLogger("NavigatorTest", logLevel));
0445
0446   BoundTrackParameters start = createStartParameters(pT, phi, theta, charge);
0447
0448   ACTS_DEBUG(">>> Run navigation tests with:\n    pT = "
0449              << pT << "\n    phi = " << phi << "\n    theta = " << theta
0450              << "\n    charge = " << charge << "\n    index = " << index);
0451
0452   ACTS_DEBUG(">>> Test reference 2 consistency slpropagator");
0453   runConsistencyTest(slpropagator, refslpropagator2, start, logger());
0454 }
0455
0456 BOOST_DATA_TEST_CASE(NavigatorRef2EigenConsistency,
0457                      eventGen ^ bdata::xrange(nTestsRefConsistency), pT, phi,
0458                      theta, charge, index) {
0459   ACTS_LOCAL_LOGGER(getDefaultLogger("NavigatorTest", logLevel));
0460
0461   BoundTrackParameters start = createStartParameters(pT, phi, theta, charge);
0462
0463   ACTS_DEBUG(">>> Run navigation tests with:\n    pT = "
0464              << pT << "\n    phi = " << phi << "\n    theta = " << theta
0465              << "\n    charge = " << charge << "\n    index = " << index);
0466
0467   ACTS_DEBUG(">>> Test reference 2 consistency epropagator");
0468   runConsistencyTest(epropagator, refepropagator2, start, logger());
0469 }
0470
0471 }  // namespace ActsTests