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