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