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/Definitions/Algebra.hpp"
#include "Acts/MagneticField/MagneticFieldContext.hpp"
#include "Acts/MagneticField/MagneticFieldError.hpp"
#include "Acts/MagneticField/MagneticFieldProvider.hpp"
#include "Acts/Utilities/Interpolation.hpp"
#include "Acts/Utilities/Result.hpp"
#include "Acts/Utilities/detail/Grid.hpp"

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

namespace Acts {

class InterpolatedMagneticField : public MagneticFieldProvider {
 public:
  virtual std::vector<size_t> getNBins() const = 0;

  virtual std::vector<double> getMin() const = 0;

  virtual std::vector<double> getMax() const = 0;

  virtual bool isInside(const Vector3& position) const = 0;

  virtual Vector3 getFieldUnchecked(const Vector3& position) const = 0;
};

template <typename grid_t>
class InterpolatedBFieldMap : public InterpolatedMagneticField {
 public:
  using Grid = grid_t;
  using FieldType = typename Grid::value_type;
  static constexpr size_t DIM_POS = Grid::DIM;

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

   public:
    FieldCell(std::array<double, DIM_POS> lowerLeft,
              std::array<double, DIM_POS> upperRight,
              std::array<Vector3, N> fieldValues)
        : m_lowerLeft(std::move(lowerLeft)),
          m_upperRight(std::move(upperRight)),
          m_fieldValues(std::move(fieldValues)) {}

    Vector3 getField(const ActsVector<DIM_POS>& position) const {
      // defined in Interpolation.hpp
      return interpolate(position, m_lowerLeft, m_upperRight, m_fieldValues);
    }

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

   private:
    std::array<double, DIM_POS> m_lowerLeft;

    std::array<double, DIM_POS> m_upperRight;

    std::array<Vector3, N> m_fieldValues;
  };

  struct Cache {
    Cache(const MagneticFieldContext& mctx) { (void)mctx; }

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

  struct Config {
    std::function<ActsVector<DIM_POS>(const Vector3&)> transformPos;

    std::function<Vector3(const FieldType&, const Vector3&)> transformBField;

    Grid grid;

    double scale = 1.;
  };

  InterpolatedBFieldMap(Config cfg) : m_cfg{std::move(cfg)} {
    typename Grid::index_t minBin{};
    minBin.fill(1);
    m_lowerLeft = m_cfg.grid.lowerLeftBinEdge(minBin);
    m_upperRight = m_cfg.grid.lowerLeftBinEdge(m_cfg.grid.numLocalBins());
  }

  Result<FieldCell> getFieldCell(const Vector3& position) const {
    const auto& gridPosition = m_cfg.transformPos(position);
    const auto& indices = m_cfg.grid.localBinsFromPosition(gridPosition);
    const auto& lowerLeft = m_cfg.grid.lowerLeftBinEdge(indices);
    const auto& upperRight = m_cfg.grid.upperRightBinEdge(indices);

    // loop through all corner points
    constexpr size_t nCorners = 1 << DIM_POS;
    std::array<Vector3, nCorners> neighbors;
    const auto& cornerIndices = m_cfg.grid.closestPointsIndices(gridPosition);

    if (!isInsideLocal(gridPosition)) {
      return MagneticFieldError::OutOfBounds;
    }

    size_t i = 0;
    for (size_t index : cornerIndices) {
      neighbors.at(i++) = m_cfg.transformBField(m_cfg.grid.at(index), position);
    }

    assert(i == nCorners);

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

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

  std::vector<double> getMin() const final {
    return std::vector<double>(m_lowerLeft.begin(), m_lowerLeft.end());
  }

  std::vector<double> getMax() const final {
    return std::vector<double>(m_upperRight.begin(), m_upperRight.end());
  }

  bool isInside(const Vector3& position) const final {
    return isInsideLocal(m_cfg.transformPos(position));
  }

  bool isInsideLocal(const ActsVector<DIM_POS>& gridPosition) const {
    for (unsigned int i = 0; i < DIM_POS; ++i) {
      if (gridPosition[i] < m_lowerLeft[i] ||
          gridPosition[i] >= m_upperRight[i]) {
        return false;
      }
    }
    return true;
  }

  const Grid& getGrid() const { return m_cfg.grid; }

  MagneticFieldProvider::Cache makeCache(
      const MagneticFieldContext& mctx) const final {
    return MagneticFieldProvider::Cache::make<Cache>(mctx);
  }

  Result<Vector3> getField(const Vector3& position) const {
    const auto gridPosition = m_cfg.transformPos(position);
    if (!isInsideLocal(gridPosition)) {
      return Result<Vector3>::failure(MagneticFieldError::OutOfBounds);
    }

    return Result<Vector3>::success(
        m_cfg.transformBField(m_cfg.grid.interpolate(gridPosition), position));
  }

  Vector3 getFieldUnchecked(const Vector3& position) const final {
    const auto gridPosition = m_cfg.transformPos(position);
    return m_cfg.transformBField(m_cfg.grid.interpolate(gridPosition),
                                 position);
  }

  Result<Vector3> getField(const Vector3& position,
                           MagneticFieldProvider::Cache& cache) const final {
    Cache& lcache = cache.get<Cache>();
    const auto gridPosition = m_cfg.transformPos(position);
    if (!lcache.fieldCell || !(*lcache.fieldCell).isInside(gridPosition)) {
      auto res = getFieldCell(position);
      if (!res.ok()) {
        return Result<Vector3>::failure(res.error());
      }
      lcache.fieldCell = *res;
    }
    return Result<Vector3>::success((*lcache.fieldCell).getField(gridPosition));
  }

  Result<Vector3> getFieldGradient(
      const Vector3& position, ActsMatrix<3, 3>& derivative,
      MagneticFieldProvider::Cache& cache) const final {
    (void)derivative;
    return getField(position, cache);
  }

 private:
  Config m_cfg;

  typename Grid::point_t m_lowerLeft;
  typename Grid::point_t m_upperRight;
};

}  // namespace Acts