Acts/MagneticField/InterpolatedBFieldMap.hpp

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 #pragma once
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0013 #include "Acts/MagneticField/MagneticFieldError.hpp"
0014 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
0015 #include "Acts/Utilities/Interpolation.hpp"
0016 #include "Acts/Utilities/Result.hpp"
0017
0018 #include <functional>
0019 #include <optional>
0020 #include <vector>
0021
0022 namespace Acts {
0023
0024 class InterpolatedMagneticField : public MagneticFieldProvider {
0025  public:
0026   /// @brief get the number of bins for all axes of the field map
0027   ///
0028   /// @return vector returning number of bins for all field map axes
0029   virtual std::vector<std::size_t> getNBins() const = 0;
0030
0031   /// @brief get the minimum value of all axes of the field map
0032   ///
0033   /// @return vector returning the minima of all field map axes
0034   virtual std::vector<double> getMin() const = 0;
0035
0036   /// @brief get the maximum value of all axes of the field map
0037   ///
0038   /// @return vector returning the maxima of all field map axes
0039   virtual std::vector<double> getMax() const = 0;
0040
0041   /// @brief check whether given 3D position is inside look-up domain
0042   ///
0043   /// @param [in] position global 3D position
0044   /// @return @c true if position is inside the defined look-up grid,
0045   ///         otherwise @c false
0046   virtual bool isInside(const Vector3& position) const = 0;
0047
0048   /// Get a field value without checking if the lookup position is within the
0049   /// interpolation domain.
0050   ///
0051   /// @param position The lookup position in 3D
0052   /// @return The field value at @p position
0053   virtual Vector3 getFieldUnchecked(const Vector3& position) const = 0;
0054 };
0055
0056 /// @ingroup MagneticField
0057 /// @brief interpolate magnetic field value from field values on a given grid
0058 ///
0059 /// This class implements a magnetic field service which is initialized by a
0060 /// field map defined by:
0061 /// - a list of field values on a regular grid in some n-Dimensional space,
0062 /// - a transformation of global 3D coordinates onto this n-Dimensional
0063 /// space.
0064 /// - a transformation of local n-Dimensional magnetic field coordinates into
0065 /// global (cartesian) 3D coordinates
0066 ///
0067 /// The magnetic field value for a given global position is then determined by:
0068 /// - mapping the position onto the grid,
0069 /// - looking up the magnetic field values on the closest grid points,
0070 /// - doing a linear interpolation of these magnetic field values.
0071 ///
0072 /// @tparam grid_t The Grid type which provides the field storage and
0073 /// interpolation
0074 template <typename grid_t>
0075 class InterpolatedBFieldMap : public InterpolatedMagneticField {
0076  public:
0077   using Grid = grid_t;
0078   using FieldType = typename Grid::value_type;
0079   static constexpr std::size_t DIM_POS = Grid::DIM;
0080
0081   /// @brief struct representing smallest grid unit in magnetic field grid
0082   ///
0083   /// This type encapsulate all required information to perform linear
0084   /// interpolation of magnetic field values within a confined 3D volume.
0085   struct FieldCell {
0086     /// number of corner points defining the confining hyper-box
0087     static constexpr unsigned int N = 1 << DIM_POS;
0088
0089    public:
0090     /// @brief default constructor
0091     ///
0092     /// @param [in] lowerLeft   generalized lower-left corner of hyper box
0093     ///                         (containing the minima of the hyper box along
0094     ///                         each Dimension)
0095     /// @param [in] upperRight  generalized upper-right corner of hyper box
0096     ///                         (containing the maxima of the hyper box along
0097     ///                         each Dimension)
0098     /// @param [in] fieldValues field values at the hyper box corners sorted in
0099     ///                         the canonical order defined in Acts::interpolate
0100     FieldCell(std::array<double, DIM_POS> lowerLeft,
0101               std::array<double, DIM_POS> upperRight,
0102               std::array<Vector3, N> fieldValues)
0103         : m_lowerLeft(std::move(lowerLeft)),
0104           m_upperRight(std::move(upperRight)),
0105           m_fieldValues(std::move(fieldValues)) {}
0106
0107     /// @brief retrieve field at given position
0108     ///
0109     /// @param [in] position global 3D position
0110     /// @return magnetic field value at the given position
0111     ///
0112     /// @pre The given @c position must lie within the current field cell.
0113     Vector3 getField(const ActsVector<DIM_POS>& position) const {
0114       // defined in Interpolation.hpp
0115       return interpolate(position, m_lowerLeft, m_upperRight, m_fieldValues);
0116     }
0117
0118     /// @brief check whether given 3D position is inside this field cell
0119     ///
0120     /// @param [in] position global 3D position
0121     /// @return @c true if position is inside the current field cell,
0122     ///         otherwise @c false
0123     bool isInside(const ActsVector<DIM_POS>& position) const {
0124       for (unsigned int i = 0; i < DIM_POS; ++i) {
0125         if (position[i] < m_lowerLeft[i] || position[i] >= m_upperRight[i]) {
0126           return false;
0127         }
0128       }
0129       return true;
0130     }
0131
0132    private:
0133     /// generalized lower-left corner of the confining hyper-box
0134     std::array<double, DIM_POS> m_lowerLeft;
0135
0136     /// generalized upper-right corner of the confining hyper-box
0137     std::array<double, DIM_POS> m_upperRight;
0138
0139     /// @brief magnetic field vectors at the hyper-box corners
0140     ///
0141     /// @note These values must be order according to the prescription detailed
0142     ///       in Acts::interpolate.
0143     std::array<Vector3, N> m_fieldValues;
0144   };
0145
0146   struct Cache {
0147     /// @brief Constructor with magnetic field context
0148     Cache(const MagneticFieldContext& /*mctx*/) {}
0149
0150     std::optional<FieldCell> fieldCell;
0151     bool initialized = false;
0152   };
0153
0154   /// @brief  Config structure for the interpolated B field map
0155   struct Config {
0156     /// @brief mapping of global 3D coordinates (cartesian)
0157     /// onto grid space
0158     std::function<ActsVector<DIM_POS>(const Vector3&)> transformPos;
0159
0160     /// @brief calculating the global 3D coordinates
0161     /// (cartesian) of the magnetic field with the local n dimensional field and
0162     /// the global 3D position as input
0163     std::function<Vector3(const FieldType&, const Vector3&)> transformBField;
0164
0165     /// @brief grid storing magnetic field values
0166     Grid grid;
0167
0168     /// @brief global B-field scaling factor
0169     ///
0170     /// @note Negative values for @p scale are accepted and will invert the
0171     ///       direction of the magnetic field.
0172     double scale = 1.;
0173   };
0174
0175   /// @brief default constructor
0176   ///
0177   InterpolatedBFieldMap(Config cfg) : m_cfg{std::move(cfg)} {
0178     typename Grid::index_t minBin{};
0179     minBin.fill(1);
0180     m_lowerLeft = m_cfg.grid.lowerLeftBinEdge(minBin);
0181     m_upperRight = m_cfg.grid.lowerLeftBinEdge(m_cfg.grid.numLocalBins());
0182   }
0183
0184   /// @brief retrieve field cell for given position
0185   ///
0186   /// @param [in] position global 3D position
0187   /// @return field cell containing the given global position
0188   ///
0189   /// @pre The given @c position must lie within the range of the underlying
0190   ///      magnetic field map.
0191   Result<FieldCell> getFieldCell(const Vector3& position) const {
0192     const auto& gridPosition = m_cfg.transformPos(position);
0193     const auto& indices = m_cfg.grid.localBinsFromPosition(gridPosition);
0194     const auto& lowerLeft = m_cfg.grid.lowerLeftBinEdge(indices);
0195     const auto& upperRight = m_cfg.grid.upperRightBinEdge(indices);
0196
0197     // loop through all corner points
0198     constexpr std::size_t nCorners = 1 << DIM_POS;
0199     std::array<Vector3, nCorners> neighbors;
0200     const auto& cornerIndices = m_cfg.grid.closestPointsIndices(gridPosition);
0201
0202     if (!isInsideLocal(gridPosition)) {
0203       return MagneticFieldError::OutOfBounds;
0204     }
0205
0206     std::size_t i = 0;
0207     for (std::size_t index : cornerIndices) {
0208       neighbors.at(i++) = m_cfg.transformBField(m_cfg.grid.at(index), position);
0209     }
0210
0211     assert(i == nCorners);
0212
0213     return FieldCell(lowerLeft, upperRight, std::move(neighbors));
0214   }
0215
0216   /// @brief get the number of bins for all axes of the field map
0217   ///
0218   /// @return vector returning number of bins for all field map axes
0219   std::vector<std::size_t> getNBins() const final {
0220     auto nBinsArray = m_cfg.grid.numLocalBins();
0221     return std::vector<std::size_t>(nBinsArray.begin(), nBinsArray.end());
0222   }
0223
0224   /// @brief get the minimum value of all axes of the field map
0225   ///
0226   /// @return vector returning the minima of all field map axes
0227   std::vector<double> getMin() const final {
0228     return std::vector<double>(m_lowerLeft.begin(), m_lowerLeft.end());
0229   }
0230
0231   /// @brief get the maximum value of all axes of the field map
0232   ///
0233   /// @return vector returning the maxima of all field map axes
0234   std::vector<double> getMax() const final {
0235     return std::vector<double>(m_upperRight.begin(), m_upperRight.end());
0236   }
0237
0238   /// @brief check whether given 3D position is inside look-up domain
0239   ///
0240   /// @param [in] position global 3D position
0241   /// @return @c true if position is inside the defined look-up grid,
0242   ///         otherwise @c false
0243   bool isInside(const Vector3& position) const final {
0244     return isInsideLocal(m_cfg.transformPos(position));
0245   }
0246
0247   /// @brief check whether given 3D position is inside look-up domain
0248   ///
0249   /// @param [in] gridPosition local N-D position
0250   /// @return @c true if position is inside the defined look-up grid,
0251   ///         otherwise @c false
0252   bool isInsideLocal(const ActsVector<DIM_POS>& gridPosition) const {
0253     for (unsigned int i = 0; i < DIM_POS; ++i) {
0254       if (gridPosition[i] < m_lowerLeft[i] ||
0255           gridPosition[i] >= m_upperRight[i]) {
0256         return false;
0257       }
0258     }
0259     return true;
0260   }
0261
0262   /// @brief Get a const reference on the underlying grid structure
0263   ///
0264   /// @return grid reference
0265   const Grid& getGrid() const { return m_cfg.grid; }
0266
0267   /// @copydoc MagneticFieldProvider::makeCache(const MagneticFieldContext&) const
0268   MagneticFieldProvider::Cache makeCache(
0269       const MagneticFieldContext& mctx) const final {
0270     return MagneticFieldProvider::Cache{std::in_place_type<Cache>, mctx};
0271   }
0272
0273   /// @brief retrieve field at given position
0274   ///
0275   /// @param [in] position global 3D position
0276   /// @return magnetic field value at the given position
0277   ///
0278   /// @pre The given @c position must lie within the range of the underlying
0279   ///      magnetic field map.
0280   Result<Vector3> getField(const Vector3& position) const {
0281     const auto gridPosition = m_cfg.transformPos(position);
0282     if (!isInsideLocal(gridPosition)) {
0283       return Result<Vector3>::failure(MagneticFieldError::OutOfBounds);
0284     }
0285
0286     return Result<Vector3>::success(
0287         m_cfg.transformBField(m_cfg.grid.interpolate(gridPosition), position));
0288   }
0289
0290   Vector3 getFieldUnchecked(const Vector3& position) const final {
0291     const auto gridPosition = m_cfg.transformPos(position);
0292     return m_cfg.transformBField(m_cfg.grid.interpolate(gridPosition),
0293                                  position);
0294   }
0295
0296   /// @copydoc MagneticFieldProvider::getField(const Vector3&,MagneticFieldProvider::Cache&) const
0297   Result<Vector3> getField(const Vector3& position,
0298                            MagneticFieldProvider::Cache& cache) const final {
0299     Cache& lcache = cache.as<Cache>();
0300     const auto gridPosition = m_cfg.transformPos(position);
0301     if (!lcache.fieldCell || !(*lcache.fieldCell).isInside(gridPosition)) {
0302       auto res = getFieldCell(position);
0303       if (!res.ok()) {
0304         return Result<Vector3>::failure(res.error());
0305       }
0306       lcache.fieldCell = *res;
0307     }
0308     return Result<Vector3>::success((*lcache.fieldCell).getField(gridPosition));
0309   }
0310
0311  private:
0312   Config m_cfg;
0313
0314   typename Grid::point_t m_lowerLeft;
0315   typename Grid::point_t m_upperRight;
0316 };
0317
0318 }  // namespace Acts