Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:10:54

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