InterpolatedBFieldMap.hpp

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

// This file is part of the Acts project.
//
// Copyright (C) 2016-2018 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/MagneticField/MagneticFieldContext.hpp"
#include "Acts/Utilities/Definitions.hpp"
#include "Acts/Utilities/Interpolation.hpp"
#include "Acts/Utilities/detail/Grid.hpp"

#include <functional>
#include <optional>
#include <vector>

namespace Acts {

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

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

   public:
    FieldCell(std::function<ActsVectorD<DIM_POS>(const Vector3D&)> transformPos,
              std::array<double, DIM_POS> lowerLeft,
              std::array<double, DIM_POS> upperRight,
              std::array<Vector3D, N> fieldValues)
        : m_transformPos(std::move(transformPos)),
          m_lowerLeft(std::move(lowerLeft)),
          m_upperRight(std::move(upperRight)),
          m_fieldValues(std::move(fieldValues)) {}

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

    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<Vector3D, N> m_fieldValues;
  };

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

  Vector3D getField(const Vector3D& position) const {
    return m_transformBField(m_grid.interpolate(m_transformPos(position)),
                             position);
  }

  FieldCell getFieldCell(const Vector3D& position) const {
    const auto& gridPosition = m_transformPos(position);
    const auto& indices = m_grid.localBinsFromPosition(gridPosition);
    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<Vector3D, nCorners> neighbors;
    const auto& cornerIndices = m_grid.closestPointsIndices(gridPosition);

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

    return FieldCell(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;
  std::function<Vector3D(const FieldType&, const Vector3D&)> m_transformBField;
  Grid_t m_grid;
};

template <typename Mapper_t>
class InterpolatedBFieldMap final {
 public:
  struct Config {
    Config(Mapper_t m) : mapper(m) {}

    double scale = 1.;

    Mapper_t mapper;
  };

  struct Cache {
    Cache(std::reference_wrapper<const MagneticFieldContext> /*mcfg*/) {}

    std::optional<typename Mapper_t::FieldCell> fieldCell;
    bool initialized = false;
  };

  InterpolatedBFieldMap(Config config) : m_config(std::move(config)) {}

  Config getConfiguration() const { return m_config; }

  Vector3D getField(const Vector3D& position) const {
    return m_config.mapper.getField(position);
  }

  Vector3D getField(const Vector3D& position, Cache& cache) const {
    if (!cache.fieldCell || !(*cache.fieldCell).isInside(position)) {
      cache.fieldCell = getFieldCell(position);
    }
    return (*cache.fieldCell).getField(position);
  }

  Vector3D getFieldGradient(const Vector3D& position,
                            ActsMatrixD<3, 3>& /*derivative*/) const {
    return m_config.mapper.getField(position);
  }

  Vector3D getFieldGradient(const Vector3D& position,
                            ActsMatrixD<3, 3>& /*derivative*/,
                            Cache& /*cache*/) const {
    return m_config.mapper.getField(position);
  }

  double getScale() const { return m_config.scale; }

  Mapper_t getMapper() const { return m_config.mapper; }

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

  void setConfiguration(const Config& config) { m_config = config; }

 private:
  typename Mapper_t::FieldCell getFieldCell(const Vector3D& position) const {
    return m_config.mapper.getFieldCell(position);
  }

  Config m_config;
};

}  // namespace Acts