Core/Propagator/KalmanExtrapolatorTests.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/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Direction.hpp"
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014 #include "Acts/Definitions/Units.hpp"
0015 #include "Acts/EventData/GenericCurvilinearTrackParameters.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/ActorList.hpp"
0021 #include "Acts/Propagator/EigenStepper.hpp"
0022 #include "Acts/Propagator/Navigator.hpp"
0023 #include "Acts/Propagator/Propagator.hpp"
0024 #include "Acts/Tests/CommonHelpers/CubicTrackingGeometry.hpp"
0025 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0026 #include "Acts/Utilities/Result.hpp"
0027
0028 #include <algorithm>
0029 #include <array>
0030 #include <map>
0031 #include <memory>
0032 #include <optional>
0033 #include <tuple>
0034 #include <utility>
0035 #include <vector>
0036
0037 namespace Acts {
0038 class Logger;
0039 struct EndOfWorldReached;
0040 }  // namespace Acts
0041
0042 using namespace Acts::UnitLiterals;
0043
0044 namespace Acts::Test {
0045
0046 using Jacobian = BoundMatrix;
0047 using Covariance = BoundSquareMatrix;
0048
0049 // Create a test context
0050 GeometryContext tgContext = GeometryContext();
0051 MagneticFieldContext mfContext = MagneticFieldContext();
0052 ///
0053 /// @brief the bound state propagation
0054 ///
0055 struct StepWiseActor {
0056   /// The result is the piece-wise jacobian
0057   struct this_result {
0058     std::vector<Jacobian> jacobians = {};
0059     std::vector<double> paths = {};
0060
0061     double fullPath = 0.;
0062
0063     bool finalized = false;
0064   };
0065   /// Broadcast the result type
0066   using result_type = this_result;
0067
0068   /// @brief Kalman sequence operation
0069   ///
0070   /// @tparam propagator_state_t is the type of Propagator state
0071   /// @tparam stepper_t Type of the stepper used for the propagation
0072   /// @tparam navigator_t Type of the navigator used for the propagation
0073   ///
0074   /// @param state is the mutable propagator state object
0075   /// @param stepper The stepper in use
0076   /// @param navigator The navigator in use
0077   /// @param result is the mutable result state object
0078   template <typename propagator_state_t, typename stepper_t,
0079             typename navigator_t>
0080   void act(propagator_state_t& state, const stepper_t& stepper,
0081            const navigator_t& navigator, result_type& result,
0082            const Logger& /*logger*/) const {
0083     // Listen to the surface and create bound state where necessary
0084     auto surface = navigator.currentSurface(state.navigation);
0085     if (surface && surface->associatedDetectorElement()) {
0086       // Create a bound state and log the jacobian
0087       auto boundState = stepper.boundState(state.stepping, *surface).value();
0088       result.jacobians.push_back(std::move(std::get<Jacobian>(boundState)));
0089       result.paths.push_back(std::get<double>(boundState));
0090     }
0091     // Also store the jacobian and full path
0092     if (state.stage == PropagatorStage::postPropagation && !result.finalized) {
0093       // Set the last stepping parameter
0094       result.paths.push_back(state.stepping.pathAccumulated);
0095       // Set the full parameter
0096       result.fullPath = state.stepping.pathAccumulated;
0097       // Remember that you finalized this
0098       result.finalized = true;
0099     }
0100   }
0101 };
0102
0103 ///
0104 /// @brief Unit test for Kalman fitter propagation
0105 ///
0106 BOOST_AUTO_TEST_CASE(kalman_extrapolator) {
0107   // Build detector, take the cubic detector
0108   CubicTrackingGeometry cGeometry(tgContext);
0109   auto detector = cGeometry();
0110
0111   // The Navigator through the detector geometry
0112   Navigator::Config cfg{detector};
0113   cfg.resolvePassive = false;
0114   cfg.resolveMaterial = true;
0115   cfg.resolveSensitive = true;
0116   Navigator navigator(cfg);
0117
0118   // Configure propagation with deactivated B-field
0119   auto bField = std::make_shared<ConstantBField>(Vector3(0., 0., 0.));
0120   using Stepper = EigenStepper<>;
0121   Stepper stepper(bField);
0122   using Propagator = Propagator<Stepper, Navigator>;
0123   Propagator propagator(stepper, navigator);
0124
0125   // Set initial parameters for the particle track
0126   Covariance cov;
0127   cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
0128       0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
0129       0, 0;
0130   // The start parameters
0131   CurvilinearTrackParameters start(Vector4(-3_m, 0, 0, 42_ns), 0_degree,
0132                                    90_degree, 1_e / 1_GeV, cov,
0133                                    ParticleHypothesis::pion());
0134
0135   // Create the ActionList and AbortList
0136   using StepWiseResult = StepWiseActor::result_type;
0137   using StepWiseActors = ActorList<StepWiseActor, EndOfWorldReached>;
0138
0139   // Create some options
0140   using StepWiseOptions = Propagator::Options<StepWiseActors>;
0141   StepWiseOptions swOptions(tgContext, mfContext);
0142
0143   using PlainActors = ActorList<EndOfWorldReached>;
0144   using PlainOptions = Propagator::Options<PlainActors>;
0145   PlainOptions pOptions(tgContext, mfContext);
0146
0147   // Run the standard propagation
0148   const auto& pResult = propagator.propagate(start, pOptions).value();
0149   // Let's get the end parameters and jacobian matrix
0150   const auto& pJacobian = *(pResult.transportJacobian);
0151
0152   // Run the stepwise propagation
0153   const auto& swResult = propagator.propagate(start, swOptions).value();
0154   auto swJacobianTest = swResult.template get<StepWiseResult>();
0155
0156   // (1) Path length test
0157   double accPath = 0.;
0158   auto swPaths = swJacobianTest.paths;
0159   // Sum up the step-wise paths, they are total though
0160   for (auto cpath = swPaths.rbegin(); cpath != swPaths.rend(); ++cpath) {
0161     if (cpath != swPaths.rend() - 1) {
0162       accPath += (*cpath) - (*(cpath + 1));
0163       continue;
0164     }
0165     accPath += (*cpath) - 0.;
0166   }
0167   CHECK_CLOSE_REL(swJacobianTest.fullPath, accPath, 1e-6);
0168
0169   // (2) Jacobian test
0170   Jacobian accJacobian = Jacobian::Identity();
0171   // The stepwise jacobians
0172   auto swJacobians = swJacobianTest.jacobians;
0173   // The last-step jacobian, needed for the full jacobian transport
0174   const auto& swlJacobian = *(swResult.transportJacobian);
0175
0176   // Build up the step-wise jacobians
0177   for (auto& j : swJacobians) {
0178     accJacobian = j * accJacobian;
0179   }
0180   accJacobian = swlJacobian * accJacobian;
0181   CHECK_CLOSE_OR_SMALL(pJacobian, accJacobian, 1e-6, 1e-9);
0182 }
0183
0184 }  // namespace Acts::Test