Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-11-07 09:15:08

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 /// @brief Struct for mapping global 3D positions to material values
0023 ///
0024 /// Global 3D positions are transformed into a @c DIM_POS Dimensional vector
0025 /// which is used to look up the material classification value in the
0026 /// underlying material map.
0027 template <typename G>
0028 struct MaterialMapper {
0029  public:
0030   /// Type alias for material grid
0031   using Grid_t = G;
0032   /// Dimensionality of the position space for material interpolation
0033   static constexpr std::size_t DIM_POS = Grid_t::DIM;
0034 
0035   /// @brief Struct representing smallest grid unit in material grid
0036   ///
0037   /// This type encapsulate all required information to perform linear
0038   /// interpolation of material classification values within a 3D volume.
0039   struct MaterialCell {
0040     /// Number of corner points defining the confining hyper-box
0041     static constexpr unsigned int N = 1 << DIM_POS;
0042 
0043     /// @brief Default constructor
0044     ///
0045     /// @param [in] transformPos   Mapping of global 3D coordinates onto grid
0046     /// space
0047     /// @param [in] lowerLeft   Generalized lower-left corner of hyper box
0048     ///                         (containing the minima of the hyper box along
0049     ///                         each Dimension)
0050     /// @param [in] upperRight  Generalized upper-right corner of hyper box
0051     ///                         (containing the maxima of the hyper box along
0052     ///                         each Dimension)
0053     /// @param [in] materialValues Material classification values at the hyper
0054     /// box corners sorted in the canonical order defined in Acts::interpolate
0055     MaterialCell(
0056         std::function<ActsVector<DIM_POS>(const Vector3&)> transformPos,
0057         std::array<double, DIM_POS> lowerLeft,
0058         std::array<double, DIM_POS> upperRight,
0059         std::array<Material::ParametersVector, N> materialValues)
0060         : m_transformPos(std::move(transformPos)),
0061           m_lowerLeft(std::move(lowerLeft)),
0062           m_upperRight(std::move(upperRight)),
0063           m_materialValues(std::move(materialValues)) {}
0064 
0065     /// @brief Retrieve material at given position
0066     ///
0067     /// @param [in] position Global 3D position
0068     /// @return Material at the given position
0069     ///
0070     /// @pre The given @c position must lie within the current cell.
0071     Material getMaterial(const Vector3& position) const {
0072       // defined in Interpolation.hpp
0073       return Material(interpolate(m_transformPos(position), m_lowerLeft,
0074                                   m_upperRight, m_materialValues));
0075     }
0076 
0077     /// @brief Check whether given 3D position is inside this cell
0078     ///
0079     /// @param [in] position Global 3D position
0080     /// @return @c true if position is inside the current cell,
0081     ///         otherwise @c false
0082     bool isInside(const Vector3& position) const {
0083       const auto& gridCoordinates = m_transformPos(position);
0084       for (unsigned int i = 0; i < DIM_POS; ++i) {
0085         if (gridCoordinates[i] < m_lowerLeft.at(i) ||
0086             gridCoordinates[i] >= m_upperRight.at(i)) {
0087           return false;
0088         }
0089       }
0090       return true;
0091     }
0092 
0093    private:
0094     /// Geometric transformation applied to global 3D positions
0095     std::function<ActsVector<DIM_POS>(const Vector3&)> m_transformPos;
0096 
0097     /// Generalized lower-left corner of the confining hyper-box
0098     std::array<double, DIM_POS> m_lowerLeft;
0099 
0100     /// Generalized upper-right corner of the confining hyper-box
0101     std::array<double, DIM_POS> m_upperRight;
0102 
0103     /// @brief Material component vectors at the hyper-box corners
0104     ///
0105     /// @note These values must be order according to the prescription detailed
0106     ///       in Acts::interpolate.
0107     std::array<Material::ParametersVector, N> m_materialValues;
0108   };
0109 
0110   /// @brief Default constructor
0111   ///
0112   /// @param [in] transformPos Mapping of global 3D coordinates (cartesian)
0113   /// onto grid space
0114   /// @param [in] grid Grid storing material classification values
0115   MaterialMapper(
0116       std::function<ActsVector<DIM_POS>(const Vector3&)> transformPos,
0117       Grid_t grid)
0118       : m_transformPos(std::move(transformPos)), m_grid(std::move(grid)) {}
0119 
0120   /// @brief Retrieve binned material at given position
0121   ///
0122   /// @param [in] position Global 3D position
0123   /// @return Material at the given position
0124   ///
0125   /// @pre The given @c position must lie within the range of the underlying
0126   /// map.
0127   Material material(const Vector3& position) const {
0128     return Material(m_grid.atLocalBins(
0129         m_grid.localBinsFromLowerLeftEdge(m_transformPos(position))));
0130   }
0131 
0132   /// @brief Retrieve interpolated material at given position
0133   ///
0134   /// @param [in] position Global 3D position
0135   /// @return Material at the given position
0136   ///
0137   /// @pre The given @c position must lie within the range of the underlying
0138   /// map.
0139   Material getMaterial(const Vector3& position) const {
0140     return Material(m_grid.interpolate(m_transformPos(position)));
0141   }
0142 
0143   /// @brief Retrieve material cell for given position
0144   ///
0145   /// @param [in] position Global 3D position
0146   /// @return material cell containing the given global position
0147   ///
0148   /// @pre The given @c position must lie within the range of the underlying
0149   /// map.
0150   MaterialCell getMaterialCell(const Vector3& position) const {
0151     const auto& gridPosition = m_transformPos(position);
0152     std::size_t bin = m_grid.globalBinFromPosition(gridPosition);
0153     const auto& indices = m_grid.localBinsFromPosition(bin);
0154     const auto& lowerLeft = m_grid.lowerLeftBinEdge(indices);
0155     const auto& upperRight = m_grid.upperRightBinEdge(indices);
0156 
0157     // Loop through all corner points
0158     constexpr std::size_t nCorners = 1 << DIM_POS;
0159     std::array<Material::ParametersVector, nCorners> neighbors;
0160     const auto& cornerIndices = m_grid.closestPointsIndices(gridPosition);
0161 
0162     std::size_t i = 0;
0163     for (std::size_t index : cornerIndices) {
0164       neighbors.at(i++) = m_grid.at(index);
0165     }
0166 
0167     return MaterialCell(m_transformPos, lowerLeft, upperRight,
0168                         std::move(neighbors));
0169   }
0170 
0171   /// @brief Get the number of bins for all axes of the map
0172   ///
0173   /// @return Vector returning number of bins for all map axes
0174   std::vector<std::size_t> getNBins() const {
0175     auto nBinsArray = m_grid.numLocalBins();
0176     return std::vector<std::size_t>(nBinsArray.begin(), nBinsArray.end());
0177   }
0178 
0179   /// @brief Get the minimum value of all axes of the map
0180   ///
0181   /// @return Vector returning the minima of all map axes
0182   std::vector<double> getMin() const {
0183     auto minArray = m_grid.minPosition();
0184     return std::vector<double>(minArray.begin(), minArray.end());
0185   }
0186 
0187   /// @brief Get the maximum value of all axes of the map
0188   ///
0189   /// @return Vector returning the maxima of all map axes
0190   std::vector<double> getMax() const {
0191     auto maxArray = m_grid.maxPosition();
0192     return std::vector<double>(maxArray.begin(), maxArray.end());
0193   }
0194 
0195   /// @brief Check whether given 3D position is inside look-up domain
0196   ///
0197   /// @param [in] position Global 3D position
0198   /// @return @c true if position is inside the defined look-up grid,
0199   ///         otherwise @c false
0200   bool isInside(const Vector3& position) const {
0201     return m_grid.isInside(m_transformPos(position));
0202   }
0203 
0204   /// @brief Get a const reference on the underlying grid structure
0205   ///
0206   /// @return Grid reference
0207   const Grid_t& getGrid() const { return m_grid; }
0208 
0209  private:
0210   /// Geometric transformation applied to global 3D positions
0211   std::function<ActsVector<DIM_POS>(const Vector3&)> m_transformPos;
0212   /// Grid storing material values
0213   Grid_t m_grid;
0214 };
0215 
0216 /// @ingroup Material
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                                ActsMatrix<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                                ActsMatrix<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 }  // namespace Acts