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 using ParamRecorder = BoundParameterRecorder<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, ParamRecorder>;
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   auto& paramRecorder = fwdOptions.actorList.template get<ParamRecorder>();
0093   paramRecorder.selector.selectSensitive = true;
0094   paramRecorder.selector.selectMaterial = true;
0095   paramRecorder.selector.selectPassive = true;
0096
0097   ACTS_DEBUG(">>> Forward Propagation : start.");
0098   auto fwdResult = prop.propagate(start, fwdOptions).value();
0099   auto fwdSurfaceHits =
0100       fwdResult.template get<TestSurfaceCollector::result_type>().collected;
0101
0102   auto fwdTrackRecords = fwdResult.template get<ParamRecorder::result_type>();
0103
0104   BOOST_CHECK_EQUAL(fwdTrackRecords.size(), fwdSurfaceHits.size());
0105   for (const auto& boundPars :
0106        fwdResult.template get<ParamRecorder::result_type>()) {
0107     BOOST_CHECK_EQUAL(boundPars.particleHypothesis(),
0108                       start.particleHypothesis());
0109     BOOST_CHECK_EQUAL(boundPars.covariance().has_value(),
0110                       start.covariance().has_value());
0111     BOOST_CHECK_EQUAL(
0112         std::ranges::any_of(
0113             fwdSurfaceHits,
0114             [&](const Acts::SurfaceHit& hit) {
0115               return hit.surface ==
0116                      boundPars.referenceSurface().getSharedPtr().get();
0117             }),
0118         true);
0119   }
0120
0121   auto fwdSurfaces = collectRelevantGeoIds(
0122       fwdResult.template get<TestSurfaceCollector::result_type>());
0123
0124   ACTS_DEBUG(">>> Surface hits found on ...");
0125   for (const auto& fwdSteps : fwdSurfaces) {
0126     ACTS_DEBUG("--> Surface with " << fwdSteps);
0127   }
0128   ACTS_DEBUG(">>> Forward Propagation : end.");
0129
0130   // backward surface test
0131   Options bwdOptions(tgContext, mfContext);
0132   bwdOptions.pathLimit = 25_cm;
0133   bwdOptions.direction = Direction::Backward();
0134
0135   // get the surface collector and configure it
0136   auto& bwdMSurfaceCollector =
0137       bwdOptions.actorList.template get<TestSurfaceCollector>();
0138   bwdMSurfaceCollector.selector.selectSensitive = true;
0139   bwdMSurfaceCollector.selector.selectMaterial = true;
0140   bwdMSurfaceCollector.selector.selectPassive = true;
0141
0142   const auto& startSurface = start.referenceSurface();
0143
0144   ACTS_DEBUG(">>> Backward Propagation : start.");
0145   auto bwdResult =
0146       prop.propagate(*fwdResult.endParameters, startSurface, bwdOptions)
0147           .value();
0148   auto bwdSurfaceHits =
0149       bwdResult.template get<TestSurfaceCollector::result_type>().collected;
0150   auto bwdSurfaces = collectRelevantGeoIds(
0151       bwdResult.template get<TestSurfaceCollector::result_type>());
0152
0153   ACTS_DEBUG(">>> Surface hits found on ...");
0154   for (auto& bwdSteps : bwdSurfaces) {
0155     ACTS_DEBUG("--> Surface with " << bwdSteps);
0156   }
0157   ACTS_DEBUG(">>> Backward Propagation : end.");
0158
0159   // forward-backward compatibility test
0160   {
0161     std::ranges::reverse(bwdSurfaces);
0162     BOOST_CHECK_EQUAL_COLLECTIONS(bwdSurfaces.begin(), bwdSurfaces.end(),
0163                                   fwdSurfaces.begin(), fwdSurfaces.end());
0164   }
0165
0166   // stepping from one surface to the next
0167   // now go from surface to surface and check
0168   Options fwdStepOptions(tgContext, mfContext);
0169
0170   // get the surface collector and configure it
0171   auto& fwdStepSurfaceCollector =
0172       fwdOptions.actorList.template get<TestSurfaceCollector>();
0173   fwdStepSurfaceCollector.selector.selectSensitive = true;
0174   fwdStepSurfaceCollector.selector.selectMaterial = true;
0175   fwdStepSurfaceCollector.selector.selectPassive = true;
0176
0177   std::vector<GeometryIdentifier> fwdStepSurfaces;
0178
0179   // move forward step by step through the surfaces
0180   BoundTrackParameters sParameters = start;
0181   std::vector<BoundTrackParameters> stepParameters;
0182   for (auto& fwdSteps : fwdSurfaceHits) {
0183     ACTS_DEBUG(">>> Forward step : "
0184                << sParameters.referenceSurface().geometryId() << " --> "
0185                << fwdSteps.surface->geometryId());
0186
0187     // make a forward step
0188     auto fwdStep =
0189         prop.propagate(sParameters, *fwdSteps.surface, fwdStepOptions).value();
0190
0191     auto fwdStepSurfacesTmp = collectRelevantGeoIds(
0192         fwdStep.template get<TestSurfaceCollector::result_type>());
0193     fwdStepSurfaces.insert(fwdStepSurfaces.end(), fwdStepSurfacesTmp.begin(),
0194                            fwdStepSurfacesTmp.end());
0195
0196     if (fwdStep.endParameters.has_value()) {
0197       // make sure the parameters do not run out of scope
0198       stepParameters.push_back(*fwdStep.endParameters);
0199       sParameters = stepParameters.back();
0200     }
0201   }
0202   // final destination surface
0203   const Surface& dSurface = fwdResult.endParameters->referenceSurface();
0204   ACTS_DEBUG(">>> Forward step : "
0205              << sParameters.referenceSurface().geometryId() << " --> "
0206              << dSurface.geometryId());
0207   auto fwdStepFinal =
0208       prop.propagate(sParameters, dSurface, fwdStepOptions).value();
0209   auto fwdStepSurfacesTmp = collectRelevantGeoIds(
0210       fwdStepFinal.template get<TestSurfaceCollector::result_type>());
0211   fwdStepSurfaces.insert(fwdStepSurfaces.end(), fwdStepSurfacesTmp.begin(),
0212                          fwdStepSurfacesTmp.end());
0213
0214   // TODO forward-forward step compatibility test
0215
0216   // stepping from one surface to the next : backwards
0217   // now go from surface to surface and check
0218   Options bwdStepOptions(tgContext, mfContext);
0219   bwdStepOptions.direction = Direction::Backward();
0220
0221   // get the surface collector and configure it
0222   auto& bwdStepSurfaceCollector =
0223       bwdOptions.actorList.template get<TestSurfaceCollector>();
0224   bwdStepSurfaceCollector.selector.selectSensitive = true;
0225   bwdStepSurfaceCollector.selector.selectMaterial = true;
0226   bwdStepSurfaceCollector.selector.selectPassive = true;
0227
0228   std::vector<GeometryIdentifier> bwdStepSurfaces;
0229
0230   // move forward step by step through the surfaces
0231   sParameters = *fwdResult.endParameters;
0232   for (auto& bwdSteps : bwdSurfaceHits) {
0233     ACTS_DEBUG(">>> Backward step : "
0234                << sParameters.referenceSurface().geometryId() << " --> "
0235                << bwdSteps.surface->geometryId());
0236
0237     // make a forward step
0238     auto bwdStep =
0239         prop.propagate(sParameters, *bwdSteps.surface, bwdStepOptions).value();
0240
0241     auto bwdStepSurfacesTmp = collectRelevantGeoIds(
0242         bwdStep.template get<TestSurfaceCollector::result_type>());
0243     bwdStepSurfaces.insert(bwdStepSurfaces.end(), bwdStepSurfacesTmp.begin(),
0244                            bwdStepSurfacesTmp.end());
0245
0246     if (bwdStep.endParameters.has_value()) {
0247       // make sure the parameters do not run out of scope
0248       stepParameters.push_back(*bwdStep.endParameters);
0249       sParameters = stepParameters.back();
0250     }
0251   }
0252   // final destination surface
0253   const Surface& dbSurface = start.referenceSurface();
0254   ACTS_DEBUG(">>> Backward step : "
0255              << sParameters.referenceSurface().geometryId() << " --> "
0256              << dSurface.geometryId());
0257   auto bwdStepFinal =
0258       prop.propagate(sParameters, dbSurface, bwdStepOptions).value();
0259   auto bwdStepSurfacesTmp = collectRelevantGeoIds(
0260       bwdStepFinal.template get<TestSurfaceCollector::result_type>());
0261   bwdStepSurfaces.insert(bwdStepSurfaces.end(), bwdStepSurfacesTmp.begin(),
0262                          bwdStepSurfacesTmp.end());
0263
0264   // TODO backward-backward step compatibility test
0265
0266   std::ranges::reverse(bwdStepSurfaces);
0267   BOOST_CHECK_EQUAL_COLLECTIONS(bwdStepSurfaces.begin(), bwdStepSurfaces.end(),
0268                                 fwdStepSurfaces.begin(), fwdStepSurfaces.end());
0269 }
0270
0271 /// the actual test method that runs the test can be used with several
0272 /// propagator types
0273 ///
0274 /// @tparam propagator_probe_t is the probe propagator type
0275 /// @tparam propagator_ref_t is the reference propagator type
0276 ///
0277 /// @param propProbe is the probe propagator instance
0278 /// @param propRef is the reference propagator instance
0279 /// @param start start parameters for propagation
0280 /// @param logger A logger instance
0281 template <typename propagator_probe_t, typename propagator_ref_t>
0282 void runConsistencyTest(const propagator_probe_t& propProbe,
0283                         const propagator_ref_t& propRef,
0284                         const BoundTrackParameters& start,
0285                         const Logger& logger) {
0286   // Action list and abort list
0287   using ActorList = ActorList<TestSurfaceCollector>;
0288
0289   auto run = [&](const auto& prop) {
0290     using propagator_t = std::decay_t<decltype(prop)>;
0291     using Options = typename propagator_t::template Options<ActorList>;
0292
0293     // forward surface test
0294     Options fwdOptions(tgContext, mfContext);
0295     fwdOptions.pathLimit = 25_cm;
0296     fwdOptions.stepping.maxStepSize = 1_cm;
0297
0298     // get the surface collector and configure it
0299     auto& fwdSurfaceCollector =
0300         fwdOptions.actorList.template get<TestSurfaceCollector>();
0301     fwdSurfaceCollector.selector.selectSensitive = true;
0302     fwdSurfaceCollector.selector.selectMaterial = true;
0303     fwdSurfaceCollector.selector.selectPassive = true;
0304
0305     auto fwdResult = prop.propagate(start, fwdOptions).value();
0306     auto fwdSurfaces = collectRelevantGeoIds(
0307         fwdResult.template get<TestSurfaceCollector::result_type>());
0308
0309     ACTS_DEBUG(">>> Surface hits found on ...");
0310     for (const auto& fwdSteps : fwdSurfaces) {
0311       ACTS_DEBUG("--> Surface with " << fwdSteps);
0312     }
0313
0314     return fwdSurfaces;
0315   };
0316
0317   ACTS_DEBUG(">>> Probe Propagation : start.");
0318   const auto& probeSurfaces = run(propProbe);
0319   ACTS_DEBUG(">>> Probe Propagation : end.");
0320
0321   ACTS_DEBUG(">>> Reference Propagation : start.");
0322   const auto& refSurfaces = run(propRef);
0323   ACTS_DEBUG(">>> Reference Propagation : end.");
0324
0325   // probe-ref compatibility test
0326   BOOST_CHECK_EQUAL_COLLECTIONS(probeSurfaces.begin(), probeSurfaces.end(),
0327                                 refSurfaces.begin(), refSurfaces.end());
0328 }
0329
0330 Logging::Level logLevel = Logging::INFO;
0331
0332 const int nTestsSelfConsistency = 500;
0333 const int nTestsRefConsistency = 500;
0334
0335 using StraightLinePropagator = Propagator<StraightLineStepper, Navigator>;
0336 using TestEigenStepper = EigenStepper<>;
0337 using EigenPropagator = Propagator<TestEigenStepper, Navigator>;
0338 using ReferenceStraightLinePropagator =
0339     Propagator<StraightLineStepper, Experimental::TryAllNavigator>;
0340 using ReferenceEigenPropagator =
0341     Propagator<TestEigenStepper, Experimental::TryAllNavigator>;
0342
0343 StraightLineStepper slstepper;
0344 TestEigenStepper estepper(bField);
0345
0346 StraightLinePropagator slpropagator(slstepper,
0347                                     Navigator({tGeometry, true, true, false},
0348                                               getDefaultLogger("sl_nav",
0349                                                                Logging::INFO)),
0350                                     getDefaultLogger("sl_prop", Logging::INFO));
0351 EigenPropagator epropagator(estepper,
0352                             Navigator({tGeometry, true, true, false},
0353                                       getDefaultLogger("e_nav", Logging::INFO)),
0354                             getDefaultLogger("e_prop", Logging::INFO));
0355
0356 ReferenceStraightLinePropagator refslpropagator(
0357     slstepper,
0358     Experimental::TryAllNavigator({tGeometry, true, true, false},
0359                                   getDefaultLogger("ref_sl_nav",
0360                                                    Logging::INFO)),
0361     getDefaultLogger("ref_sl_prop", Logging::INFO));
0362 ReferenceEigenPropagator refepropagator(
0363     estepper,
0364     Experimental::TryAllNavigator({tGeometry, true, true, false,
0365                                    BoundaryTolerance::Infinite()},
0366                                   getDefaultLogger("ref_e_nav", Logging::INFO)),
0367     getDefaultLogger("ref_e_prop", Logging::INFO));
0368
0369 auto eventGen =
0370     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 20,
0371                    bdata::distribution = std::uniform_real_distribution<double>(
0372                        0.5_GeV, 10_GeV))) ^
0373     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 21,
0374                    bdata::distribution = std::uniform_real_distribution<double>(
0375                        -std::numbers::pi, std::numbers::pi))) ^
0376     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 22,
0377                    bdata::distribution = std::uniform_real_distribution<double>(
0378                        1., std::numbers::pi - 1.))) ^
0379     bdata::random(
0380         (bdata::engine = std::mt19937(), bdata::seed = 23,
0381          bdata::distribution = std::uniform_int_distribution<int>(0, 1)));
0382
0383 BoundTrackParameters createStartParameters(double pT, double phi, double theta,
0384                                            int charge) {
0385   double p = pT / std::sin(theta);
0386   double q = -1 + 2 * charge;
0387   return BoundTrackParameters::createCurvilinear(Vector4::Zero(), phi, theta,
0388                                                  q / p, std::nullopt,
0389                                                  ParticleHypothesis::pion());
0390 }
0391
0392 BOOST_DATA_TEST_CASE(NavigatorStraightLineSelfConsistency,
0393                      eventGen ^ bdata::xrange(nTestsSelfConsistency), pT, phi,
0394                      theta, charge, index) {
0395   ACTS_LOCAL_LOGGER(getDefaultLogger("NavigatorTest", logLevel));
0396
0397   BoundTrackParameters start = createStartParameters(pT, phi, theta, charge);
0398
0399   ACTS_DEBUG(">>> Run navigation tests with:\n    pT = "
0400              << pT << "\n    phi = " << phi << "\n    theta = " << theta
0401              << "\n    charge = " << charge << "\n    index = " << index);
0402
0403   ACTS_DEBUG(">>> Test self consistency slpropagator");
0404   runSelfConsistencyTest(slpropagator, start, logger());
0405 }
0406
0407 BOOST_DATA_TEST_CASE(NavigatorEigenSelfConsistency,
0408                      eventGen ^ bdata::xrange(nTestsSelfConsistency), pT, phi,
0409                      theta, charge, index) {
0410   ACTS_LOCAL_LOGGER(getDefaultLogger("NavigatorTest", logLevel));
0411
0412   BoundTrackParameters start = createStartParameters(pT, phi, theta, charge);
0413
0414   ACTS_DEBUG(">>> Run navigation tests with:\n    pT = "
0415              << pT << "\n    phi = " << phi << "\n    theta = " << theta
0416              << "\n    charge = " << charge << "\n    index = " << index);
0417
0418   ACTS_DEBUG(">>> Test self consistency epropagator");
0419   runSelfConsistencyTest(epropagator, start, logger());
0420 }
0421
0422 BOOST_DATA_TEST_CASE(NavigatorRefStraightLineConsistency,
0423                      eventGen ^ bdata::xrange(nTestsRefConsistency), pT, phi,
0424                      theta, charge, index) {
0425   ACTS_LOCAL_LOGGER(getDefaultLogger("NavigatorTest", logLevel));
0426
0427   BoundTrackParameters start = createStartParameters(pT, phi, theta, charge);
0428
0429   ACTS_DEBUG(">>> Run navigation tests with:\n    pT = "
0430              << pT << "\n    phi = " << phi << "\n    theta = " << theta
0431              << "\n    charge = " << charge << "\n    index = " << index);
0432
0433   ACTS_DEBUG(">>> Test reference consistency slpropagator");
0434   runConsistencyTest(slpropagator, refslpropagator, start, logger());
0435 }
0436
0437 BOOST_DATA_TEST_CASE(NavigatorRefEigenConsistency,
0438                      eventGen ^ bdata::xrange(nTestsRefConsistency), pT, phi,
0439                      theta, charge, index) {
0440   ACTS_LOCAL_LOGGER(getDefaultLogger("NavigatorTest", logLevel));
0441
0442   BoundTrackParameters start = createStartParameters(pT, phi, theta, charge);
0443
0444   ACTS_DEBUG(">>> Run navigation tests with:\n    pT = "
0445              << pT << "\n    phi = " << phi << "\n    theta = " << theta
0446              << "\n    charge = " << charge << "\n    index = " << index);
0447
0448   ACTS_DEBUG(">>> Test reference consistency epropagator");
0449   runConsistencyTest(epropagator, refepropagator, start, logger());
0450 }
0451
0452 }  // namespace ActsTests