Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-01-10 09:16:04

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