From e9781f4cc05dfac309f2ed1593ef75c2aaf6d08f Mon Sep 17 00:00:00 2001 From: Andreas Salzburger Date: Wed, 11 Oct 2017 14:56:36 +0200 Subject: [PATCH 01/11] Introduction of new bField, implementation of field cell caching and cache_call signature for the propagation --- .../ACTS/MagneticField/SharedBFieldMap.hpp | 75 +++++++++++++++ Core/include/ACTS/Propagator/AtlasStepper.hpp | 31 ++++++- Core/include/ACTS/Propagator/EigenStepper.hpp | 51 ++++++---- Core/include/ACTS/Propagator/Propagator.hpp | 92 +++++++++++++------ 4 files changed, 200 insertions(+), 49 deletions(-) create mode 100644 Core/include/ACTS/MagneticField/SharedBFieldMap.hpp diff --git a/Core/include/ACTS/MagneticField/SharedBFieldMap.hpp b/Core/include/ACTS/MagneticField/SharedBFieldMap.hpp new file mode 100644 index 00000000000..25aa8541864 --- /dev/null +++ b/Core/include/ACTS/MagneticField/SharedBFieldMap.hpp @@ -0,0 +1,75 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 ACTS project team +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef ACTS_MAGNETICFIELD_SHAREDBFIELD_H +#define ACTS_MAGNETICFIELD_SHAREDBFIELD_H + +#include "ACTS/MagneticField/concept/AnyFieldLookup.hpp" +#include "ACTS/Utilities/Definitions.hpp" + +namespace Acts { + +/// @ingroup MagneticField +/// +/// @brief allows to use a shared magnetic field +/// in several places and with multiple steppers +/// mainly targeted to save memory +template class SharedBField { + public: + /// @brief the constructur with a shared pointer + /// @tparam bField is the shared BField to be stored + SharedBField(std::shared_ptr bField) : + m_bField(bField) + {} + + /// @brief retrieve magnetic field value + /// + /// @param [in] position global 3D position + /// + /// @return magnetic field vector at given position + Vector3D + getField(const Vector3D& position) const + { + return m_bField->getField(position); + } + + /// @brief retrieve field cell for given position + /// + /// @param [in] position global 3D position + /// @return field cell containing the given global position + /// + /// @pre The given @c position must lie within the range of the underlying + /// magnetic field map. + concept::AnyFieldCell<> + getFieldCell(const Vector3D& position) const + { + return m_bField->getFieldCell(position); + } + + /// @brief retrieve magnetic field value & its gradient + /// + /// @param [in] position global 3D position + /// @param [out] derivative gradient of magnetic field vector as (3x3) matrix + /// @return magnetic field vector + /// + /// @note currently the derivative is not calculated + /// @todo return derivative + Vector3D + getFieldGradient(const Vector3D& position, + ActsMatrixD<3, 3>& derivative) const + { + return m_bField->getField(position); + } + private: + std::shared_ptr m_bField; + +}; + +} // namespace Acts + +#endif // ACTS_MAGNETICFIELD_SHAREDBFIELD_H diff --git a/Core/include/ACTS/Propagator/AtlasStepper.hpp b/Core/include/ACTS/Propagator/AtlasStepper.hpp index 5489af64c9e..823c8ad8ceb 100644 --- a/Core/include/ACTS/Propagator/AtlasStepper.hpp +++ b/Core/include/ACTS/Propagator/AtlasStepper.hpp @@ -13,6 +13,7 @@ #include "ACTS/EventData/TrackParameters.hpp" #include "ACTS/Surfaces/Surface.hpp" #include "ACTS/Utilities/Units.hpp" +#include "ACTS/MagneticField/concept/AnyFieldLookup.hpp" namespace Acts { @@ -36,6 +37,12 @@ class AtlasStepper double parameters[NGlobalPars] = {0., 0., 0., 0., 0.}; const ActsSymMatrixD* covariance; double jacobian[NGlobalPars * NGlobalPars]; + + /// Lazily initialized cache + /// It caches the current magneticl field cell and stays interpolates within + /// as long as this is valid. See step() code for details. + bool field_cache_ready = false; + concept::AnyFieldCell<> field_cache; Vector3D position() const @@ -446,6 +453,24 @@ class AtlasStepper return i.pathLength; } + /// Get the field for the stepping + /// It checks first if the access is still within the Cell, + /// and updates the cell if necessary, then it takes the field + /// from the cell + /// @param [in,out] cache is the propagation cache associated with the track + /// the magnetic field cell is used (and potentially updated) + /// @param [in] pos is the field position + Vector3D getField(Cache& cache, const Vector3D& pos) const + { + if (!cache.field_cache_ready || !cache.field_cache.isInside(pos)){ + cache.field_cache = m_bField.getFieldCell(pos); + } + // get the field from the cell + cache.field = cache.field_cache.getField(pos); + return cache.field; + } + + double step(Cache& cache, double& h) const { @@ -463,7 +488,7 @@ class AtlasStepper // if new field is required get it if (cache.newfield) { const Vector3D pos(R[0], R[1], R[2]); - f0 = m_bField.getField(pos); + f0 = getField(cache, pos); } else { f0 = cache.field; } @@ -491,7 +516,7 @@ class AtlasStepper // if (!Helix) { const Vector3D pos(R[0] + A1 * S4, R[1] + B1 * S4, R[2] + C1 * S4); - f = m_bField.getField(pos); + f = getField(cache, pos); } else { f = f0; } @@ -511,7 +536,7 @@ class AtlasStepper // if (!Helix) { const Vector3D pos(R[0] + h * A4, R[1] + h * B4, R[2] + h * C4); - f = m_bField.getField(pos); + f = getField(cache, pos); } else { f = f0; } diff --git a/Core/include/ACTS/Propagator/EigenStepper.hpp b/Core/include/ACTS/Propagator/EigenStepper.hpp index 0db363ab148..bb11116af33 100644 --- a/Core/include/ACTS/Propagator/EigenStepper.hpp +++ b/Core/include/ACTS/Propagator/EigenStepper.hpp @@ -14,6 +14,7 @@ #include "ACTS/Surfaces/Surface.hpp" #include "ACTS/Utilities/Definitions.hpp" #include "ACTS/Utilities/Units.hpp" +#include "ACTS/MagneticField/concept/AnyFieldLookup.hpp" namespace Acts { @@ -54,10 +55,10 @@ class EigenStepper , qop(par.charge() / par.momentum().norm()) , cov(par.covariance()) { - const double phi = dir.phi(); - const double theta = dir.theta(); - if (cov) { + const double phi = dir.phi(); + const double theta = dir.theta(); + const auto transform = par.referenceSurface().transform().matrix(); jacobian(0, eLOC_0) = transform(0, eLOC_0); jacobian(0, eLOC_1) = transform(0, eLOC_1); @@ -100,16 +101,18 @@ class EigenStepper /// Jacobian used to transport the covariance matrix ActsMatrixD<7, 5> jacobian = ActsMatrixD<7, 5>::Zero(); - /// ??? + /// The propagation derivative ActsVectorD<7> derivative = ActsVectorD<7>::Zero(); /// Covariance matrix assocated with the initial error on track parameters const ActsSymMatrixD<5>* cov = nullptr; - /// Lazily initialized cache which contains an approximation of the - /// magnetic field at the current position. See step() code for details. + /// Lazily initialized cache + /// It caches the current magneticl field cell and stays interpolates within + /// as long as this is valid. See step() code for details. bool field_cache_ready = false; - Vector3D field_cache = Vector3D(0, 0, 0); + concept::AnyFieldCell<> field_cache; + }; // This struct is a meta-function which normally maps to BoundParameters... @@ -223,6 +226,8 @@ class EigenStepper } /// Convert the propagation cache to track parameters at a certain surface + /// @param [in] cache Propagation cache used + /// @param [in] surface Destination surface to which the conversion is done template static BoundParameters convert(Cache& cache, const S& surface) @@ -281,6 +286,22 @@ class EigenStepper return i.pathLength; } + /// Get the field for the stepping + /// It checks first if the access is still within the Cell, + /// and updates the cell if necessary, then it takes the field + /// from the cell + /// @param [in,out] cache is the propagation cache associated with the track + /// the magnetic field cell is used (and potentially updated) + /// @param [in] pos is the field position + Vector3D getField(Cache& cache, const Vector3D& pos) const + { + if (!cache.field_cache_ready || !cache.field_cache.isInside(pos)){ + cache.field_cache = m_bField.getFieldCell(pos); + } + // get the field from the cell + return std::move(cache.field_cache.getField(pos)); + } + /// Perform a Runge-Kutta track parameter propagation step /// /// @param[in,out] cache is the propagation cache associated with the track @@ -301,14 +322,8 @@ class EigenStepper double h2, half_h; Vector3D B_middle, B_last, k2, k3, k4; - // Initialize the magnetic field cache on the first run - if (!cache.field_cache_ready) { - cache.field_cache = m_bField.getField(cache.pos); - cache.field_cache_ready = true; - } - // First Runge-Kutta point (at current position) - const Vector3D B_first = cache.field_cache; + const Vector3D B_first = getField(cache, cache.pos); const Vector3D k1 = qop * cache.dir.cross(B_first); // The following functor starts to perform a Runge-Kutta step of a certain @@ -322,7 +337,7 @@ class EigenStepper // Second Runge-Kutta point const Vector3D pos1 = cache.pos + half_h * cache.dir + h2 / 8 * k1; - B_middle = m_bField.getField(pos1); + B_middle = getField(cache, pos1);; k2 = qop * (cache.dir + half_h * k1).cross(B_middle); // Third Runge-Kutta point @@ -330,7 +345,7 @@ class EigenStepper // Last Runge-Kutta point const Vector3D pos2 = cache.pos + h * cache.dir + h2 / 2 * k3; - B_last = m_bField.getField(pos2); + B_last = getField(cache, pos2); k4 = qop * (cache.dir + h * k3).cross(B_last); // Return an estimate of the local integration error @@ -418,10 +433,6 @@ class EigenStepper cache.derivative.template head<3>() = cache.dir; cache.derivative.template segment<3>(3) = k4; - // We approximate the magnetic field at our new position as being equal to - // the last field measurement, thusly avoiding one field lookup per step. - cache.field_cache = B_last; - // Return the updated step size return h; } diff --git a/Core/include/ACTS/Propagator/Propagator.hpp b/Core/include/ACTS/Propagator/Propagator.hpp index fadffb1d8a8..eef0e4f7f23 100644 --- a/Core/include/ACTS/Propagator/Propagator.hpp +++ b/Core/include/ACTS/Propagator/Propagator.hpp @@ -91,6 +91,10 @@ namespace propagation { class Propagator final { public: + + // Type of cache object used by the propagation implementation + typedef typename Impl::template cache_type cache_type; + /// @brief Options for propagate() call /// /// @tparam ObserverList List of observer types called after each @@ -111,7 +115,7 @@ namespace propagation { /// Maximum number of steps for one propagate() call unsigned int max_steps = 1000; - /// Required distance to surface + /// Required distance to surface double target_surface_distance = 1 * units::_um; /// Absolute minimum step size @@ -167,9 +171,10 @@ namespace propagation { /// template using obs_list_result_t = typename result_type_helper::type; - + public: - /// @brief Propagate track parameters + + /// @brief Propagate track parameters - User method /// /// This function performs the propagation of the track parameters using the /// internal implementation object, until at least one abort condition is @@ -195,8 +200,6 @@ namespace propagation { propagate(const TrackParameters& start, const Options& options) const { - // Type of internal cache object used by the propagation implementation - typedef typename Impl::template cache_type cache_type; // Type of track parameters produced by the propagation typedef typename Impl::template return_parameter_type @@ -210,7 +213,7 @@ namespace propagation { "return track parameter type must be copy-constructible"); // Initialize the propagation result object - result_type r(Status::INPROGRESS); + result_type r(Status::INPROGRESS); // Compute the signed path limit and maximum step size const double signed_pathLimit @@ -218,21 +221,21 @@ namespace propagation { double stepMax = options.direction * options.max_step_size; // Initialize the internal propagation cache - cache_type propagation_cache(start); + cache_type cache(start); // Propagation loop for (; r.steps < options.max_steps; ++r.steps) { // Perform a propagation step (which can alter the step size) - r.pathLength += m_impl.step(propagation_cache, stepMax); + r.pathLength += m_impl.step(cache, stepMax); // Call the observers with the current and previous track parameters, // and let them fill in some propagation results - options.observer_list(propagation_cache, r); + options.observer_list(cache, r); // Is it time to stop the propagation? if (std::abs(r.pathLength) >= options.max_path_length - || options.stop_conditions(r, propagation_cache, stepMax)) { + || options.stop_conditions(r, cache, stepMax)) { // Compute the final results and mark the propagation as successful r.endParameters = std::make_unique( - m_impl.convert(propagation_cache)); + m_impl.convert(cache)); r.status = Status::SUCCESS; break; } @@ -240,22 +243,26 @@ namespace propagation { if (std::abs(stepMax) > std::abs(signed_pathLimit - r.pathLength)) stepMax = signed_pathLimit - r.pathLength; } - return r; + return r; } - /// @brief Propagate track parameters + /// @brief Propagate track parameters - Expert method with propagation cache /// /// This function performs the propagation of the track parameters according /// to the internal implementation object until at least one abort condition /// is fulfilled, the destination surface is hit or the maximum number of /// steps/path length as given in the propagation options is reached. /// + /// It does check/re-use the propgation cache as much as possible for + /// performance reasons + /// /// @tparam TrackParameters Type of initial track parameters to propagate /// @tparam Surface Type of target surface /// @tparam ObserverList Type list of observers /// @tparam AbortList Type list of abort conditions /// - /// @param [in] start Initial track parameters to propagate + /// @param [in] cache Stepper cache built/updated from the start parameters + /// @param [in] target Target surface of to propagate to /// @param [in] options Propagation options /// /// @return Propagation result containing the propagation status, final @@ -268,14 +275,12 @@ namespace propagation { obs_list_result_t< typename Impl::template return_parameter_type, ObserverList> - propagate(const TrackParameters& start, - const Surface& target, - const Options& options) const + propagate_with_cache(cache_type& cache, + const TrackParameters& start, + const Surface& target, + const Options& options) const { - // Type of internal cache object used by the propagation implementation - typedef typename Impl::template cache_type - cache_type; - + // Type of track parameters produced at the end of the propagation typedef typename Impl::template return_parameter_type @@ -289,16 +294,13 @@ namespace propagation { "return track parameter type must be copy-constructible"); // Initialize the propagation result object - result_type r(Status::INPROGRESS); + result_type r(Status::INPROGRESS); // Compute the signed path limit and maximum step size const double signed_pathLimit = options.direction * options.max_path_length; double stepMax = options.direction * options.max_step_size; - // Initialize the internal propagation cache - cache_type cache(start); - // Compute the distance to the target surface double distance = m_impl.distance(target, cache.position(), cache.direction()); @@ -325,7 +327,7 @@ namespace propagation { distance = m_impl.distance(target, cache.position(), cache.direction()); // Is it time to stop the propagation ? - if (std::abs(distance) < options.target_surface_distance + if (std::abs(distance) < options.target_surface_distance || std::abs(r.pathLength) >= options.max_path_length || options.stop_conditions(r, cache, stepMax)) { // Compute the final results and mark the propagation as successful @@ -343,6 +345,44 @@ namespace propagation { return r; } + /// @brief Propagate track parameters - User method + /// + /// This function performs the propagation of the track parameters according + /// to the internal implementation object until at least one abort condition + /// is fulfilled, the destination surface is hit or the maximum number of + /// steps/path length as given in the propagation options is reached. + /// + /// A stepper cache object is built internally for this call and the + /// Expert method with the cache call signature is called. + /// + /// @tparam TrackParameters Type of initial track parameters to propagate + /// @tparam Surface Type of target surface + /// @tparam ObserverList Type list of observers + /// @tparam AbortList Type list of abort conditions + /// + /// @param [in] start Initial track parameters to propagate + /// @param [in] target Target surface of to propagate to + /// @param [in] options Propagation options + /// + /// @return Propagation result containing the propagation status, final + /// track parameters, and output of observers (if they produce any) + /// + template + obs_list_result_t< + typename Impl::template return_parameter_type, + ObserverList> + propagate(const TrackParameters& start, + const Surface& target, + const Options& options) const + { + // Initialize the internal propagation cache + cache_type cache(start); + return propagate_with_cache(cache, start, target, options); + } + private: /// implementation of propagation algorithm Impl m_impl; From 7f756af1b6ecf771a3f88849055f059cbc0235e9 Mon Sep 17 00:00:00 2001 From: Andreas Salzburger Date: Wed, 11 Oct 2017 15:01:05 +0200 Subject: [PATCH 02/11] Documentation updates and todo adding --- Core/include/ACTS/Propagator/EigenStepper.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Core/include/ACTS/Propagator/EigenStepper.hpp b/Core/include/ACTS/Propagator/EigenStepper.hpp index bb11116af33..eacc97d7d0b 100644 --- a/Core/include/ACTS/Propagator/EigenStepper.hpp +++ b/Core/include/ACTS/Propagator/EigenStepper.hpp @@ -48,6 +48,7 @@ class EigenStepper struct Cache { /// Constructor from the initial track parameters + /// @tparam [in] par are the track parameters template explicit Cache(const T& par) : pos(par.position()) @@ -162,6 +163,8 @@ class EigenStepper EigenStepper(BField bField = BField()) : m_bField(std::move(bField)){}; /// Convert the propagation cache to curvilinear parameters + /// @param cache is the stepper cache + /// @todo check: what if cache is already in courvilinear ? is this caught ? static CurvilinearParameters convert(const Cache& cache) { From 53c2daf3124f9f2a3401d27ae5ae4e6bb215e428 Mon Sep 17 00:00:00 2001 From: Andreas Salzburger Date: Wed, 11 Oct 2017 15:04:01 +0200 Subject: [PATCH 03/11] More documentation update - this fixes ACTS-385 and ACTS-379 --- Core/include/ACTS/Propagator/AtlasStepper.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Core/include/ACTS/Propagator/AtlasStepper.hpp b/Core/include/ACTS/Propagator/AtlasStepper.hpp index 823c8ad8ceb..cc558d523cc 100644 --- a/Core/include/ACTS/Propagator/AtlasStepper.hpp +++ b/Core/include/ACTS/Propagator/AtlasStepper.hpp @@ -294,6 +294,9 @@ class AtlasStepper return CurvilinearParameters(std::move(cov), gp, mom, charge); } + /// convert method into bound parameters + /// @param cache is the propagation cache to be converted + /// @param s the target surface static BoundParameters convert(Cache& cache, const Surface& s) { From 6c1d94efcd84879abf516a09b536d6c4ac49c32a Mon Sep 17 00:00:00 2001 From: Andreas Salzburger Date: Thu, 12 Oct 2017 13:51:36 +0200 Subject: [PATCH 04/11] set the cache flag to true once it is prepared --- Core/include/ACTS/Propagator/AtlasStepper.hpp | 1 + Core/include/ACTS/Propagator/EigenStepper.hpp | 1 + 2 files changed, 2 insertions(+) diff --git a/Core/include/ACTS/Propagator/AtlasStepper.hpp b/Core/include/ACTS/Propagator/AtlasStepper.hpp index cc558d523cc..60817423660 100644 --- a/Core/include/ACTS/Propagator/AtlasStepper.hpp +++ b/Core/include/ACTS/Propagator/AtlasStepper.hpp @@ -466,6 +466,7 @@ class AtlasStepper Vector3D getField(Cache& cache, const Vector3D& pos) const { if (!cache.field_cache_ready || !cache.field_cache.isInside(pos)){ + cache.field_cache_ready = true; cache.field_cache = m_bField.getFieldCell(pos); } // get the field from the cell diff --git a/Core/include/ACTS/Propagator/EigenStepper.hpp b/Core/include/ACTS/Propagator/EigenStepper.hpp index eacc97d7d0b..21a14161530 100644 --- a/Core/include/ACTS/Propagator/EigenStepper.hpp +++ b/Core/include/ACTS/Propagator/EigenStepper.hpp @@ -299,6 +299,7 @@ class EigenStepper Vector3D getField(Cache& cache, const Vector3D& pos) const { if (!cache.field_cache_ready || !cache.field_cache.isInside(pos)){ + cache.field_cache_ready = true; cache.field_cache = m_bField.getFieldCell(pos); } // get the field from the cell From db39cee201a0f85676eeb4c9c2354605dc0af44d Mon Sep 17 00:00:00 2001 From: Andreas Salzburger Date: Fri, 13 Oct 2017 12:11:54 +0200 Subject: [PATCH 05/11] copy covariance (not just a plain pointer caching that could run out of scope) explicitly calls the method propagate_with_cache in order not to make clear what it does prepare interface for internal cache_updating (to be implemented) fix integration test Fixes ACTS-385 - further updates will be done in a new branch/MR --- Core/include/ACTS/Propagator/EigenStepper.hpp | 25 ++++++++++++------- .../covariance_validation_fixture.hpp | 2 ++ 2 files changed, 18 insertions(+), 9 deletions(-) diff --git a/Core/include/ACTS/Propagator/EigenStepper.hpp b/Core/include/ACTS/Propagator/EigenStepper.hpp index 21a14161530..1008ca490b7 100644 --- a/Core/include/ACTS/Propagator/EigenStepper.hpp +++ b/Core/include/ACTS/Propagator/EigenStepper.hpp @@ -49,14 +49,19 @@ class EigenStepper { /// Constructor from the initial track parameters /// @tparam [in] par are the track parameters + /// + /// @note the covariance matrix is copied when needed template explicit Cache(const T& par) : pos(par.position()) , dir(par.momentum().normalized()) , qop(par.charge() / par.momentum().norm()) - , cov(par.covariance()) + , cov_transport(false) { - if (cov) { + if (par.covariance()) { + cov_transport = true; + cov = ActsSymMatrixD<5>(*par.covariance()); + const double phi = dir.phi(); const double theta = dir.theta(); @@ -105,8 +110,10 @@ class EigenStepper /// The propagation derivative ActsVectorD<7> derivative = ActsVectorD<7>::Zero(); - /// Covariance matrix assocated with the initial error on track parameters - const ActsSymMatrixD<5>* cov = nullptr; + /// Covariance matrix (and indicator) + //// assocated with the initial error on track parameters + bool cov_transport = false; + ActsSymMatrixD<5> cov = ActsSymMatrixD<5>::Zero(); /// Lazily initialized cache /// It caches the current magneticl field cell and stays interpolates within @@ -172,7 +179,7 @@ class EigenStepper std::unique_ptr> cov = nullptr; // Perform error propagation if an initial covariance matrix was provided - if (cache.cov) { + if (cache.cov_transport) { const double phi = cache.dir.phi(); const double theta = cache.dir.theta(); @@ -220,7 +227,7 @@ class EigenStepper jacobian -= tmp; auto jac = J * jacobian; - cov = std::make_unique>(jac * (*cache.cov) + cov = std::make_unique>(jac * cache.cov * jac.transpose()); } @@ -239,7 +246,7 @@ class EigenStepper std::unique_ptr> cov = nullptr; // Perform error propagation if an initial covariance matrix was provided - if (cache.cov) { + if (cache.cov_transport) { const double phi = cache.dir.phi(); const double theta = cache.dir.theta(); @@ -270,7 +277,7 @@ class EigenStepper jacobian -= tmp; auto jac = J * cache.jacobian; - cov = std::make_unique>(jac * (*cache.cov) + cov = std::make_unique>(jac * cache.cov * jac.transpose()); } @@ -364,7 +371,7 @@ class EigenStepper } // When doing error propagation, update the associated Jacobian matrix - if (cache.cov) { + if (cache.cov_transport) { ActsMatrixD<7, 7> D = ActsMatrixD<7, 7>::Zero(); D(6, 6) = 1; diff --git a/IntegrationTests/covariance_validation_fixture.hpp b/IntegrationTests/covariance_validation_fixture.hpp index 1eb848c105b..07b4f10f152 100644 --- a/IntegrationTests/covariance_validation_fixture.hpp +++ b/IntegrationTests/covariance_validation_fixture.hpp @@ -26,6 +26,8 @@ namespace IntegrationTest { { } + /// Numerical transport of covariance using the ridder's algorithm + /// this is for covariance propagation validation template ActsSymMatrixD<5> calculateCovariance(const CurvilinearParameters& trackPars, From 2fecc1e126fa64d47ca558029e398563fe4c8c7a Mon Sep 17 00:00:00 2001 From: Andreas Salzburger Date: Fri, 13 Oct 2017 14:22:03 +0200 Subject: [PATCH 06/11] applying clang-format --- .../ACTS/MagneticField/SharedBFieldMap.hpp | 100 +++++++++--------- Core/include/ACTS/Propagator/AtlasStepper.hpp | 18 ++-- Core/include/ACTS/Propagator/EigenStepper.hpp | 39 +++---- Core/include/ACTS/Propagator/Propagator.hpp | 32 +++--- 4 files changed, 94 insertions(+), 95 deletions(-) diff --git a/Core/include/ACTS/MagneticField/SharedBFieldMap.hpp b/Core/include/ACTS/MagneticField/SharedBFieldMap.hpp index 25aa8541864..de00c868897 100644 --- a/Core/include/ACTS/MagneticField/SharedBFieldMap.hpp +++ b/Core/include/ACTS/MagneticField/SharedBFieldMap.hpp @@ -18,58 +18,58 @@ namespace Acts { /// /// @brief allows to use a shared magnetic field /// in several places and with multiple steppers -/// mainly targeted to save memory -template class SharedBField { - public: - /// @brief the constructur with a shared pointer - /// @tparam bField is the shared BField to be stored - SharedBField(std::shared_ptr bField) : - m_bField(bField) - {} - - /// @brief retrieve magnetic field value - /// - /// @param [in] position global 3D position - /// - /// @return magnetic field vector at given position - Vector3D - getField(const Vector3D& position) const - { - return m_bField->getField(position); - } +/// mainly targeted to save memory +template +class SharedBField +{ +public: + /// @brief the constructur with a shared pointer + /// @tparam bField is the shared BField to be stored + SharedBField(std::shared_ptr bField) : m_bField(bField) {} - /// @brief retrieve field cell for given position - /// - /// @param [in] position global 3D position - /// @return field cell containing the given global position - /// - /// @pre The given @c position must lie within the range of the underlying - /// magnetic field map. - concept::AnyFieldCell<> - getFieldCell(const Vector3D& position) const - { - return m_bField->getFieldCell(position); - } + /// @brief retrieve magnetic field value + /// + /// @param [in] position global 3D position + /// + /// @return magnetic field vector at given position + Vector3D + getField(const Vector3D& position) const + { + return m_bField->getField(position); + } - /// @brief retrieve magnetic field value & its gradient - /// - /// @param [in] position global 3D position - /// @param [out] derivative gradient of magnetic field vector as (3x3) matrix - /// @return magnetic field vector - /// - /// @note currently the derivative is not calculated - /// @todo return derivative - Vector3D - getFieldGradient(const Vector3D& position, - ActsMatrixD<3, 3>& derivative) const - { - return m_bField->getField(position); - } - private: - std::shared_ptr m_bField; - + /// @brief retrieve field cell for given position + /// + /// @param [in] position global 3D position + /// @return field cell containing the given global position + /// + /// @pre The given @c position must lie within the range of the underlying + /// magnetic field map. + concept::AnyFieldCell<> + getFieldCell(const Vector3D& position) const + { + return m_bField->getFieldCell(position); + } + + /// @brief retrieve magnetic field value & its gradient + /// + /// @param [in] position global 3D position + /// @param [out] derivative gradient of magnetic field vector as (3x3) matrix + /// @return magnetic field vector + /// + /// @note currently the derivative is not calculated + /// @todo return derivative + Vector3D + getFieldGradient(const Vector3D& position, + ActsMatrixD<3, 3>& derivative) const + { + return m_bField->getField(position); + } + +private: + std::shared_ptr m_bField; }; -} // namespace Acts +} // namespace Acts -#endif // ACTS_MAGNETICFIELD_SHAREDBFIELD_H +#endif // ACTS_MAGNETICFIELD_SHAREDBFIELD_H diff --git a/Core/include/ACTS/Propagator/AtlasStepper.hpp b/Core/include/ACTS/Propagator/AtlasStepper.hpp index 60817423660..339703c25cf 100644 --- a/Core/include/ACTS/Propagator/AtlasStepper.hpp +++ b/Core/include/ACTS/Propagator/AtlasStepper.hpp @@ -11,9 +11,9 @@ #include #include "ACTS/EventData/TrackParameters.hpp" +#include "ACTS/MagneticField/concept/AnyFieldLookup.hpp" #include "ACTS/Surfaces/Surface.hpp" #include "ACTS/Utilities/Units.hpp" -#include "ACTS/MagneticField/concept/AnyFieldLookup.hpp" namespace Acts { @@ -37,11 +37,11 @@ class AtlasStepper double parameters[NGlobalPars] = {0., 0., 0., 0., 0.}; const ActsSymMatrixD* covariance; double jacobian[NGlobalPars * NGlobalPars]; - - /// Lazily initialized cache - /// It caches the current magneticl field cell and stays interpolates within + + /// Lazily initialized cache + /// It caches the current magneticl field cell and stays interpolates within /// as long as this is valid. See step() code for details. - bool field_cache_ready = false; + bool field_cache_ready = false; concept::AnyFieldCell<> field_cache; Vector3D @@ -463,18 +463,18 @@ class AtlasStepper /// @param [in,out] cache is the propagation cache associated with the track /// the magnetic field cell is used (and potentially updated) /// @param [in] pos is the field position - Vector3D getField(Cache& cache, const Vector3D& pos) const + Vector3D + getField(Cache& cache, const Vector3D& pos) const { - if (!cache.field_cache_ready || !cache.field_cache.isInside(pos)){ + if (!cache.field_cache_ready || !cache.field_cache.isInside(pos)) { cache.field_cache_ready = true; - cache.field_cache = m_bField.getFieldCell(pos); + cache.field_cache = m_bField.getFieldCell(pos); } // get the field from the cell cache.field = cache.field_cache.getField(pos); return cache.field; } - double step(Cache& cache, double& h) const { diff --git a/Core/include/ACTS/Propagator/EigenStepper.hpp b/Core/include/ACTS/Propagator/EigenStepper.hpp index 1008ca490b7..74a78983ea2 100644 --- a/Core/include/ACTS/Propagator/EigenStepper.hpp +++ b/Core/include/ACTS/Propagator/EigenStepper.hpp @@ -11,10 +11,10 @@ #include #include "ACTS/EventData/TrackParameters.hpp" +#include "ACTS/MagneticField/concept/AnyFieldLookup.hpp" #include "ACTS/Surfaces/Surface.hpp" #include "ACTS/Utilities/Definitions.hpp" #include "ACTS/Utilities/Units.hpp" -#include "ACTS/MagneticField/concept/AnyFieldLookup.hpp" namespace Acts { @@ -49,7 +49,7 @@ class EigenStepper { /// Constructor from the initial track parameters /// @tparam [in] par are the track parameters - /// + /// /// @note the covariance matrix is copied when needed template explicit Cache(const T& par) @@ -60,11 +60,11 @@ class EigenStepper { if (par.covariance()) { cov_transport = true; - cov = ActsSymMatrixD<5>(*par.covariance()); - + cov = ActsSymMatrixD<5>(*par.covariance()); + const double phi = dir.phi(); const double theta = dir.theta(); - + const auto transform = par.referenceSurface().transform().matrix(); jacobian(0, eLOC_0) = transform(0, eLOC_0); jacobian(0, eLOC_1) = transform(0, eLOC_1); @@ -107,20 +107,19 @@ class EigenStepper /// Jacobian used to transport the covariance matrix ActsMatrixD<7, 5> jacobian = ActsMatrixD<7, 5>::Zero(); - /// The propagation derivative + /// The propagation derivative ActsVectorD<7> derivative = ActsVectorD<7>::Zero(); /// Covariance matrix (and indicator) //// assocated with the initial error on track parameters - bool cov_transport = false; - ActsSymMatrixD<5> cov = ActsSymMatrixD<5>::Zero(); + bool cov_transport = false; + ActsSymMatrixD<5> cov = ActsSymMatrixD<5>::Zero(); - /// Lazily initialized cache - /// It caches the current magneticl field cell and stays interpolates within + /// Lazily initialized cache + /// It caches the current magneticl field cell and stays interpolates within /// as long as this is valid. See step() code for details. - bool field_cache_ready = false; + bool field_cache_ready = false; concept::AnyFieldCell<> field_cache; - }; // This struct is a meta-function which normally maps to BoundParameters... @@ -170,7 +169,7 @@ class EigenStepper EigenStepper(BField bField = BField()) : m_bField(std::move(bField)){}; /// Convert the propagation cache to curvilinear parameters - /// @param cache is the stepper cache + /// @param cache is the stepper cache /// @todo check: what if cache is already in courvilinear ? is this caught ? static CurvilinearParameters convert(const Cache& cache) @@ -236,7 +235,7 @@ class EigenStepper } /// Convert the propagation cache to track parameters at a certain surface - /// @param [in] cache Propagation cache used + /// @param [in] cache Propagation cache used /// @param [in] surface Destination surface to which the conversion is done template static BoundParameters @@ -303,11 +302,12 @@ class EigenStepper /// @param [in,out] cache is the propagation cache associated with the track /// the magnetic field cell is used (and potentially updated) /// @param [in] pos is the field position - Vector3D getField(Cache& cache, const Vector3D& pos) const + Vector3D + getField(Cache& cache, const Vector3D& pos) const { - if (!cache.field_cache_ready || !cache.field_cache.isInside(pos)){ + if (!cache.field_cache_ready || !cache.field_cache.isInside(pos)) { cache.field_cache_ready = true; - cache.field_cache = m_bField.getFieldCell(pos); + cache.field_cache = m_bField.getFieldCell(pos); } // get the field from the cell return std::move(cache.field_cache.getField(pos)); @@ -348,8 +348,9 @@ class EigenStepper // Second Runge-Kutta point const Vector3D pos1 = cache.pos + half_h * cache.dir + h2 / 8 * k1; - B_middle = getField(cache, pos1);; - k2 = qop * (cache.dir + half_h * k1).cross(B_middle); + B_middle = getField(cache, pos1); + ; + k2 = qop * (cache.dir + half_h * k1).cross(B_middle); // Third Runge-Kutta point k3 = qop * (cache.dir + half_h * k2).cross(B_middle); diff --git a/Core/include/ACTS/Propagator/Propagator.hpp b/Core/include/ACTS/Propagator/Propagator.hpp index eef0e4f7f23..b260219e6ad 100644 --- a/Core/include/ACTS/Propagator/Propagator.hpp +++ b/Core/include/ACTS/Propagator/Propagator.hpp @@ -91,10 +91,9 @@ namespace propagation { class Propagator final { public: - // Type of cache object used by the propagation implementation typedef typename Impl::template cache_type cache_type; - + /// @brief Options for propagate() call /// /// @tparam ObserverList List of observer types called after each @@ -115,7 +114,7 @@ namespace propagation { /// Maximum number of steps for one propagate() call unsigned int max_steps = 1000; - /// Required distance to surface + /// Required distance to surface double target_surface_distance = 1 * units::_um; /// Absolute minimum step size @@ -171,9 +170,8 @@ namespace propagation { /// template using obs_list_result_t = typename result_type_helper::type; - + public: - /// @brief Propagate track parameters - User method /// /// This function performs the propagation of the track parameters using the @@ -213,7 +211,7 @@ namespace propagation { "return track parameter type must be copy-constructible"); // Initialize the propagation result object - result_type r(Status::INPROGRESS); + result_type r(Status::INPROGRESS); // Compute the signed path limit and maximum step size const double signed_pathLimit @@ -243,7 +241,7 @@ namespace propagation { if (std::abs(stepMax) > std::abs(signed_pathLimit - r.pathLength)) stepMax = signed_pathLimit - r.pathLength; } - return r; + return r; } /// @brief Propagate track parameters - Expert method with propagation cache @@ -253,7 +251,7 @@ namespace propagation { /// is fulfilled, the destination surface is hit or the maximum number of /// steps/path length as given in the propagation options is reached. /// - /// It does check/re-use the propgation cache as much as possible for + /// It does check/re-use the propgation cache as much as possible for /// performance reasons /// /// @tparam TrackParameters Type of initial track parameters to propagate @@ -262,7 +260,7 @@ namespace propagation { /// @tparam AbortList Type list of abort conditions /// /// @param [in] cache Stepper cache built/updated from the start parameters - /// @param [in] target Target surface of to propagate to + /// @param [in] target Target surface of to propagate to /// @param [in] options Propagation options /// /// @return Propagation result containing the propagation status, final @@ -275,12 +273,12 @@ namespace propagation { obs_list_result_t< typename Impl::template return_parameter_type, ObserverList> - propagate_with_cache(cache_type& cache, - const TrackParameters& start, - const Surface& target, + propagate_with_cache(cache_type& cache, + const TrackParameters& start, + const Surface& target, const Options& options) const { - + // Type of track parameters produced at the end of the propagation typedef typename Impl::template return_parameter_type @@ -294,7 +292,7 @@ namespace propagation { "return track parameter type must be copy-constructible"); // Initialize the propagation result object - result_type r(Status::INPROGRESS); + result_type r(Status::INPROGRESS); // Compute the signed path limit and maximum step size const double signed_pathLimit @@ -327,7 +325,7 @@ namespace propagation { distance = m_impl.distance(target, cache.position(), cache.direction()); // Is it time to stop the propagation ? - if (std::abs(distance) < options.target_surface_distance + if (std::abs(distance) < options.target_surface_distance || std::abs(r.pathLength) >= options.max_path_length || options.stop_conditions(r, cache, stepMax)) { // Compute the final results and mark the propagation as successful @@ -361,7 +359,7 @@ namespace propagation { /// @tparam AbortList Type list of abort conditions /// /// @param [in] start Initial track parameters to propagate - /// @param [in] target Target surface of to propagate to + /// @param [in] target Target surface of to propagate to /// @param [in] options Propagation options /// /// @return Propagation result containing the propagation status, final @@ -379,7 +377,7 @@ namespace propagation { const Options& options) const { // Initialize the internal propagation cache - cache_type cache(start); + cache_type cache(start); return propagate_with_cache(cache, start, target, options); } From 7c1d62f5224eaa1de869095ea61f1fe0972844dc Mon Sep 17 00:00:00 2001 From: Andreas Salzburger Date: Tue, 17 Oct 2017 10:46:30 +0200 Subject: [PATCH 07/11] Bring RungeKuttaEngine and Atlas/EigenStepper into same unit system using units::Nat2Si Publish EigenStepper::Cache for new propagate_with_cache interface --- .../Extrapolation/detail/RungeKuttaEngine.ipp | 82 +++++++++---------- Core/include/ACTS/Propagator/EigenStepper.hpp | 67 ++++++++------- Core/src/Extrapolation/RungeKuttaUtils.cpp | 4 +- 3 files changed, 77 insertions(+), 76 deletions(-) diff --git a/Core/include/ACTS/Extrapolation/detail/RungeKuttaEngine.ipp b/Core/include/ACTS/Extrapolation/detail/RungeKuttaEngine.ipp index aa4fd95bf6c..d0543f85bae 100644 --- a/Core/include/ACTS/Extrapolation/detail/RungeKuttaEngine.ipp +++ b/Core/include/ACTS/Extrapolation/detail/RungeKuttaEngine.ipp @@ -500,8 +500,9 @@ Acts::RungeKuttaEngine::propagateWithJacobian( SA[0] = SA[1] = SA[2] = 0.; pCache.maxPathLimit = false; - if (pCache.mcondition && std::abs(pCache.pVector[6]) > 100.) return false; - + // @todo the inverse momentum condition is hard-coded + if (pCache.mcondition && std::abs(pCache.pVector[6]) > 100000.) return false; + // Step estimation until surface bool Q; double S, step = m_rkUtils.stepEstimator(kind, Su, pCache.pVector, Q); @@ -646,10 +647,8 @@ Acts::RungeKuttaEngine::rungeKuttaStep(int navigationStep, double* R = &(pCache.pVector[0]); // Coordinates double* A = &(pCache.pVector[3]); // Directions double* sA = &(pCache.pVector[42]); - - // @todo bring the numbers into initialize - double scaleM = 1e-3 / units::_T; - double Pi = 149.89626 * pCache.pVector[6] * units::_MeV; // Invert mometum/2. + // invert momentum in right units + double Pi = 0.5 / units::Nat2SI(1. / pCache.pVector[6]); double dltm = m_cfg.dlt * .03; double f0[3], f[3]; @@ -657,9 +656,9 @@ Acts::RungeKuttaEngine::rungeKuttaStep(int navigationStep, // if new field is required get it if (pCache.newfield) { Vector3D bField = m_cfg.fieldService->getField(Vector3D(R[0], R[1], R[2])); - f0[0] = scaleM * bField.x(); - f0[1] = scaleM * bField.y(); - f0[2] = scaleM * bField.z(); + f0[0] = bField.x(); + f0[1] = bField.y(); + f0[2] = bField.z(); } else { f0[0] = pCache.field[0]; f0[1] = pCache.field[1]; @@ -690,9 +689,9 @@ Acts::RungeKuttaEngine::rungeKuttaStep(int navigationStep, double gP[3] = {R[0] + A1 * S4, R[1] + B1 * S4, R[2] + C1 * S4}; Vector3D bField = m_cfg.fieldService->getField(Vector3D(gP[0], gP[1], gP[2])); - f[0] = scaleM * bField.x(); - f[1] = scaleM * bField.y(); - f[2] = scaleM * bField.z(); + f[0] = bField.x(); + f[1] = bField.y(); + f[2] = bField.z(); } else { f[0] = f0[0]; f[1] = f0[1]; @@ -716,9 +715,9 @@ Acts::RungeKuttaEngine::rungeKuttaStep(int navigationStep, double gP[3] = {R[0] + S * A4, R[1] + S * B4, R[2] + S * C4}; Vector3D bField = m_cfg.fieldService->getField(Vector3D(gP[0], gP[1], gP[2])); - f[0] = scaleM * bField.x(); - f[1] = scaleM * bField.y(); - f[2] = scaleM * bField.z(); + f[0] = bField.x(); + f[1] = bField.y(); + f[2] = bField.z(); } else { f[0] = f0[0]; f[1] = f0[1]; @@ -890,10 +889,7 @@ Acts::RungeKuttaEngine::rungeKuttaStepWithGradient( double* A = &(pCache.pVector[3]); // Directions double* sA = &(pCache.pVector[42]); // express the conversion factor with units - // @todo bring the conversion into initialize - double scaleM = 1e-3 / units::_T; - double Pi = 149.89626 * pCache.pVector[6] * units::_MeV; // Invert mometum/2. - + double Pi = 0.5 / units::Nat2SI(1. / pCache.pVector[6]); double dltm = m_cfg.dlt * .03; double f0[3], f1[3], f2[3], g0[9], g1[9], g2[9], H0[12], H1[12], H2[12]; @@ -902,18 +898,18 @@ Acts::RungeKuttaEngine::rungeKuttaStepWithGradient( deriv0 << 0., 0., 0., 0., 0., 0., 0., 0., 0.; Vector3D bField0 = m_cfg.fieldService->getFieldGradient( Vector3D(R[0], R[1], R[2]), deriv0); - f0[0] = scaleM * bField0.x(); - f0[1] = scaleM * bField0.y(); - f0[2] = scaleM * bField0.z(); - g0[0] = scaleM * deriv0(0, 0); - g0[1] = scaleM * deriv0(0, 1); - g0[2] = scaleM * deriv0(0, 2); - g0[3] = scaleM * deriv0(1, 0); - g0[4] = scaleM * deriv0(1, 1); - g0[5] = scaleM * deriv0(1, 2); - g0[6] = scaleM * deriv0(2, 0); - g0[7] = scaleM * deriv0(2, 1); - g0[8] = scaleM * deriv0(2, 2); + f0[0] = bField0.x(); + f0[1] = bField0.y(); + f0[2] = bField0.z(); + g0[0] = deriv0(0, 0); + g0[1] = deriv0(0, 1); + g0[2] = deriv0(0, 2); + g0[3] = deriv0(1, 0); + g0[4] = deriv0(1, 1); + g0[5] = deriv0(1, 2); + g0[6] = deriv0(2, 0); + g0[7] = deriv0(2, 1); + g0[8] = deriv0(2, 2); while (S != 0.) { double S3 = C33 * S, S4 = .25 * S, PS2 = Pi * S; @@ -940,18 +936,18 @@ Acts::RungeKuttaEngine::rungeKuttaStepWithGradient( deriv1 << 0., 0., 0., 0., 0., 0., 0., 0., 0.; Vector3D bField1 = m_cfg.fieldService->getFieldGradient( Vector3D(gP1[0], gP1[1], gP1[2]), deriv1); - f1[0] = scaleM * bField1.x(); - f1[1] = scaleM * bField1.y(); - f1[2] = scaleM * bField1.z(); - g1[0] = scaleM * deriv1(0, 0); - g1[1] = scaleM * deriv1(0, 1); - g1[2] = scaleM * deriv1(0, 2); - g1[3] = scaleM * deriv1(1, 0); - g1[4] = scaleM * deriv1(1, 1); - g1[5] = scaleM * deriv1(1, 2); - g1[6] = scaleM * deriv1(2, 0); - g1[7] = scaleM * deriv1(2, 1); - g1[8] = scaleM * deriv1(2, 2); + f1[0] = bField1.x(); + f1[1] = bField1.y(); + f1[2] = bField1.z(); + g1[0] = deriv1(0, 0); + g1[1] = deriv1(0, 1); + g1[2] = deriv1(0, 2); + g1[3] = deriv1(1, 0); + g1[4] = deriv1(1, 1); + g1[5] = deriv1(1, 2); + g1[6] = deriv1(2, 0); + g1[7] = deriv1(2, 1); + g1[8] = deriv1(2, 2); H1[0] = f1[0] * PS2; H1[1] = f1[1] * PS2; H1[2] = f1[2] * PS2; diff --git a/Core/include/ACTS/Propagator/EigenStepper.hpp b/Core/include/ACTS/Propagator/EigenStepper.hpp index 74a78983ea2..16ead69aad1 100644 --- a/Core/include/ACTS/Propagator/EigenStepper.hpp +++ b/Core/include/ACTS/Propagator/EigenStepper.hpp @@ -43,8 +43,41 @@ cross(const ActsMatrixD<3, 3>& m, const Vector3D& v) template class EigenStepper { + private: - /// Internal cache for track parameter propagation + // This struct is a meta-function which normally maps to BoundParameters... + template + struct s + { + typedef BoundParameters type; + }; + + // ...unless type S is int, in which case it maps to Curvilinear parameters + template + struct s + { + typedef CurvilinearParameters type; + }; + + /// Rotation matrix going from global coordinates to local surface coordinates + static ActsMatrixD<3, 3> + dLocaldGlobal(const Surface& p, const Vector3D& gpos) + { + // A surface's associated transform is a 4x4 affine transformation matrix, + // whose top-left corner is the 3x3 linear local-to-global rotation matrix, + // and whose right column is a translation vector, with a trailing 1. + // + // So to get the global-to-local rotation matrix, we only need to take the + // transpose of the top-left corner of the surface's associated transform. + // + return p.transform().matrix().topLeftCorner<3, 3>().transpose(); + } + +public: + /// Cache for track parameter propagation + /// + /// it is exposed to public for use of the expert-only + /// propagate_with_cache method of the propagator struct Cache { /// Constructor from the initial track parameters @@ -121,36 +154,7 @@ class EigenStepper bool field_cache_ready = false; concept::AnyFieldCell<> field_cache; }; - - // This struct is a meta-function which normally maps to BoundParameters... - template - struct s - { - typedef BoundParameters type; - }; - - // ...unless type S is int, in which case it maps to Curvilinear parameters - template - struct s - { - typedef CurvilinearParameters type; - }; - - /// Rotation matrix going from global coordinates to local surface coordinates - static ActsMatrixD<3, 3> - dLocaldGlobal(const Surface& p, const Vector3D& gpos) - { - // A surface's associated transform is a 4x4 affine transformation matrix, - // whose top-left corner is the 3x3 linear local-to-global rotation matrix, - // and whose right column is a translation vector, with a trailing 1. - // - // So to get the global-to-local rotation matrix, we only need to take the - // transpose of the top-left corner of the surface's associated transform. - // - return p.transform().matrix().topLeftCorner<3, 3>().transpose(); - } - -public: + /// Always use the same propagation cache type, independently of the initial /// track parameter type and of the target surface template @@ -449,6 +453,7 @@ class EigenStepper return h; } + private: /// Magnetic field inside of the detector BField m_bField; diff --git a/Core/src/Extrapolation/RungeKuttaUtils.cpp b/Core/src/Extrapolation/RungeKuttaUtils.cpp index b92f00b8a9a..ce80dc6727b 100644 --- a/Core/src/Extrapolation/RungeKuttaUtils.cpp +++ b/Core/src/Extrapolation/RungeKuttaUtils.cpp @@ -576,7 +576,7 @@ Acts::RungeKuttaUtils::stepEstimator(int kind, if (kind == 0) return stepEstimatorToStraw(Su, P, Q); if (kind == 2) return stepEstimatorToCylinder(Su, P, Q); if (kind == 3) return stepEstimatorToCone(Su, P, Q); - return 1000000.; + return 1000000.; // @todo maximum not hard-coded } ///////////////////////////////////////////////////////////////////////////////// @@ -594,7 +594,7 @@ Acts::RungeKuttaUtils::stepEstimatorToPlane(double* S, double A = a[0] * S[0] + a[1] * S[1] + a[2] * S[2]; if (A == 0.) { Q = false; - return 1000000.; + return 1000000.; // @todo maximum not hard-coded } double D = (S[3] - r[0] * S[0]) - (r[1] * S[1] + r[2] * S[2]); Q = true; From caed76a6cb7cc60d41571a4c1a1a0f3c52055979 Mon Sep 17 00:00:00 2001 From: Andreas Salzburger Date: Tue, 17 Oct 2017 11:49:50 +0200 Subject: [PATCH 08/11] more info for integration test --- .../Extrapolation/detail/RungeKuttaEngine.ipp | 6 ++-- Core/include/ACTS/Propagator/EigenStepper.hpp | 35 +++++++++++++++---- Core/include/ACTS/Propagator/Propagator.hpp | 3 ++ Core/src/Extrapolation/RungeKuttaUtils.cpp | 4 +-- IntegrationTests/PropagationTests.cpp | 1 + 5 files changed, 38 insertions(+), 11 deletions(-) diff --git a/Core/include/ACTS/Extrapolation/detail/RungeKuttaEngine.ipp b/Core/include/ACTS/Extrapolation/detail/RungeKuttaEngine.ipp index d0543f85bae..c244c8f54ce 100644 --- a/Core/include/ACTS/Extrapolation/detail/RungeKuttaEngine.ipp +++ b/Core/include/ACTS/Extrapolation/detail/RungeKuttaEngine.ipp @@ -502,7 +502,7 @@ Acts::RungeKuttaEngine::propagateWithJacobian( // @todo the inverse momentum condition is hard-coded if (pCache.mcondition && std::abs(pCache.pVector[6]) > 100000.) return false; - + // Step estimation until surface bool Q; double S, step = m_rkUtils.stepEstimator(kind, Su, pCache.pVector, Q); @@ -648,7 +648,7 @@ Acts::RungeKuttaEngine::rungeKuttaStep(int navigationStep, double* A = &(pCache.pVector[3]); // Directions double* sA = &(pCache.pVector[42]); // invert momentum in right units - double Pi = 0.5 / units::Nat2SI(1. / pCache.pVector[6]); + double Pi = 0.5 / units::Nat2SI(1. / pCache.pVector[6]); double dltm = m_cfg.dlt * .03; double f0[3], f[3]; @@ -889,7 +889,7 @@ Acts::RungeKuttaEngine::rungeKuttaStepWithGradient( double* A = &(pCache.pVector[3]); // Directions double* sA = &(pCache.pVector[42]); // express the conversion factor with units - double Pi = 0.5 / units::Nat2SI(1. / pCache.pVector[6]); + double Pi = 0.5 / units::Nat2SI(1. / pCache.pVector[6]); double dltm = m_cfg.dlt * .03; double f0[3], f1[3], f2[3], g0[9], g1[9], g2[9], H0[12], H1[12], H2[12]; diff --git a/Core/include/ACTS/Propagator/EigenStepper.hpp b/Core/include/ACTS/Propagator/EigenStepper.hpp index 16ead69aad1..218651a8025 100644 --- a/Core/include/ACTS/Propagator/EigenStepper.hpp +++ b/Core/include/ACTS/Propagator/EigenStepper.hpp @@ -43,7 +43,7 @@ cross(const ActsMatrixD<3, 3>& m, const Vector3D& v) template class EigenStepper { - + private: // This struct is a meta-function which normally maps to BoundParameters... template @@ -72,16 +72,16 @@ class EigenStepper // return p.transform().matrix().topLeftCorner<3, 3>().transpose(); } - + public: /// Cache for track parameter propagation - /// + /// /// it is exposed to public for use of the expert-only /// propagate_with_cache method of the propagator struct Cache { /// Constructor from the initial track parameters - /// @tparam [in] par are the track parameters + /// @tparam [in] par The track parameters at start /// /// @note the covariance matrix is copied when needed template @@ -90,6 +90,30 @@ class EigenStepper , dir(par.momentum().normalized()) , qop(par.charge() / par.momentum().norm()) , cov_transport(false) + { + update_covariance(par); + } + + /// The cache update for optimal performance + /// @tparam [in] par The new track parameters at start + /// + /// @todo check to identify an reuse of start/cache + template + void + update(const T& par) + { + pos = par.position(); + dir = par.momentum().normalized(); + qop = (par.charge() / par.momentum().norm()); + cov_transport = false; + update_covariance(par); + } + + /// The covariance update + /// @tparam [in] par The (new) track parameters at start + template + void + update_covariance(const T& par) { if (par.covariance()) { cov_transport = true; @@ -154,7 +178,7 @@ class EigenStepper bool field_cache_ready = false; concept::AnyFieldCell<> field_cache; }; - + /// Always use the same propagation cache type, independently of the initial /// track parameter type and of the target surface template @@ -453,7 +477,6 @@ class EigenStepper return h; } - private: /// Magnetic field inside of the detector BField m_bField; diff --git a/Core/include/ACTS/Propagator/Propagator.hpp b/Core/include/ACTS/Propagator/Propagator.hpp index b260219e6ad..6cd1716db1f 100644 --- a/Core/include/ACTS/Propagator/Propagator.hpp +++ b/Core/include/ACTS/Propagator/Propagator.hpp @@ -284,6 +284,9 @@ namespace propagation { Surface> return_parameter_type; + // update the cache if necessary + cache.update(start); + // Type of the full propagation result, including output from observers typedef obs_list_result_t result_type; diff --git a/Core/src/Extrapolation/RungeKuttaUtils.cpp b/Core/src/Extrapolation/RungeKuttaUtils.cpp index ce80dc6727b..e380fd823cb 100644 --- a/Core/src/Extrapolation/RungeKuttaUtils.cpp +++ b/Core/src/Extrapolation/RungeKuttaUtils.cpp @@ -576,7 +576,7 @@ Acts::RungeKuttaUtils::stepEstimator(int kind, if (kind == 0) return stepEstimatorToStraw(Su, P, Q); if (kind == 2) return stepEstimatorToCylinder(Su, P, Q); if (kind == 3) return stepEstimatorToCone(Su, P, Q); - return 1000000.; // @todo maximum not hard-coded + return 1000000.; // @todo maximum not hard-coded } ///////////////////////////////////////////////////////////////////////////////// @@ -594,7 +594,7 @@ Acts::RungeKuttaUtils::stepEstimatorToPlane(double* S, double A = a[0] * S[0] + a[1] * S[1] + a[2] * S[2]; if (A == 0.) { Q = false; - return 1000000.; // @todo maximum not hard-coded + return 1000000.; // @todo maximum not hard-coded } double D = (S[3] - r[0] * S[0]) - (r[1] * S[1] + r[2] * S[2]); Q = true; diff --git a/IntegrationTests/PropagationTests.cpp b/IntegrationTests/PropagationTests.cpp index 3230c1862f0..8cae814e8c6 100644 --- a/IntegrationTests/PropagationTests.cpp +++ b/IntegrationTests/PropagationTests.cpp @@ -239,6 +239,7 @@ namespace IntegrationTest { if ((calculated_cov - *tp->covariance()).norm() / std::min(calculated_cov.norm(), tp->covariance()->norm()) > 2e-7) { + std::cout << "initial parameters = " << tp->parameters() << std::endl; std::cout << "calculated = " << calculated_cov << std::endl << std::endl; std::cout << "obtained = " << *tp->covariance() << std::endl; } From cde6074016840cefcd049058c5ce9757992a4b18 Mon Sep 17 00:00:00 2001 From: Andreas Salzburger Date: Wed, 18 Oct 2017 10:56:55 +0200 Subject: [PATCH 09/11] update() method for Stepper --- Core/include/ACTS/Propagator/EigenStepper.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Core/include/ACTS/Propagator/EigenStepper.hpp b/Core/include/ACTS/Propagator/EigenStepper.hpp index 218651a8025..0cd1ce4da13 100644 --- a/Core/include/ACTS/Propagator/EigenStepper.hpp +++ b/Core/include/ACTS/Propagator/EigenStepper.hpp @@ -118,10 +118,9 @@ class EigenStepper if (par.covariance()) { cov_transport = true; cov = ActsSymMatrixD<5>(*par.covariance()); - const double phi = dir.phi(); const double theta = dir.theta(); - + // @todo - check this might have to be the measuremetn frame ... const auto transform = par.referenceSurface().transform().matrix(); jacobian(0, eLOC_0) = transform(0, eLOC_0); jacobian(0, eLOC_1) = transform(0, eLOC_1); From d349cde7f0ab6509ce8a89e62563ae069d0b2d3e Mon Sep 17 00:00:00 2001 From: Paul Gessinger Date: Thu, 19 Oct 2017 13:19:38 +0200 Subject: [PATCH 10/11] small format fix --- Core/include/ACTS/Propagator/EigenStepper.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Core/include/ACTS/Propagator/EigenStepper.hpp b/Core/include/ACTS/Propagator/EigenStepper.hpp index 0cd1ce4da13..234759fc3ca 100644 --- a/Core/include/ACTS/Propagator/EigenStepper.hpp +++ b/Core/include/ACTS/Propagator/EigenStepper.hpp @@ -116,11 +116,11 @@ class EigenStepper update_covariance(const T& par) { if (par.covariance()) { - cov_transport = true; - cov = ActsSymMatrixD<5>(*par.covariance()); + cov_transport = true; + cov = ActsSymMatrixD<5>(*par.covariance()); const double phi = dir.phi(); const double theta = dir.theta(); - // @todo - check this might have to be the measuremetn frame ... + // @todo - check this might have to be the measuremetn frame ... const auto transform = par.referenceSurface().transform().matrix(); jacobian(0, eLOC_0) = transform(0, eLOC_0); jacobian(0, eLOC_1) = transform(0, eLOC_1); From 0f308223f5750446ab4ad23516a5b2f2530ac107 Mon Sep 17 00:00:00 2001 From: Andreas Salzburger Date: Fri, 20 Oct 2017 10:34:04 +0200 Subject: [PATCH 11/11] Enforce shared BField to be const --- Core/include/ACTS/MagneticField/SharedBFieldMap.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/Core/include/ACTS/MagneticField/SharedBFieldMap.hpp b/Core/include/ACTS/MagneticField/SharedBFieldMap.hpp index de00c868897..9a4ba99e2da 100644 --- a/Core/include/ACTS/MagneticField/SharedBFieldMap.hpp +++ b/Core/include/ACTS/MagneticField/SharedBFieldMap.hpp @@ -24,8 +24,9 @@ class SharedBField { public: /// @brief the constructur with a shared pointer + /// @note since it is a shared field, we enforce it to be const /// @tparam bField is the shared BField to be stored - SharedBField(std::shared_ptr bField) : m_bField(bField) {} + SharedBField(std::shared_ptr bField) : m_bField(bField) {} /// @brief retrieve magnetic field value /// @@ -67,7 +68,7 @@ class SharedBField } private: - std::shared_ptr m_bField; + std::shared_ptr m_bField; }; } // namespace Acts