ECCE @ EIC 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 
11 #include <functional>
12 #include <optional>
18 namespace Acts {
19 
25 template <typename G>
27  public:
28  using Grid_t = G;
29  static constexpr size_t DIM_POS = Grid_t::DIM;
30 
35  struct MaterialCell {
37  static constexpr unsigned int N = 1 << DIM_POS;
38 
52  std::function<ActsVectorD<DIM_POS>(const Vector3D&)> transformPos,
53  std::array<double, DIM_POS> lowerLeft,
54  std::array<double, DIM_POS> upperRight,
55  std::array<ActsVectorF<5>, N> materialValues)
56  : m_transformPos(std::move(transformPos)),
57  m_lowerLeft(std::move(lowerLeft)),
58  m_upperRight(std::move(upperRight)),
59  m_materialValues(std::move(materialValues)) {}
60 
68  // defined in Interpolation.hpp
71  }
72 
78  bool isInside(const Vector3D& position) const {
79  const auto& gridCoordinates = m_transformPos(position);
80  for (unsigned int i = 0; i < DIM_POS; ++i) {
81  if (gridCoordinates[i] < m_lowerLeft.at(i) ||
82  gridCoordinates[i] >= m_upperRight.at(i)) {
83  return false;
84  }
85  }
86  return true;
87  }
88 
89  private:
91  std::function<ActsVectorD<DIM_POS>(const Vector3D&)> m_transformPos;
92 
94  std::array<double, DIM_POS> m_lowerLeft;
95 
97  std::array<double, DIM_POS> m_upperRight;
98 
103  std::array<ActsVectorF<5>, N> m_materialValues;
104  };
105 
112  std::function<ActsVectorD<DIM_POS>(const Vector3D&)> transformPos,
113  Grid_t grid)
114  : m_transformPos(std::move(transformPos)), m_grid(std::move(grid)) {}
115 
124  return Material(m_grid.atLocalBins(
125  m_grid.localBinsFromLowerLeftEdge(m_transformPos(position))));
126  }
127 
136  return Material(m_grid.interpolate(m_transformPos(position)));
137  }
138 
147  const auto& gridPosition = m_transformPos(position);
148  size_t bin = m_grid.globalBinFromPosition(gridPosition);
149  const auto& indices = m_grid.localBinsFromPosition(bin);
150  const auto& lowerLeft = m_grid.lowerLeftBinEdge(indices);
151  const auto& upperRight = m_grid.upperRightBinEdge(indices);
152 
153  // Loop through all corner points
154  constexpr size_t nCorners = 1 << DIM_POS;
155  std::array<ActsVectorF<5>, nCorners> neighbors;
156  const auto& cornerIndices = m_grid.closestPointsIndices(gridPosition);
157 
158  size_t i = 0;
159  for (size_t index : cornerIndices) {
160  neighbors.at(i++) = m_grid.at(index);
161  }
162 
163  return MaterialCell(m_transformPos, lowerLeft, upperRight,
164  std::move(neighbors));
165  }
166 
170  std::vector<size_t> getNBins() const {
171  auto nBinsArray = m_grid.numLocalBins();
172  return std::vector<size_t>(nBinsArray.begin(), nBinsArray.end());
173  }
174 
178  std::vector<double> getMin() const {
179  auto minArray = m_grid.minPosition();
180  return std::vector<double>(minArray.begin(), minArray.end());
181  }
182 
186  std::vector<double> getMax() const {
187  auto maxArray = m_grid.maxPosition();
188  return std::vector<double>(maxArray.begin(), maxArray.end());
189  }
190 
196  bool isInside(const Vector3D& position) const {
197  return m_grid.isInside(m_transformPos(position));
198  }
199 
203  const Grid_t& getGrid() const { return m_grid; }
204 
205  private:
207  std::function<ActsVectorD<DIM_POS>(const Vector3D&)> m_transformPos;
210 };
211 
234 template <typename Mapper_t>
236  public:
238  struct Cache {
240  std::optional<typename Mapper_t::MaterialCell> matCell;
242  bool initialized = false;
243  };
244 
248  InterpolatedMaterialMap(Mapper_t&& mapper) : m_mapper(std::move(mapper)) {}
249 
254  InterpolatedMaterialMap(Mapper_t&& mapper, BinUtility bu)
255  : m_mapper(std::move(mapper)), m_binUtility(bu) {}
256 
262  const Material material(const Vector3D& position) const {
263  return m_mapper.material(position);
264  }
265 
272  return m_mapper.getMaterial(position);
273  }
274 
282  Material getMaterial(const Vector3D& position, Cache& cache) const {
283  if (!cache.initialized || !(*cache.matCell).isInside(position)) {
284  cache.matCell = getMaterialCell(position);
285  cache.initialized = true;
286  }
287  return (*cache.matCell).getMaterial(position);
288  }
289 
299  ActsMatrixD<5, 5>& /*derivative*/) const {
300  return m_mapper.getMaterial(position);
301  }
302 
314  ActsMatrixD<5, 5>& /*derivative*/,
315  Cache& /*cache*/) const {
316  return m_mapper.getMaterial(position);
317  }
318 
322  const Mapper_t& getMapper() const { return m_mapper; }
323 
328  bool isInside(const Vector3D& position) const {
329  return m_mapper.isInside(position);
330  }
331 
333  const BinUtility& binUtility() const { return m_binUtility; }
334 
338  std::ostream& toStream(std::ostream& sl) const {
339  sl << "Acts::InterpolatedMaterialMap : " << std::endl;
340  sl << " - Number of Material bins [0,1] : " << m_binUtility.max(0) + 1
341  << " / " << m_binUtility.max(1) + 1 << std::endl;
342  sl << " - Parse full update material : " << std::endl; //
343  }
344 
345  private:
353  typename Mapper_t::MaterialCell getMaterialCell(
354  const Vector3D& position) const {
355  return m_mapper.getMaterialCell(position);
356  }
357 
363  Mapper_t m_mapper;
364 
366 };
367 } // namespace Acts