Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
InterpolatedMaterialMap.hpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file InterpolatedMaterialMap.hpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019-2020 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #pragma once
10 
16 
17 #include <functional>
18 #include <optional>
19 
20 namespace Acts {
21 
27 template <typename G>
29  public:
30  using Grid_t = G;
31  static constexpr size_t DIM_POS = Grid_t::DIM;
32 
37  struct MaterialCell {
39  static constexpr unsigned int N = 1 << DIM_POS;
40 
54  std::function<ActsVector<DIM_POS>(const Vector3&)> transformPos,
55  std::array<double, DIM_POS> lowerLeft,
56  std::array<double, DIM_POS> upperRight,
57  std::array<Material::ParametersVector, N> materialValues)
58  : m_transformPos(std::move(transformPos)),
59  m_lowerLeft(std::move(lowerLeft)),
60  m_upperRight(std::move(upperRight)),
61  m_materialValues(std::move(materialValues)) {}
62 
70  // defined in Interpolation.hpp
73  }
74 
80  bool isInside(const Vector3& position) const {
81  const auto& gridCoordinates = m_transformPos(position);
82  for (unsigned int i = 0; i < DIM_POS; ++i) {
83  if (gridCoordinates[i] < m_lowerLeft.at(i) ||
84  gridCoordinates[i] >= m_upperRight.at(i)) {
85  return false;
86  }
87  }
88  return true;
89  }
90 
91  private:
93  std::function<ActsVector<DIM_POS>(const Vector3&)> m_transformPos;
94 
96  std::array<double, DIM_POS> m_lowerLeft;
97 
99  std::array<double, DIM_POS> m_upperRight;
100 
105  std::array<Material::ParametersVector, N> m_materialValues;
106  };
107 
114  std::function<ActsVector<DIM_POS>(const Vector3&)> transformPos,
115  Grid_t grid)
116  : m_transformPos(std::move(transformPos)), m_grid(std::move(grid)) {}
117 
126  return Material(m_grid.atLocalBins(
127  m_grid.localBinsFromLowerLeftEdge(m_transformPos(position))));
128  }
129 
138  return Material(m_grid.interpolate(m_transformPos(position)));
139  }
140 
149  const auto& gridPosition = m_transformPos(position);
150  size_t bin = m_grid.globalBinFromPosition(gridPosition);
151  const auto& indices = m_grid.localBinsFromPosition(bin);
152  const auto& lowerLeft = m_grid.lowerLeftBinEdge(indices);
153  const auto& upperRight = m_grid.upperRightBinEdge(indices);
154 
155  // Loop through all corner points
156  constexpr size_t nCorners = 1 << DIM_POS;
157  std::array<Material::ParametersVector, nCorners> neighbors;
158  const auto& cornerIndices = m_grid.closestPointsIndices(gridPosition);
159 
160  size_t i = 0;
161  for (size_t index : cornerIndices) {
162  neighbors.at(i++) = m_grid.at(index);
163  }
164 
165  return MaterialCell(m_transformPos, lowerLeft, upperRight,
166  std::move(neighbors));
167  }
168 
172  std::vector<size_t> getNBins() const {
173  auto nBinsArray = m_grid.numLocalBins();
174  return std::vector<size_t>(nBinsArray.begin(), nBinsArray.end());
175  }
176 
180  std::vector<double> getMin() const {
181  auto minArray = m_grid.minPosition();
182  return std::vector<double>(minArray.begin(), minArray.end());
183  }
184 
188  std::vector<double> getMax() const {
189  auto maxArray = m_grid.maxPosition();
190  return std::vector<double>(maxArray.begin(), maxArray.end());
191  }
192 
198  bool isInside(const Vector3& position) const {
199  return m_grid.isInside(m_transformPos(position));
200  }
201 
205  const Grid_t& getGrid() const { return m_grid; }
206 
207  private:
209  std::function<ActsVector<DIM_POS>(const Vector3&)> m_transformPos;
212 };
213 
236 template <typename Mapper_t>
238  public:
240  struct Cache {
242  std::optional<typename Mapper_t::MaterialCell> matCell;
244  bool initialized = false;
245  };
246 
251 
257  : m_mapper(std::move(mapper)), m_binUtility(std::move(bu)) {}
258 
264  const Material material(const Vector3& position) const override {
265  return m_mapper.material(position);
266  }
267 
274  return m_mapper.getMaterial(position);
275  }
276 
284  Material getMaterial(const Vector3& position, Cache& cache) const {
285  if (!cache.initialized || !(*cache.matCell).isInside(position)) {
286  cache.matCell = getMaterialCell(position);
287  cache.initialized = true;
288  }
289  return (*cache.matCell).getMaterial(position);
290  }
291 
301  ActsMatrix<5, 5>& derivative) const {
302  (void)derivative;
303  return m_mapper.getMaterial(position);
304  }
305 
317  ActsMatrix<5, 5>& derivative,
318  Cache& cache) const {
319  (void)derivative;
320  (void)cache;
321  return m_mapper.getMaterial(position);
322  }
323 
327  const Mapper_t& getMapper() const { return m_mapper; }
328 
333  bool isInside(const Vector3& position) const {
334  return m_mapper.isInside(position);
335  }
336 
338  const BinUtility& binUtility() const { return m_binUtility; }
339 
343  std::ostream& toStream(std::ostream& sl) const override {
344  sl << "Acts::InterpolatedMaterialMap : " << std::endl;
345  sl << " - Number of Material bins [0,1] : " << m_binUtility.max(0) + 1
346  << " / " << m_binUtility.max(1) + 1 << std::endl;
347  sl << " - Parse full update material : " << std::endl;
348  return sl;
349  }
350 
351  private:
359  typename Mapper_t::MaterialCell getMaterialCell(
360  const Vector3& position) const {
361  return m_mapper.getMaterialCell(position);
362  }
363 
369  Mapper_t m_mapper;
370 
372 };
373 } // namespace Acts