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