Core/Propagator/KalmanExtrapolatorTests.cpp

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