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