InterpolatedMaterialMap.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file InterpolatedMaterialMap.hpp

// This file is part of the Acts project.
//
// Copyright (C) 2019-2020 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once

#include "Acts/Material/IVolumeMaterial.hpp"
#include "Acts/Material/Material.hpp"
#include "Acts/Utilities/BinUtility.hpp"
#include "Acts/Utilities/Definitions.hpp"
#include "Acts/Utilities/Interpolation.hpp"

#include <functional>
#include <optional>

namespace Acts {

template <typename G>
struct MaterialMapper {
 public:
  using Grid_t = G;
  static constexpr size_t DIM_POS = Grid_t::DIM;

  struct MaterialCell {
    static constexpr unsigned int N = 1 << DIM_POS;

    MaterialCell(
        std::function<ActsVectorD<DIM_POS>(const Vector3D&)> transformPos,
        std::array<double, DIM_POS> lowerLeft,
        std::array<double, DIM_POS> upperRight,
        std::array<ActsVectorF<5>, N> materialValues)
        : m_transformPos(std::move(transformPos)),
          m_lowerLeft(std::move(lowerLeft)),
          m_upperRight(std::move(upperRight)),
          m_materialValues(std::move(materialValues)) {}

    Material getMaterial(const Vector3D& position) const {
      // defined in Interpolation.hpp
      return Material(interpolate(m_transformPos(position), m_lowerLeft,
                                  m_upperRight, m_materialValues));
    }

    bool isInside(const Vector3D& position) const {
      const auto& gridCoordinates = m_transformPos(position);
      for (unsigned int i = 0; i < DIM_POS; ++i) {
        if (gridCoordinates[i] < m_lowerLeft.at(i) ||
            gridCoordinates[i] >= m_upperRight.at(i)) {
          return false;
        }
      }
      return true;
    }

   private:
    std::function<ActsVectorD<DIM_POS>(const Vector3D&)> m_transformPos;

    std::array<double, DIM_POS> m_lowerLeft;

    std::array<double, DIM_POS> m_upperRight;

    std::array<ActsVectorF<5>, N> m_materialValues;
  };

  MaterialMapper(
      std::function<ActsVectorD<DIM_POS>(const Vector3D&)> transformPos,
      Grid_t grid)
      : m_transformPos(std::move(transformPos)), m_grid(std::move(grid)) {}

  Material material(const Vector3D& position) const {
    return Material(m_grid.atLocalBins(
        m_grid.localBinsFromLowerLeftEdge(m_transformPos(position))));
  }

  Material getMaterial(const Vector3D& position) const {
    return Material(m_grid.interpolate(m_transformPos(position)));
  }

  MaterialCell getMaterialCell(const Vector3D& position) const {
    const auto& gridPosition = m_transformPos(position);
    size_t bin = m_grid.globalBinFromPosition(gridPosition);
    const auto& indices = m_grid.localBinsFromPosition(bin);
    const auto& lowerLeft = m_grid.lowerLeftBinEdge(indices);
    const auto& upperRight = m_grid.upperRightBinEdge(indices);

    // Loop through all corner points
    constexpr size_t nCorners = 1 << DIM_POS;
    std::array<ActsVectorF<5>, nCorners> neighbors;
    const auto& cornerIndices = m_grid.closestPointsIndices(gridPosition);

    size_t i = 0;
    for (size_t index : cornerIndices) {
      neighbors.at(i++) = m_grid.at(index);
    }

    return MaterialCell(m_transformPos, lowerLeft, upperRight,
                        std::move(neighbors));
  }

  std::vector<size_t> getNBins() const {
    auto nBinsArray = m_grid.numLocalBins();
    return std::vector<size_t>(nBinsArray.begin(), nBinsArray.end());
  }

  std::vector<double> getMin() const {
    auto minArray = m_grid.minPosition();
    return std::vector<double>(minArray.begin(), minArray.end());
  }

  std::vector<double> getMax() const {
    auto maxArray = m_grid.maxPosition();
    return std::vector<double>(maxArray.begin(), maxArray.end());
  }

  bool isInside(const Vector3D& position) const {
    return m_grid.isInside(m_transformPos(position));
  }

  const Grid_t& getGrid() const { return m_grid; }

 private:
  std::function<ActsVectorD<DIM_POS>(const Vector3D&)> m_transformPos;
  Grid_t m_grid;
};

template <typename Mapper_t>
class InterpolatedMaterialMap : public IVolumeMaterial {
 public:
  struct Cache {
    std::optional<typename Mapper_t::MaterialCell> matCell;
    bool initialized = false;
  };

  InterpolatedMaterialMap(Mapper_t&& mapper) : m_mapper(std::move(mapper)) {}

  InterpolatedMaterialMap(Mapper_t&& mapper, BinUtility bu)
      : m_mapper(std::move(mapper)), m_binUtility(bu) {}

  const Material material(const Vector3D& position) const {
    return m_mapper.material(position);
  }

  Material getMaterial(const Vector3D& position) const {
    return m_mapper.getMaterial(position);
  }

  Material getMaterial(const Vector3D& position, Cache& cache) const {
    if (!cache.initialized || !(*cache.matCell).isInside(position)) {
      cache.matCell = getMaterialCell(position);
      cache.initialized = true;
    }
    return (*cache.matCell).getMaterial(position);
  }

  Material getMaterialGradient(const Vector3D& position,
                               ActsMatrixD<5, 5>& /*derivative*/) const {
    return m_mapper.getMaterial(position);
  }

  Material getMaterialGradient(const Vector3D& position,
                               ActsMatrixD<5, 5>& /*derivative*/,
                               Cache& /*cache*/) const {
    return m_mapper.getMaterial(position);
  }

  const Mapper_t& getMapper() const { return m_mapper; }

  bool isInside(const Vector3D& position) const {
    return m_mapper.isInside(position);
  }

  const BinUtility& binUtility() const { return m_binUtility; }

  std::ostream& toStream(std::ostream& sl) const {
    sl << "Acts::InterpolatedMaterialMap : " << std::endl;
    sl << "   - Number of Material bins [0,1] : " << m_binUtility.max(0) + 1
       << " / " << m_binUtility.max(1) + 1 << std::endl;
    sl << "   - Parse full update material    : " << std::endl;
    return sl;
  }

 private:
  typename Mapper_t::MaterialCell getMaterialCell(
      const Vector3D& position) const {
    return m_mapper.getMaterialCell(position);
  }

  Mapper_t m_mapper;

  BinUtility m_binUtility{};
};
}  // namespace Acts