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