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