Core/Propagator/DirectNavigatorTests.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/MagneticField/MagneticFieldContext.hpp"
0018 #include "Acts/Propagator/ActorList.hpp"
0019 #include "Acts/Propagator/DirectNavigator.hpp"
0020 #include "Acts/Propagator/EigenStepper.hpp"
0021 #include "Acts/Propagator/MaterialInteractor.hpp"
0022 #include "Acts/Propagator/Navigator.hpp"
0023 #include "Acts/Propagator/Propagator.hpp"
0024 #include "Acts/Propagator/StandardAborters.hpp"
0025 #include "Acts/Propagator/StraightLineStepper.hpp"
0026 #include "Acts/Propagator/SurfaceCollector.hpp"
0027 #include "Acts/Surfaces/PlaneSurface.hpp"
0028 #include "Acts/Tests/CommonHelpers/CylindricalTrackingGeometry.hpp"
0029 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0030
0031 #include <cmath>
0032 #include <iostream>
0033 #include <memory>
0034 #include <numbers>
0035 #include <random>
0036 #include <utility>
0037 #include <vector>
0038
0039 namespace Acts {
0040 class Surface;
0041 }  // namespace Acts
0042
0043 namespace bdata = boost::unit_test::data;
0044 using namespace Acts::UnitLiterals;
0045
0046 namespace Acts::Test {
0047
0048 // Create a test context
0049 GeometryContext tgContext = GeometryContext();
0050 MagneticFieldContext mfContext = MagneticFieldContext();
0051
0052 CylindricalTrackingGeometry cGeometry(tgContext);
0053 auto tGeometry = cGeometry();
0054
0055 // Create a navigator for this tracking geometry
0056 Navigator navigator({tGeometry});
0057 DirectNavigator dnavigator;
0058
0059 using BField = ConstantBField;
0060 using Stepper = EigenStepper<>;
0061 using ReferencePropagator = Propagator<Stepper, Navigator>;
0062 using DirectPropagator = Propagator<Stepper, DirectNavigator>;
0063
0064 const double Bz = 2_T;
0065 auto bField = std::make_shared<BField>(Vector3{0, 0, Bz});
0066 Stepper estepper(bField);
0067 Stepper dstepper(bField);
0068
0069 ReferencePropagator rpropagator(std::move(estepper), std::move(navigator));
0070 DirectPropagator dpropagator(std::move(dstepper), std::move(dnavigator));
0071
0072 const int ntests = 1000;
0073 const int skip = 0;
0074 bool referenceTiming = false;
0075 bool oversteppingTest = false;
0076 double oversteppingMaxStepSize = 1_mm;
0077
0078 /// The actual test method that runs the test
0079 /// can be used with several propagator types
0080 ///
0081 /// @tparam rpropagator_t is the reference propagator type
0082 /// @tparam dpropagator_t is the direct propagator type
0083 ///
0084 /// @param rprop is the reference propagator instance
0085 /// @param dprop is the direct propagator instance
0086 /// @param pT the transverse momentum
0087 /// @param phi the azimuthal angle of the track at creation
0088 /// @param theta the polar angle of the track at creation
0089 /// @param charge is the charge of the particle
0090 /// @param index is the run index from the test
0091 template <typename rpropagator_t, typename dpropagator_t>
0092 void runTest(const rpropagator_t& rprop, const dpropagator_t& dprop, double pT,
0093              double phi, double theta, int charge, int index) {
0094   double dcharge = -1 + 2 * charge;
0095
0096   if (index < skip) {
0097     return;
0098   }
0099
0100   // Define start parameters from ranom input
0101   double p = pT / sin(theta);
0102   BoundTrackParameters start = BoundTrackParameters::createCurvilinear(
0103       Vector4::Zero(), phi, theta, dcharge / p, std::nullopt,
0104       ParticleHypothesis::pion());
0105
0106   using EndOfWorld = EndOfWorldReached;
0107
0108   // Action list and abort list
0109   using ReferenceActorList =
0110       ActorList<MaterialInteractor, SurfaceCollector<>, EndOfWorld>;
0111
0112   // Options definition
0113   using Options = typename rpropagator_t::template Options<ReferenceActorList>;
0114   Options pOptions(tgContext, mfContext);
0115   if (oversteppingTest) {
0116     pOptions.stepping.maxStepSize = oversteppingMaxStepSize;
0117   }
0118
0119   // Surface collector configuration
0120   auto& sCollector = pOptions.actorList.template get<SurfaceCollector<>>();
0121   sCollector.selector.selectSensitive = true;
0122   sCollector.selector.selectMaterial = true;
0123
0124   // Result is immediately used, non-valid result would indicate failure
0125   const auto& pResult = rprop.propagate(start, pOptions).value();
0126   auto& cSurfaces = pResult.template get<SurfaceCollector<>::result_type>();
0127   auto& cMaterial = pResult.template get<MaterialInteractor::result_type>();
0128   const Surface& destination = pResult.endParameters->referenceSurface();
0129
0130   std::cout << " - the standard navigator yielded "
0131             << cSurfaces.collected.size() << " collected surfaces" << std::endl;
0132
0133   if (!referenceTiming) {
0134     // Create the surface sequence
0135     std::vector<const Surface*> surfaceSequence;
0136     surfaceSequence.reserve(cSurfaces.collected.size());
0137     for (auto& cs : cSurfaces.collected) {
0138       surfaceSequence.push_back(cs.surface);
0139     }
0140
0141     // Action list for direct navigator with its initializer
0142     using DirectActorList = ActorList<MaterialInteractor, SurfaceCollector<>>;
0143
0144     // Direct options definition
0145     using DirectOptions =
0146         typename dpropagator_t::template Options<DirectActorList>;
0147     DirectOptions dOptions(tgContext, mfContext);
0148     // Set the surface sequence
0149     dOptions.navigation.surfaces = surfaceSequence;
0150     // Surface collector configuration
0151     auto& dCollector = dOptions.actorList.template get<SurfaceCollector<>>();
0152     dCollector.selector.selectSensitive = true;
0153     dCollector.selector.selectMaterial = true;
0154
0155     // Now redo the propagation with the direct propagator
0156     const auto& ddResult =
0157         dprop.propagate(start, destination, dOptions).value();
0158     auto& ddSurfaces = ddResult.template get<SurfaceCollector<>::result_type>();
0159     auto& ddMaterial = ddResult.template get<MaterialInteractor::result_type>();
0160
0161     // CHECK if you have as many surfaces collected as the default navigator
0162     BOOST_CHECK_EQUAL(cSurfaces.collected.size(), ddSurfaces.collected.size());
0163     CHECK_CLOSE_REL(cMaterial.materialInX0, ddMaterial.materialInX0, 1e-3);
0164
0165     // Now redo the propagation with the direct propagator - without destination
0166     const auto& dwResult = dprop.propagate(start, dOptions).value();
0167     auto& dwSurfaces = dwResult.template get<SurfaceCollector<>::result_type>();
0168
0169     // CHECK if you have as many surfaces collected as the default navigator
0170     BOOST_CHECK_EQUAL(cSurfaces.collected.size(), dwSurfaces.collected.size());
0171   }
0172 }
0173
0174 // This test case checks that no segmentation fault appears
0175 // - this tests the collection of surfaces
0176 BOOST_DATA_TEST_CASE(
0177     test_direct_navigator,
0178     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 20,
0179                    bdata::distribution = std::uniform_real_distribution<double>(
0180                        0.15_GeV, 10_GeV))) ^
0181         bdata::random(
0182             (bdata::engine = std::mt19937(), bdata::seed = 21,
0183              bdata::distribution = std::uniform_real_distribution<double>(
0184                  -std::numbers::pi, std::numbers::pi))) ^
0185         bdata::random(
0186             (bdata::engine = std::mt19937(), bdata::seed = 22,
0187              bdata::distribution = std::uniform_real_distribution<double>(
0188                  1., std::numbers::pi - 1.))) ^
0189         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 23,
0190                        bdata::distribution =
0191                            std::uniform_int_distribution<std::uint8_t>(0, 1))) ^
0192         bdata::xrange(ntests),
0193     pT, phi, theta, charge, index) {
0194   // Run the test
0195   runTest(rpropagator, dpropagator, pT, phi, theta, charge, index);
0196 }
0197
0198 struct NavigationBreakAborter {
0199   bool checkAbort(const auto& state, const auto& /*stepper*/, const auto& nav,
0200                   const auto& /*logger*/) const {
0201     return nav.navigationBreak(state.navigation);
0202   }
0203 };
0204
0205 // According to stackoverflow, clang before version 16 cannot handle
0206 // std::ranges::subrange See
0207 // https://stackoverflow.com/questions/64300832/why-does-clang-think-gccs-subrange-does-not-satisfy-gccs-ranges-begin-functi
0208 #if defined(__clang__) && __clang_major__ < 16
0209 #define CLANG_RANGE_BUG_WORKAROUND
0210 #endif
0211
0212 #ifdef CLANG_RANGE_BUG_WORKAROUND
0213 template <typename It>
0214 struct Subrange {
0215   It b, e;
0216   Subrange(It b_, It e_) : b(b_), e(e_) {}
0217   auto begin() const { return b; }
0218   auto end() const { return e; }
0219 };
0220 #endif
0221
0222 /// Run a simple test with a sequence of surfaces to check if fwd and backward
0223 /// navigation works
0224 #ifdef CLANG_RANGE_BUG_WORKAROUND
0225 template <typename ref_surfaces_t>
0226 #else
0227 template <std::ranges::range ref_surfaces_t>
0228 #endif
0229 void runSimpleTest(const std::vector<const Surface*>& surfaces,
0230                    Direction direction, const Surface* startSurface,
0231                    ref_surfaces_t expectedSurfaces) {
0232   Propagator<StraightLineStepper, DirectNavigator> prop(StraightLineStepper{},
0233                                                         DirectNavigator{});
0234
0235   using DirectActorList = ActorList<SurfaceCollector<>, NavigationBreakAborter>;
0236   using DirectOptions =
0237       typename Propagator<StraightLineStepper,
0238                           DirectNavigator>::template Options<DirectActorList>;
0239   DirectOptions pOptions(tgContext, mfContext);
0240   pOptions.direction = direction;
0241   pOptions.navigation.surfaces = surfaces;
0242   pOptions.navigation.startSurface = startSurface;
0243   auto& dCollector = pOptions.actorList.template get<SurfaceCollector<>>();
0244   dCollector.selector.selectSensitive = true;
0245   dCollector.selector.selectMaterial = true;
0246   dCollector.selector.selectPassive = true;
0247
0248   // Create the start parameters in the middle of the start surface
0249   BoundTrackParameters startParameters = BoundTrackParameters(
0250       startSurface->getSharedPtr(),
0251       {0.0_mm, 0.0_mm, 0.0_rad, 0.0_rad, 1.0 / 1.0_GeV, 0.0_ns}, std::nullopt,
0252       ParticleHypothesis::muon());
0253
0254   // Propagate the track
0255   auto result = prop.propagate(startParameters, pOptions);
0256
0257   // Check if the result is valid
0258   BOOST_REQUIRE(result.ok());
0259
0260   // Check if the surfaces are the same
0261   const auto& collectedSurfaceHits =
0262       result->get<SurfaceCollector<>::result_type>().collected;
0263   std::vector<const Surface*> collectedSurfaces;
0264   std::ranges::transform(collectedSurfaceHits,
0265                          std::back_inserter(collectedSurfaces),
0266                          [](const auto& hit) { return hit.surface; });
0267   // the initial surface is twice in the collection
0268   collectedSurfaces.erase(
0269       std::unique(collectedSurfaces.begin(), collectedSurfaces.end()),
0270       collectedSurfaces.end());
0271   BOOST_CHECK_EQUAL_COLLECTIONS(
0272       collectedSurfaces.begin(), collectedSurfaces.end(),
0273       expectedSurfaces.begin(), expectedSurfaces.end());
0274 }
0275
0276 BOOST_AUTO_TEST_CASE(test_direct_navigator_fwd_bwd) {
0277   // Create 10 surfaces at z = 0, 100, 200, ..., 900
0278   std::vector<std::shared_ptr<const Acts::Surface>> surfaces;
0279   for (int i = 0; i < 10; i++) {
0280     Transform3 transform = Transform3::Identity();
0281     transform.translate(Vector3{0.0_mm, 0.0_mm, i * 100.0_mm});
0282     auto surface = Surface::makeShared<PlaneSurface>(transform, nullptr);
0283     surface->assignGeometryId(
0284         Acts::GeometryIdentifier().withVolume(1).withLayer(1).withSensitive(i +
0285                                                                             1));
0286     surfaces.push_back(surface);
0287   }
0288
0289   // Create vector of pointers to the surfaces
0290   std::vector<const Acts::Surface*> surfacePointers;
0291   std::ranges::transform(surfaces, std::back_inserter(surfacePointers),
0292                          [](const auto& s) { return s.get(); });
0293
0294   for (auto it = surfacePointers.begin(); it != surfacePointers.end(); ++it) {
0295     runSimpleTest(surfacePointers, Direction::Forward(), *it,
0296 #ifndef CLANG_RANGE_BUG_WORKAROUND
0297                   std::ranges::subrange{it, surfacePointers.end()});
0298 #else
0299                   Subrange{it, surfacePointers.end()});
0300 #endif
0301   }
0302   for (auto it = surfacePointers.rbegin(); it != surfacePointers.rend(); ++it) {
0303     runSimpleTest(surfacePointers, Direction::Backward(), *it,
0304 #ifndef CLANG_RANGE_BUG_WORKAROUND
0305                   std::ranges::subrange{it, surfacePointers.rend()});
0306 #else
0307                   Subrange{it, surfacePointers.rend()});
0308 #endif
0309   }
0310 }
0311
0312 }  // namespace Acts::Test