Core/Propagator/DirectNavigatorTests.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2019 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/data/test_case.hpp>
0010 #include <boost/test/tools/output_test_stream.hpp>
0011 #include <boost/test/unit_test.hpp>
0012
0013 #include "Acts/Definitions/Algebra.hpp"
0014 #include "Acts/Definitions/Direction.hpp"
0015 #include "Acts/Definitions/Units.hpp"
0016 #include "Acts/EventData/TrackParameters.hpp"
0017 #include "Acts/Geometry/GeometryContext.hpp"
0018 #include "Acts/MagneticField/ConstantBField.hpp"
0019 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0020 #include "Acts/Propagator/AbortList.hpp"
0021 #include "Acts/Propagator/ActionList.hpp"
0022 #include "Acts/Propagator/DirectNavigator.hpp"
0023 #include "Acts/Propagator/EigenStepper.hpp"
0024 #include "Acts/Propagator/MaterialInteractor.hpp"
0025 #include "Acts/Propagator/Navigator.hpp"
0026 #include "Acts/Propagator/Propagator.hpp"
0027 #include "Acts/Propagator/StandardAborters.hpp"
0028 #include "Acts/Propagator/SurfaceCollector.hpp"
0029 #include "Acts/Tests/CommonHelpers/CylindricalTrackingGeometry.hpp"
0030 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0031
0032 #include <algorithm>
0033 #include <array>
0034 #include <cmath>
0035 #include <iostream>
0036 #include <memory>
0037 #include <random>
0038 #include <tuple>
0039 #include <utility>
0040 #include <vector>
0041
0042 namespace Acts {
0043 class Surface;
0044 }  // namespace Acts
0045
0046 namespace bdata = boost::unit_test::data;
0047 using namespace Acts::UnitLiterals;
0048
0049 namespace Acts::Test {
0050
0051 // Create a test context
0052 GeometryContext tgContext = GeometryContext();
0053 MagneticFieldContext mfContext = MagneticFieldContext();
0054
0055 CylindricalTrackingGeometry cGeometry(tgContext);
0056 auto tGeometry = cGeometry();
0057
0058 // Create a navigator for this tracking geometry
0059 Navigator navigator({tGeometry});
0060 DirectNavigator dnavigator;
0061
0062 using BField = ConstantBField;
0063 using Stepper = EigenStepper<>;
0064 using ReferencePropagator = Propagator<Stepper, Navigator>;
0065 using DirectPropagator = Propagator<Stepper, DirectNavigator>;
0066
0067 const double Bz = 2_T;
0068 auto bField = std::make_shared<BField>(Vector3{0, 0, Bz});
0069 Stepper estepper(bField);
0070 Stepper dstepper(bField);
0071
0072 ReferencePropagator rpropagator(std::move(estepper), std::move(navigator));
0073 DirectPropagator dpropagator(std::move(dstepper), std::move(dnavigator));
0074
0075 const int ntests = 1000;
0076 const int skip = 0;
0077 bool referenceTiming = false;
0078 bool oversteppingTest = false;
0079 double oversteppingMaxStepSize = 1_mm;
0080
0081 /// The actual test nethod that runs the test
0082 /// can be used with several propagator types
0083 ///
0084 /// @tparam rpropagator_t is the reference propagator type
0085 /// @tparam dpropagator_t is the direct propagator type
0086 ///
0087 /// @param rprop is the reference propagator instance
0088 /// @param dprop is the direct propagator instance
0089 /// @param pT the transverse momentum
0090 /// @param phi the azimuthal angle of the track at creation
0091 /// @param theta the polar angle of the track at creation
0092 /// @param charge is the charge of the particle
0093 /// @param index is the run index from the test
0094 template <typename rpropagator_t, typename dpropagator_t>
0095 void runTest(const rpropagator_t& rprop, const dpropagator_t& dprop, double pT,
0096              double phi, double theta, int charge, int index) {
0097   double dcharge = -1 + 2 * charge;
0098
0099   if (index < skip) {
0100     return;
0101   }
0102
0103   // Define start parameters from ranom input
0104   double p = pT / sin(theta);
0105   CurvilinearTrackParameters start(Vector4(0, 0, 0, 0), phi, theta, dcharge / p,
0106                                    std::nullopt, ParticleHypothesis::pion());
0107
0108   using EndOfWorld = EndOfWorldReached;
0109
0110   // Action list and abort list
0111   using RefereceActionList = ActionList<MaterialInteractor, SurfaceCollector<>>;
0112   using ReferenceAbortList = AbortList<EndOfWorld>;
0113
0114   // Options definition
0115   using Options = PropagatorOptions<RefereceActionList, ReferenceAbortList>;
0116   Options pOptions(tgContext, mfContext);
0117   if (oversteppingTest) {
0118     pOptions.maxStepSize = oversteppingMaxStepSize;
0119   }
0120
0121   // Surface collector configuration
0122   auto& sCollector = pOptions.actionList.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 DirectActionList = ActionList<DirectNavigator::Initializer,
0145                                         MaterialInteractor, SurfaceCollector<>>;
0146
0147     // Direct options definition
0148     using DirectOptions = PropagatorOptions<DirectActionList, AbortList<>>;
0149     DirectOptions dOptions(tgContext, mfContext);
0150     // Set the surface sequence
0151     auto& dInitializer =
0152         dOptions.actionList.get<DirectNavigator::Initializer>();
0153     dInitializer.navSurfaces = surfaceSequence;
0154     // Surface collector configuration
0155     auto& dCollector = dOptions.actionList.template get<SurfaceCollector<>>();
0156     dCollector.selector.selectSensitive = true;
0157     dCollector.selector.selectMaterial = true;
0158
0159     // Now redo the propagation with the direct propagator
0160     const auto& ddResult =
0161         dprop.propagate(start, destination, dOptions).value();
0162     auto& ddSurfaces = ddResult.template get<SurfaceCollector<>::result_type>();
0163     auto& ddMaterial = ddResult.template get<MaterialInteractor::result_type>();
0164
0165     // CHECK if you have as many surfaces collected as the default navigator
0166     BOOST_CHECK_EQUAL(cSurfaces.collected.size(), ddSurfaces.collected.size());
0167     CHECK_CLOSE_REL(cMaterial.materialInX0, ddMaterial.materialInX0, 1e-3);
0168
0169     // Now redo the propagation with the direct propagator - without destination
0170     const auto& dwResult = dprop.propagate(start, dOptions).value();
0171     auto& dwSurfaces = dwResult.template get<SurfaceCollector<>::result_type>();
0172
0173     // CHECK if you have as many surfaces collected as the default navigator
0174     BOOST_CHECK_EQUAL(cSurfaces.collected.size(), dwSurfaces.collected.size());
0175   }
0176 }
0177
0178 // This test case checks that no segmentation fault appears
0179 // - this tests the collection of surfaces
0180 BOOST_DATA_TEST_CASE(
0181     test_direct_navigator,
0182     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 20,
0183                    bdata::distribution = std::uniform_real_distribution<double>(
0184                        0.15_GeV, 10_GeV))) ^
0185         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 21,
0186                        bdata::distribution =
0187                            std::uniform_real_distribution<double>(-M_PI,
0188                                                                   M_PI))) ^
0189         bdata::random(
0190             (bdata::engine = std::mt19937(), bdata::seed = 22,
0191              bdata::distribution =
0192                  std::uniform_real_distribution<double>(1.0, M_PI - 1.0))) ^
0193         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 23,
0194                        bdata::distribution =
0195                            std::uniform_int_distribution<std::uint8_t>(0, 1))) ^
0196         bdata::xrange(ntests),
0197     pT, phi, theta, charge, index) {
0198   // Run the test
0199   runTest(rpropagator, dpropagator, pT, phi, theta, charge, index);
0200 }
0201
0202 }  // namespace Acts::Test