refactor!: B field access returns Result (#825)

Previously, we didn't have a mechanism to communicate that a field query failed. This is problematic, because it's not clear what to do when for instance our interpolated field goes out of bounds (see #784). This PR makes `MagneticFieldProvider` require the `getField` and `getFieldGradient` methods to return `Result<Vector3>` instead of the bare vector. This changes the public interface slightly, and all clients need to be adjusted.

Additionally, the `InterpolatedBFieldMap` now checks each time a field cell is created or if the field is queried directly, whether the query position is out of bounds of the map, and will return a new `MagneticFieldError::OutOfBounds` if that happens.

BREAKING CHANGE:
- The signature of field query methods in `MagneticFieldProvider` changes from
   ```cpp
   virtual Vector3 getField(const Vector3& position, Cache& cache) const = 0;
   virtual Vector3 getFieldGradient(const Vector3& position, ActsMatrix<3, 3>& derivative, Cache& cache) const = 0;
   ```
   to
   ```cpp
   virtual Result<Vector3> getField(const Vector3& position, Cache& cache) const = 0;
   virtual Result<Vector3> getFieldGradient(const Vector3& position, ActsMatrix<3, 3>& derivative, Cache& cache) const = 0;
   ```
- `InterpolatedBFieldMap::getMin` and `InterpolatedBFieldMap::getMax` now return the extent of the valid **interpolation domain**, rather than the raw grid extent.

Core/include/Acts/MagneticField/ConstantBField.hpp

@@ -38,9 +38,10 @@ class ConstantBField final : public MagneticFieldProvider {
   ///
   /// @note The @p position is ignored and only kept as argument to provide
   ///       a consistent interface with other magnetic field services.
-  Vector3 getField(const Vector3& /*position*/,
-                   MagneticFieldProvider::Cache& /*cache*/) const override {
-    return m_BField;
+  Result<Vector3> getField(
+      const Vector3& /*position*/,
+      MagneticFieldProvider::Cache& /*cache*/) const override {
+    return Result<Vector3>::success(m_BField);
   }
 
   /// @copydoc MagneticFieldProvider::getFieldGradient(const
@@ -50,10 +51,10 @@ class ConstantBField final : public MagneticFieldProvider {
   ///       a consistent interface with other magnetic field services.
   /// @note currently the derivative is not calculated
   /// @todo return derivative
-  Vector3 getFieldGradient(
+  Result<Vector3> getFieldGradient(
       const Vector3& /*position*/, ActsMatrix<3, 3>& /*derivative*/,
       MagneticFieldProvider::Cache& /*cache*/) const override {
-    return m_BField;
+    return Result<Vector3>::success(m_BField);
   }
 
   /// @copydoc MagneticFieldProvider::makeCache(const MagneticFieldContext&)

Core/include/Acts/MagneticField/InterpolatedBFieldMap.hpp

@@ -10,8 +10,10 @@
 
 #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>
@@ -84,8 +86,8 @@ struct InterpolatedBFieldMapper {
     bool isInside(const Vector3& 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)) {
+        if (gridCoordinates[i] < m_lowerLeft[i] ||
+            gridCoordinates[i] >= m_upperRight[i]) {
           return false;
         }
       }
@@ -123,7 +125,12 @@ struct InterpolatedBFieldMapper {
       Grid_t grid)
       : m_transformPos(std::move(transformPos)),
         m_transformBField(std::move(transformBField)),
-        m_grid(std::move(grid)) {}
+        m_grid(std::move(grid)) {
+    typename Grid_t::index_t minBin{};
+    minBin.fill(1);
+    m_lowerLeft = m_grid.lowerLeftBinEdge(minBin);
+    m_upperRight = m_grid.lowerLeftBinEdge(m_grid.numLocalBins());
+  }
 
   /// @brief retrieve field at given position
   ///
@@ -132,9 +139,13 @@ struct InterpolatedBFieldMapper {
   ///
   /// @pre The given @c position must lie within the range of the underlying
   ///      magnetic field map.
-  Vector3 getField(const Vector3& position) const {
-    return m_transformBField(m_grid.interpolate(m_transformPos(position)),
-                             position);
+  Result<Vector3> getField(const Vector3& position) const {
+    if (!isInside(position)) {
+      return Result<Vector3>::failure(MagneticFieldError::OutOfBounds);
+    }
+
+    return Result<Vector3>::success(m_transformBField(
+        m_grid.interpolate(m_transformPos(position)), position));
   }
 
   /// @brief retrieve field cell for given position
@@ -144,17 +155,21 @@ struct InterpolatedBFieldMapper {
   ///
   /// @pre The given @c position must lie within the range of the underlying
   ///      magnetic field map.
-  FieldCell getFieldCell(const Vector3& position) const {
+  Result<FieldCell> getFieldCell(const Vector3& 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<Vector3, nCorners> neighbors;
+    std::array<Vector3, nCorners> neighbors{Vector3{0.0, 0.0, 0.0}};
     const auto& cornerIndices = m_grid.closestPointsIndices(gridPosition);
 
+    if (!isInside(position)) {
+      return MagneticFieldError::OutOfBounds;
+    }
+
     size_t i = 0;
     for (size_t index : cornerIndices) {
       neighbors.at(i++) = m_transformBField(m_grid.at(index), position);
@@ -176,16 +191,14 @@ struct InterpolatedBFieldMapper {
   ///
   /// @return vector returning the minima of all field map axes
   std::vector<double> getMin() const {
-    auto minArray = m_grid.minPosition();
-    return std::vector<double>(minArray.begin(), minArray.end());
+    return std::vector<double>(m_lowerLeft.begin(), m_lowerLeft.end());
   }
 
   /// @brief get the maximum value of all axes of the field map
   ///
   /// @return vector returning the maxima of all field map axes
   std::vector<double> getMax() const {
-    auto maxArray = m_grid.maxPosition();
-    return std::vector<double>(maxArray.begin(), maxArray.end());
+    return std::vector<double>(m_upperRight.begin(), m_upperRight.end());
   }
 
   /// @brief check whether given 3D position is inside look-up domain
@@ -194,7 +207,14 @@ struct InterpolatedBFieldMapper {
   /// @return @c true if position is inside the defined look-up grid,
   ///         otherwise @c false
   bool isInside(const Vector3& position) const {
-    return m_grid.isInside(m_transformPos(position));
+    const auto& gridPosition = m_transformPos(position);
+    for (unsigned int i = 0; i < DIM_POS; ++i) {
+      if (gridPosition[i] < m_lowerLeft[i] ||
+          gridPosition[i] >= m_upperRight[i]) {
+        return false;
+      }
+    }
+    return true;
   }
 
   /// @brief Get a const reference on the underlying grid structure
@@ -211,6 +231,9 @@ struct InterpolatedBFieldMapper {
   std::function<Vector3(const FieldType&, const Vector3&)> m_transformBField;
   /// grid storing magnetic field values
   Grid_t m_grid;
+
+  typename Grid_t::point_t m_lowerLeft;
+  typename Grid_t::point_t m_upperRight;
 };
 
 /// @ingroup MagneticField
@@ -300,22 +323,25 @@ class InterpolatedBFieldMap final : public MagneticFieldProvider {
     return MagneticFieldProvider::Cache::make<Cache>(mctx);
   }
 
-  /// @brief Get the B field at a position
-  ///
-  /// @param position The position to query at
-  Vector3 getField(const Vector3& position) const {
+  /// @copydoc MagneticFieldProvider::getField(const Vector3&)
+  Result<Vector3> getField(const Vector3& position) const {
     return m_config.mapper.getField(position);
   }
 
   /// @copydoc MagneticFieldProvider::getField(const
   /// Vector3&,MagneticFieldProvider::Cache&)
-  Vector3 getField(const Vector3& position,
-                   MagneticFieldProvider::Cache& gcache) const override {
+  Result<Vector3> getField(
+      const Vector3& position,
+      MagneticFieldProvider::Cache& gcache) const override {
     Cache& cache = gcache.get<Cache>();
     if (!cache.fieldCell || !(*cache.fieldCell).isInside(position)) {
-      cache.fieldCell = getFieldCell(position);
+      auto res = getFieldCell(position);
+      if (!res.ok()) {
+        return Result<Vector3>::failure(res.error());
+      }
+      cache.fieldCell = *res;
     }
-    return (*cache.fieldCell).getField(position);
+    return Result<Vector3>::success((*cache.fieldCell).getField(position));
   }
 
   /// @copydoc MagneticFieldProvider::getFieldGradient(const
@@ -324,10 +350,10 @@ class InterpolatedBFieldMap final : public MagneticFieldProvider {
   /// @note currently the derivative is not calculated
   /// @note Cache is not used currently
   /// @todo return derivative
-  Vector3 getFieldGradient(
+  Result<Vector3> getFieldGradient(
       const Vector3& position, ActsMatrix<3, 3>& /*derivative*/,
-      MagneticFieldProvider::Cache& /*cache*/) const override {
-    return m_config.mapper.getField(position);
+      MagneticFieldProvider::Cache& cache) const override {
+    return getField(position, cache);
   }
 
  private:
@@ -338,7 +364,8 @@ class InterpolatedBFieldMap final : public MagneticFieldProvider {
   ///
   /// @pre The given @c position must lie within the range of the underlying
   ///      magnetic field map.
-  typename Mapper_t::FieldCell getFieldCell(const Vector3& position) const {
+  Result<typename Mapper_t::FieldCell> getFieldCell(
+      const Vector3& position) const {
     return m_config.mapper.getFieldCell(position);
   }
 

Core/include/Acts/MagneticField/MagneticFieldError.hpp

@@ -0,0 +1,27 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2021 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 <system_error>
+
+namespace Acts {
+
+enum class MagneticFieldError {
+  OutOfBounds = 1,
+};
+
+std::error_code make_error_code(Acts::MagneticFieldError e);
+
+}  // namespace Acts
+
+namespace std {
+// register with STL
+template <>
+struct is_error_code_enum<Acts::MagneticFieldError> : std::true_type {};
+}  // namespace std

Core/include/Acts/MagneticField/MagneticFieldProvider.hpp

@@ -11,6 +11,7 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/MagneticField/MagneticFieldContext.hpp"
 #include "Acts/MagneticField/detail/SmallObjectCache.hpp"
+#include "Acts/Utilities/Result.hpp"
 
 #include <array>
 #include <memory>
@@ -33,17 +34,18 @@ class MagneticFieldProvider {
   /// @param [in] position global 3D position
   ///
   /// @return magnetic field vector at given position
-  virtual Vector3 getField(const Vector3& position, Cache& cache) const = 0;
+  virtual Result<Vector3> getField(const Vector3& position,
+                                   Cache& cache) const = 0;
 
   /// @brief retrieve magnetic field value & its gradient
   ///
   /// @param [in]  position   global 3D position
   /// @param [out] derivative gradient of magnetic field vector as (3x3) matrix
   /// @param [in,out] cache Cache object. Contains field cell used for
   /// @return magnetic field vector
-  virtual Vector3 getFieldGradient(const Vector3& position,
-                                   ActsMatrix<3, 3>& derivative,
-                                   Cache& cache) const = 0;
+  virtual Result<Vector3> getFieldGradient(const Vector3& position,
+                                           ActsMatrix<3, 3>& derivative,
+                                           Cache& cache) const = 0;
 };
 
 }  // namespace Acts

Core/include/Acts/MagneticField/NullBField.hpp

@@ -30,9 +30,10 @@ class NullBField final : public MagneticFieldProvider {
   ///
   /// @note The @p position is ignored and only kept as argument to provide
   ///       a consistent interface with other magnetic field services.
-  Vector3 getField(const Vector3& /*position*/,
-                   MagneticFieldProvider::Cache& /*cache*/) const override {
-    return m_BField;
+  Result<Vector3> getField(
+      const Vector3& /*position*/,
+      MagneticFieldProvider::Cache& /*cache*/) const override {
+    return Result<Vector3>::success(m_BField);
   }
 
   /// @copydoc MagneticFieldProvider::getFieldGradient(const
@@ -42,10 +43,10 @@ class NullBField final : public MagneticFieldProvider {
   ///       a consistent interface with other magnetic field services.
   /// @note currently the derivative is not calculated
   /// @todo return derivative
-  Vector3 getFieldGradient(
+  Result<Vector3> getFieldGradient(
       const Vector3& /*position*/, ActsMatrix<3, 3>& /*derivative*/,
       MagneticFieldProvider::Cache& /*cache*/) const override {
-    return m_BField;
+    return Result<Vector3>::success(m_BField);
   }
 
   /// @copydoc MagneticFieldProvider::makeCache(const MagneticFieldContext&)

Core/include/Acts/MagneticField/SharedBField.hpp

@@ -33,32 +33,18 @@ class SharedBField final : public MagneticFieldProvider {
   /// @tparam bField is the shared BField to be stored
   SharedBField(std::shared_ptr<const BField> bField) : m_bField(bField) {}
 
-  /// @brief Get the B field at a position
-  ///
-  /// @param position The position to query at
-  Vector3 getField(const Vector3& position) const {
-    return m_bField->getField(position);
-  }
-
   /// @copydoc MagneticFieldProvider::getField(const
   /// Vector3&,MagneticFieldProvider::Cache&)
-  Vector3 getField(const Vector3& position,
-                   MagneticFieldProvider::Cache& cache) const override {
+  Result<Vector3> getField(const Vector3& position,
+                           MagneticFieldProvider::Cache& cache) const override {
     return m_bField->getField(position, cache);
   }
 
-  /// @copydoc MagneticFieldProvider::getFieldGradient(const
-  /// Vector3&,ActsMatrix<3,3>&)
-  Vector3 getFieldGradient(const Vector3& position,
-                           ActsMatrix<3, 3>& derivative) const {
-    return m_bField->getFieldGradient(position, derivative);
-  }
-
   /// @copydoc MagneticFieldProvider::getFieldGradient(const
   /// Vector3&,ActsMatrix<3,3>&,MagneticFieldProvider::Cache&)
-  Vector3 getFieldGradient(const Vector3& position,
-                           ActsMatrix<3, 3>& derivative,
-                           MagneticFieldProvider::Cache& cache) const override {
+  Result<Vector3> getFieldGradient(
+      const Vector3& position, ActsMatrix<3, 3>& derivative,
+      MagneticFieldProvider::Cache& cache) const override {
     return m_bField->getFieldGradient(position, derivative, cache);
   }
 

Core/include/Acts/MagneticField/SolenoidBField.hpp

@@ -109,15 +109,16 @@ class SolenoidBField final : public MagneticFieldProvider {
 
   /// @copydoc MagneticFieldProvider::getField(const
   /// Vector3&,MagneticFieldProvider::Cache&)
-  Vector3 getField(const Vector3& position,
-                   MagneticFieldProvider::Cache& /*cache*/) const override;
+  Result<Vector3> getField(
+      const Vector3& position,
+      MagneticFieldProvider::Cache& /*cache*/) const override;
 
   /// @copydoc MagneticFieldProvider::getFieldGradient(const
   /// Vector3&,ActsMatrix<3,3>&,MagneticFieldProvider::Cache&)
   ///
   /// @note currently the derivative is not calculated
   /// @todo return derivative
-  Vector3 getFieldGradient(
+  Result<Vector3> getFieldGradient(
       const Vector3& position, ActsMatrix<3, 3>& /*derivative*/,
       MagneticFieldProvider::Cache& /*cache*/) const override;