Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-04-04 07:47:35

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/Material/IVolumeMaterial.hpp"
0013 #include "Acts/Material/Material.hpp"
0014 #include "Acts/Utilities/BinUtility.hpp"
0015 #include "Acts/Utilities/Interpolation.hpp"
0016 
0017 #include <functional>
0018 #include <optional>
0019 
0020 namespace Acts {
0021 
0022 /// @addtogroup material
0023 /// @{
0024 
0025 /// @brief Struct for mapping global 3D positions to material values
0026 ///
0027 /// Global 3D positions are transformed into a @c DIM_POS Dimensional vector
0028 /// which is used to look up the material classification value in the
0029 /// underlying material map.
0030 template <typename G>
0031 struct MaterialMapLookup {
0032  public:
0033   /// Type alias for material grid
0034   using Grid_t = G;
0035   /// Dimensionality of the position space for material interpolation
0036   static constexpr std::size_t DIM_POS = Grid_t::DIM;
0037 
0038   /// @brief Struct representing smallest grid unit in material grid
0039   ///
0040   /// This type encapsulate all required information to perform linear
0041   /// interpolation of material classification values within a 3D volume.
0042   struct MaterialCell {
0043     /// Number of corner points defining the confining hyper-box
0044     static constexpr unsigned int N = 1 << DIM_POS;
0045 
0046     /// @brief Default constructor
0047     ///
0048     /// @param [in] transformPos   Mapping of global 3D coordinates onto grid
0049     /// space
0050     /// @param [in] lowerLeft   Generalized lower-left corner of hyper box
0051     ///                         (containing the minima of the hyper box along
0052     ///                         each Dimension)
0053     /// @param [in] upperRight  Generalized upper-right corner of hyper box
0054     ///                         (containing the maxima of the hyper box along
0055     ///                         each Dimension)
0056     /// @param [in] materialValues Material classification values at the hyper
0057     /// box corners sorted in the canonical order defined in Acts::interpolate
0058     MaterialCell(std::function<Vector<DIM_POS>(const Vector3&)> transformPos,
0059                  std::array<double, DIM_POS> lowerLeft,
0060                  std::array<double, DIM_POS> upperRight,
0061                  std::array<Material::ParametersVector, N> materialValues)
0062         : m_transformPos(std::move(transformPos)),
0063           m_lowerLeft(std::move(lowerLeft)),
0064           m_upperRight(std::move(upperRight)),
0065           m_materialValues(std::move(materialValues)) {}
0066 
0067     /// @brief Retrieve material at given position
0068     ///
0069     /// @param [in] position Global 3D position
0070     /// @return Material at the given position
0071     ///
0072     /// @pre The given @c position must lie within the current cell.
0073     Material getMaterial(const Vector3& position) const {
0074       // defined in Interpolation.hpp
0075       return Material(interpolate(m_transformPos(position), m_lowerLeft,
0076                                   m_upperRight, m_materialValues));
0077     }
0078 
0079     /// @brief Check whether given 3D position is inside this cell
0080     ///
0081     /// @param [in] position Global 3D position
0082     /// @return @c true if position is inside the current cell,
0083     ///         otherwise @c false
0084     bool isInside(const Vector3& position) const {
0085       const auto& gridCoordinates = m_transformPos(position);
0086       for (unsigned int i = 0; i < DIM_POS; ++i) {
0087         if (gridCoordinates[i] < m_lowerLeft.at(i) ||
0088             gridCoordinates[i] >= m_upperRight.at(i)) {
0089           return false;
0090         }
0091       }
0092       return true;
0093     }
0094 
0095    private:
0096     /// Geometric transformation applied to global 3D positions
0097     std::function<Vector<DIM_POS>(const Vector3&)> m_transformPos;
0098 
0099     /// Generalized lower-left corner of the confining hyper-box
0100     std::array<double, DIM_POS> m_lowerLeft;
0101 
0102     /// Generalized upper-right corner of the confining hyper-box
0103     std::array<double, DIM_POS> m_upperRight;
0104 
0105     /// @brief Material component vectors at the hyper-box corners
0106     ///
0107     /// @note These values must be order according to the prescription detailed
0108     ///       in Acts::interpolate.
0109     std::array<Material::ParametersVector, N> m_materialValues;
0110   };
0111 
0112   /// @brief Default constructor
0113   ///
0114   /// @param [in] transformPos Mapping of global 3D coordinates (cartesian)
0115   /// onto grid space
0116   /// @param [in] grid Grid storing material classification values
0117   MaterialMapLookup(std::function<Vector<DIM_POS>(const Vector3&)> transformPos,
0118                     Grid_t grid)
0119       : m_transformPos(std::move(transformPos)), m_grid(std::move(grid)) {}
0120 
0121   /// @brief Retrieve binned material at given position
0122   ///
0123   /// @param [in] position Global 3D position
0124   /// @return Material at the given position
0125   ///
0126   /// @pre The given @c position must lie within the range of the underlying
0127   /// map.
0128   Material material(const Vector3& position) const {
0129     return Material(m_grid.atLocalBins(
0130         m_grid.localBinsFromLowerLeftEdge(m_transformPos(position))));
0131   }
0132 
0133   /// @brief Retrieve interpolated material at given position
0134   ///
0135   /// @param [in] position Global 3D position
0136   /// @return Material at the given position
0137   ///
0138   /// @pre The given @c position must lie within the range of the underlying
0139   /// map.
0140   Material getMaterial(const Vector3& position) const {
0141     return Material(m_grid.interpolate(m_transformPos(position)));
0142   }
0143 
0144   /// @brief Retrieve material cell for given position
0145   ///
0146   /// @param [in] position Global 3D position
0147   /// @return material cell containing the given global position
0148   ///
0149   /// @pre The given @c position must lie within the range of the underlying
0150   /// map.
0151   MaterialCell getMaterialCell(const Vector3& position) const {
0152     const auto& gridPosition = m_transformPos(position);
0153     std::size_t bin = m_grid.globalBinFromPosition(gridPosition);
0154     const auto& indices = m_grid.localBinsFromPosition(bin);
0155     const auto& lowerLeft = m_grid.lowerLeftBinEdge(indices);
0156     const auto& upperRight = m_grid.upperRightBinEdge(indices);
0157 
0158     // Loop through all corner points
0159     constexpr std::size_t nCorners = 1 << DIM_POS;
0160     std::array<Material::ParametersVector, nCorners> neighbors{};
0161     const auto& cornerIndices = m_grid.closestPointsIndices(gridPosition);
0162 
0163     std::size_t i = 0;
0164     for (std::size_t index : cornerIndices) {
0165       neighbors.at(i++) = m_grid.at(index);
0166     }
0167 
0168     return MaterialCell(m_transformPos, lowerLeft, upperRight,
0169                         std::move(neighbors));
0170   }
0171 
0172   /// @brief Get the number of bins for all axes of the map
0173   ///
0174   /// @return Vector returning number of bins for all map axes
0175   std::vector<std::size_t> getNBins() const {
0176     auto nBinsArray = m_grid.numLocalBins();
0177     return std::vector<std::size_t>(nBinsArray.begin(), nBinsArray.end());
0178   }
0179 
0180   /// @brief Get the minimum value of all axes of the map
0181   ///
0182   /// @return Vector returning the minima of all map axes
0183   std::vector<double> getMin() const {
0184     auto minArray = m_grid.minPosition();
0185     return std::vector<double>(minArray.begin(), minArray.end());
0186   }
0187 
0188   /// @brief Get the maximum value of all axes of the map
0189   ///
0190   /// @return Vector returning the maxima of all map axes
0191   std::vector<double> getMax() const {
0192     auto maxArray = m_grid.maxPosition();
0193     return std::vector<double>(maxArray.begin(), maxArray.end());
0194   }
0195 
0196   /// @brief Check whether given 3D position is inside look-up domain
0197   ///
0198   /// @param [in] position Global 3D position
0199   /// @return @c true if position is inside the defined look-up grid,
0200   ///         otherwise @c false
0201   bool isInside(const Vector3& position) const {
0202     return m_grid.isInside(m_transformPos(position));
0203   }
0204 
0205   /// @brief Get a const reference on the underlying grid structure
0206   ///
0207   /// @return Grid reference
0208   const Grid_t& getGrid() const { return m_grid; }
0209 
0210  private:
0211   /// Geometric transformation applied to global 3D positions
0212   std::function<Vector<DIM_POS>(const Vector3&)> m_transformPos;
0213   /// Grid storing material values
0214   Grid_t m_grid;
0215 };
0216 
0217 /// @brief Interpolate material classification values from material values on a
0218 /// given grid
0219 ///
0220 /// This class implements a material service which is initialized by a
0221 /// material map defined by:
0222 /// - a list of material values on a regular grid in some n-Dimensional space,
0223 /// - a transformation of global 3D coordinates onto this n-Dimensional
0224 /// space.
0225 /// - a transformation of local n-Dimensional material coordinates into
0226 /// global (cartesian) 3D coordinates
0227 ///
0228 /// The material value for a given global position is then determined by:
0229 /// - mapping the position onto the grid,
0230 /// - looking up the material classification values on the closest grid points,
0231 /// - doing a linear interpolation of these values.
0232 /// @warning Each classification number of the material is interpolated
0233 /// independently and thus does not consider any correlations that exists
0234 /// between these values. This might work out since the used material is already
0235 /// a mean of the materials in a certain bin and can therewith be treated as a
0236 /// collection of numbers.
0237 /// @tparam G Type of the grid
0238 template <typename Mapper_t>
0239 class InterpolatedMaterialMap : public IVolumeMaterial {
0240  public:
0241   /// @brief Temporary storage of a certain cell to improve material access
0242   struct Cache {
0243     /// Stored material cell
0244     std::optional<typename Mapper_t::MaterialCell> matCell;
0245     /// Boolean statement if the cell is initialized
0246     bool initialized = false;
0247   };
0248 
0249   /// @brief Create interpolated map
0250   ///
0251   /// @param [in] mapper Material map
0252   explicit InterpolatedMaterialMap(Mapper_t&& mapper)
0253       : m_mapper(std::move(mapper)) {}
0254 
0255   /// @brief Create interpolated map
0256   ///
0257   /// @param [in] mapper Material map
0258   /// @param [in] bu @c BinUtility for build from
0259   InterpolatedMaterialMap(Mapper_t&& mapper, BinUtility bu)
0260       : m_mapper(std::move(mapper)), m_binUtility(std::move(bu)) {}
0261 
0262   /// @brief Retrieve the binned material
0263   ///
0264   /// @param [in] position Global 3D position
0265   ///
0266   /// @return Material at given position
0267   const Material material(const Vector3& position) const override {
0268     return m_mapper.material(position);
0269   }
0270 
0271   /// @brief Retrieve the interpolated material
0272   ///
0273   /// @param [in] position Global 3D position
0274   ///
0275   /// @return material at given position
0276   Material getMaterial(const Vector3& position) const {
0277     return m_mapper.getMaterial(position);
0278   }
0279 
0280   /// @brief Retrieve material
0281   ///
0282   /// @param [in] position Global 3D position
0283   /// @param [in,out] cache Cache object. Contains material cell used for
0284   /// interpolation
0285   ///
0286   /// @return material at given position
0287   Material getMaterial(const Vector3& position, Cache& cache) const {
0288     if (!cache.initialized || !(*cache.matCell).isInside(position)) {
0289       cache.matCell = getMaterialCell(position);
0290       cache.initialized = true;
0291     }
0292     return (*cache.matCell).getMaterial(position);
0293   }
0294 
0295   /// @brief Retrieve material value & its "gradient"
0296   ///
0297   /// @param [in]  position   Global 3D position
0298   /// @return Material
0299   ///
0300   /// @note Currently the derivative is not calculated
0301   /// @todo return derivative
0302   Material getMaterialGradient(const Vector3& position,
0303                                Matrix<5, 5>& /*derivative*/) const {
0304     return m_mapper.getMaterial(position);
0305   }
0306 
0307   /// @brief Retrieve material value & its "gradient"
0308   ///
0309   /// @param [in]  position   Global 3D position
0310   /// @return Material
0311   ///
0312   /// @note Currently the derivative is not calculated
0313   /// @note Cache is not used currently
0314   /// @todo return derivative
0315   Material getMaterialGradient(const Vector3& position,
0316                                Matrix<5, 5>& /*derivative*/,
0317                                Cache& /*cache*/) const {
0318     return m_mapper.getMaterial(position);
0319   }
0320 
0321   /// @brief Convenience method to access underlying material mapper
0322   ///
0323   /// @return The material mapper
0324   const Mapper_t& getMapper() const { return m_mapper; }
0325 
0326   /// @brief Check whether given 3D position is inside look-up domain
0327   ///
0328   /// @param [in] position Global 3D position
0329   /// @return @c true if position is inside the defined map, otherwise @c false
0330   bool isInside(const Vector3& position) const {
0331     return m_mapper.isInside(position);
0332   }
0333 
0334   /// Return the BinUtility
0335   /// @return Const reference to the bin utility for the material map
0336   const BinUtility& binUtility() const { return m_binUtility; }
0337 
0338   /// Output Method for std::ostream
0339   ///
0340   /// @param sl The outoput stream
0341   /// @return Reference to the output stream for method chaining
0342   std::ostream& toStream(std::ostream& sl) const override {
0343     sl << "Acts::InterpolatedMaterialMap : " << std::endl;
0344     sl << "   - Number of Material bins [0,1] : " << m_binUtility.max(0) + 1
0345        << " / " << m_binUtility.max(1) + 1 << std::endl;
0346     sl << "   - Parse full update material    : " << std::endl;
0347     return sl;
0348   }
0349 
0350  private:
0351   /// @brief Retrieve cell for given position
0352   ///
0353   /// @param [in] position Global 3D position
0354   /// @return Material cell containing the given global position
0355   ///
0356   /// @pre The given @c position must lie within the range of the underlying
0357   /// map.
0358   typename Mapper_t::MaterialCell getMaterialCell(
0359       const Vector3& position) const {
0360     return m_mapper.getMaterialCell(position);
0361   }
0362 
0363   /// @brief object for global coordinate transformation and interpolation
0364   ///
0365   /// This object performs the mapping of the global 3D coordinates onto the
0366   /// material grid and the interpolation of the material component values on
0367   /// close-by grid points.
0368   Mapper_t m_mapper;
0369 
0370   BinUtility m_binUtility{};
0371 };
0372 
0373 /// @}
0374 }  // namespace Acts