feat: Add a multi-region constant B-field provider (acts-project#3927)

In order to facilitate the ALICE3 use-case (as described to me by @Cas1997) we need to support magnetic fields that are constant in multiple regions, but which may have different field values in each of these regions. I implemented a new magnetic field provider based on a mapping of 3D ranges to field vectors which should support this use case.

The core idea is that the user provides a list of regions and the provider will check (with some simple caching) which region matches. The mechanism implemented is that -- in case of overlap -- the region later in the list takes precedence.

Also implements C++ and Python tests for this new field provider.

<!-- This is an auto-generated comment: release notes by coderabbit.ai -->
## Summary by CodeRabbit

## Release Notes

- **New Features**
	- Introduced the `MultiRangeBField` class for modeling magnetic fields with multiple overlapping regions.
	- Added the `RangeXDDim3` class for handling three-dimensional ranges in the geometry module.

- **Tests**
	- Implemented comprehensive tests for the `MultiRangeBField` class, covering various scenarios to validate its behavior.
	- Added a new test function for the `acts.MultiRangeBField` class in Python.

- **Documentation**
	- Updated documentation to include details about the `Acts::MultiRangeBField` class and its functionality, including caching mechanisms.

- **Chores**
	- Expanded the unit test suite to include tests for the `MultiRangeBField` class.
	- Added a new source file for the `MultiRangeBField` class to the project.
<!-- end of auto-generated comment: release notes by coderabbit.ai -->

Core/include/Acts/MagneticField/MultiRangeBField.hpp

@@ -0,0 +1,67 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 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 https://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/RangeXD.hpp"
+
+namespace Acts {
+
+/// @ingroup MagneticField
+///
+/// @brief Magnetic field provider modelling a magnetic field consisting of
+/// several (potentially overlapping) regions of constant values.
+class MultiRangeBField final : public MagneticFieldProvider {
+ private:
+  struct Cache {
+    explicit Cache(const MagneticFieldContext& /*unused*/);
+
+    std::optional<std::size_t> index = {};
+  };
+
+  using BFieldRange = std::pair<RangeXD<3, double>, Vector3>;
+
+  // The different ranges and their corresponding field vectors. Note that
+  // regions positioned _later_ in this vector take priority over earlier
+  // regions.
+  std::vector<BFieldRange> fieldRanges;
+
+ public:
+  /// @brief Construct a magnetic field from a vector of ranges.
+  ///
+  /// @warning These ranges are listed in increasing order of precedence,
+  /// i.e. ranges further along the vector have higher priority.
+  explicit MultiRangeBField(const std::vector<BFieldRange>& ranges);
+
+  explicit MultiRangeBField(std::vector<BFieldRange>&& ranges);
+
+  /// @brief Construct a cache object.
+  MagneticFieldProvider::Cache makeCache(
+      const MagneticFieldContext& mctx) const override;
+
+  /// @brief Request the value of the magnetic field at a given position.
+  ///
+  /// @param [in] position Global 3D position for the lookup.
+  /// @param [in, out] cache Cache object.
+  /// @returns A successful value containing a field vector if the given
+  /// location is contained inside any of the regions, or a failure value
+  /// otherwise.
+  Result<Vector3> getField(const Vector3& position,
+                           MagneticFieldProvider::Cache& cache) const override;
+
+  /// @brief Get the field gradient at a given position.
+  ///
+  /// @warning This is not currently implemented.
+  Result<Vector3> getFieldGradient(
+      const Vector3& position, ActsMatrix<3, 3>& /*unused*/,
+      MagneticFieldProvider::Cache& cache) const override;
+};
+}  // namespace Acts

Core/src/MagneticField/CMakeLists.txt

@@ -1,4 +1,8 @@
 target_sources(
     ActsCore
-    PRIVATE BFieldMapUtils.cpp SolenoidBField.cpp MagneticFieldError.cpp
+    PRIVATE
+        BFieldMapUtils.cpp
+        SolenoidBField.cpp
+        MagneticFieldError.cpp
+        MultiRangeBField.cpp
 )

Core/src/MagneticField/MultiRangeBField.cpp

@@ -0,0 +1,78 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
+
+#include "Acts/MagneticField/MultiRangeBField.hpp"
+
+namespace Acts {
+
+MultiRangeBField::Cache::Cache(const MagneticFieldContext& /*unused*/) {}
+
+MultiRangeBField::MultiRangeBField(const std::vector<BFieldRange>& ranges)
+    : fieldRanges(ranges) {}
+
+MultiRangeBField::MultiRangeBField(
+    std::vector<MultiRangeBField::BFieldRange>&& ranges)
+    : fieldRanges(std::move(ranges)) {}
+
+MagneticFieldProvider::Cache MultiRangeBField::makeCache(
+    const MagneticFieldContext& mctx) const {
+  return MagneticFieldProvider::Cache(std::in_place_type<Cache>, mctx);
+}
+
+Result<Vector3> MultiRangeBField::getField(
+    const Vector3& position, MagneticFieldProvider::Cache& cache) const {
+  // Because we assume that the regions are in increasing order of
+  // precedence, we can iterate over the array, remembering the _last_
+  // region that contained the given point. At the end of the loop, this
+  // region will then be the one we query for its field vector.
+  std::optional<std::size_t> foundRange = {};
+
+  // The cache object for this type contains an optional integer; if the
+  // integer is set, it indicates the index of the region in which the last
+  // access succeeded. This acts as a cache because if the current access is
+  // still within that region, it is guaranteed that none of the regions
+  // with lower priority -- i.e. that are earlier in the region vector --
+  // can be relevant to the current access. Thus, we request the cache index
+  // if it exists and perform a membership check on it; if that succeeds, we
+  // remember the corresponding region as a candidate.
+  if (Cache& lCache = cache.as<Cache>();
+      lCache.index.has_value() &&
+      std::get<0>(fieldRanges[*lCache.index])
+          .contains({position[0], position[1], position[2]})) {
+    foundRange = *lCache.index;
+  }
+
+  // Now, we iterate over the ranges. If we already have a range candidate,
+  // we start iteration just after the existing candidate; otherwise we
+  // start from the beginning.
+  for (std::size_t i = (foundRange.has_value() ? (*foundRange) + 1 : 0);
+       i < fieldRanges.size(); ++i) {
+    if (std::get<0>(fieldRanges[i])
+            .contains({position[0], position[1], position[2]})) {
+      foundRange = i;
+    }
+  }
+
+  // Update the cache using the result of this access.
+  cache.as<Cache>().index = foundRange;
+
+  // If we found a valid range, return the corresponding vector; otherwise
+  // return an out-of-bounds error.
+  if (foundRange.has_value()) {
+    return Result<Vector3>::success(std::get<1>(fieldRanges[*foundRange]));
+  } else {
+    return Result<Vector3>::failure(MagneticFieldError::OutOfBounds);
+  }
+}
+
+Result<Vector3> MultiRangeBField::getFieldGradient(
+    const Vector3& position, ActsMatrix<3, 3>& /*unused*/,
+    MagneticFieldProvider::Cache& cache) const {
+  return getField(position, cache);
+}
+}  // namespace Acts

Examples/Python/src/Geometry.cpp

@@ -492,6 +492,21 @@ void addExperimentalGeometry(Context& ctx) {
         .def(py::init<std::shared_ptr<KdtSurfacesDim2Bin100>, const Extent&>());
   }
 
+  {
+    using RangeXDDim3 = Acts::RangeXD<3u, double>;
+
+    py::class_<RangeXDDim3>(m, "RangeXDDim3")
+        .def(py::init([](const std::array<double, 2u>& range0,
+                         const std::array<double, 2u>& range1,
+                         const std::array<double, 2u>& range2) {
+          RangeXDDim3 range;
+          range[0].shrink(range0[0], range0[1]);
+          range[1].shrink(range1[0], range1[1]);
+          range[2].shrink(range2[0], range2[1]);
+          return range;
+        }));
+  }
+
   {
     // The external volume structure builder
     py::class_<Acts::Experimental::IExternalStructureBuilder,

Examples/Python/src/MagneticField.cpp

@@ -12,6 +12,7 @@
 #include "Acts/MagneticField/BFieldMapUtils.hpp"
 #include "Acts/MagneticField/ConstantBField.hpp"
 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
+#include "Acts/MagneticField/MultiRangeBField.hpp"
 #include "Acts/MagneticField/NullBField.hpp"
 #include "Acts/MagneticField/SolenoidBField.hpp"
 #include "Acts/Plugins/Python/Utilities.hpp"
@@ -86,6 +87,11 @@ void addMagneticField(Context& ctx) {
              std::shared_ptr<Acts::NullBField>>(m, "NullBField")
       .def(py::init<>());
 
+  py::class_<Acts::MultiRangeBField, Acts::MagneticFieldProvider,
+             std::shared_ptr<Acts::MultiRangeBField>>(m, "MultiRangeBField")
+      .def(py::init<
+           std::vector<std::pair<Acts::RangeXD<3, double>, Acts::Vector3>>>());
+
   {
     using Config = Acts::SolenoidBField::Config;
 

Examples/Python/tests/test_magnetic_field.py

@@ -72,3 +72,41 @@ def test_solenoid(conf_const):
     )
 
     assert isinstance(field, acts.examples.InterpolatedMagneticField2)
+
+
+def test_multiregion_bfield():
+    with pytest.raises(TypeError):
+        acts.MultiRangeBField()
+
+    rs = [
+        (acts.RangeXDDim3((0, 3), (0, 3), (0, 3)), acts.Vector3(0.0, 0.0, 2.0)),
+        (acts.RangeXDDim3((1, 2), (1, 2), (1, 10)), acts.Vector3(2.0, 0.0, 0.0)),
+    ]
+    f = acts.MultiRangeBField(rs)
+    assert f
+
+    ctx = acts.MagneticFieldContext()
+    assert ctx
+
+    fc = f.makeCache(ctx)
+    assert fc
+
+    rv = f.getField(acts.Vector3(0.5, 0.5, 0.5), fc)
+    assert rv[0] == pytest.approx(0.0)
+    assert rv[1] == pytest.approx(0.0)
+    assert rv[2] == pytest.approx(2.0)
+
+    rv = f.getField(acts.Vector3(1.5, 1.5, 5.0), fc)
+    assert rv[0] == pytest.approx(2.0)
+    assert rv[1] == pytest.approx(0.0)
+    assert rv[2] == pytest.approx(0.0)
+
+    rv = f.getField(acts.Vector3(2.5, 2.5, 2.5), fc)
+    assert rv[0] == pytest.approx(0.0)
+    assert rv[1] == pytest.approx(0.0)
+    assert rv[2] == pytest.approx(2.0)
+
+    rv = f.getField(acts.Vector3(1.5, 1.5, 1.5), fc)
+    assert rv[0] == pytest.approx(2.0)
+    assert rv[1] == pytest.approx(0.0)
+    assert rv[2] == pytest.approx(0.0)

Tests/UnitTests/Core/MagneticField/CMakeLists.txt

@@ -1,4 +1,5 @@
 add_unittest(ConstantBField ConstantBFieldTests.cpp)
 add_unittest(InterpolatedBFieldMap InterpolatedBFieldMapTests.cpp)
 add_unittest(SolenoidBField SolenoidBFieldTests.cpp)
+add_unittest(MultiRangeBField MultiRangeBFieldTests.cpp)
 add_unittest(MagneticFieldProvider MagneticFieldProviderTests.cpp)

Tests/UnitTests/Core/MagneticField/MultiRangeBFieldTests.cpp

@@ -0,0 +1,74 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
+
+#include <boost/test/floating_point_comparison.hpp>
+#include <boost/test/unit_test.hpp>
+
+#include "Acts/MagneticField/MagneticFieldContext.hpp"
+#include "Acts/MagneticField/MultiRangeBField.hpp"
+#include "Acts/Utilities/Result.hpp"
+
+namespace Acts::Test {
+
+MagneticFieldContext mfContext = MagneticFieldContext();
+
+BOOST_AUTO_TEST_CASE(TestMultiRangeBField) {
+  std::vector<std::pair<RangeXD<3, double>, Vector3>> inputs;
+
+  inputs.emplace_back(RangeXD<3, double>{{0., 0., 0.}, {3., 3., 3.}},
+                      Vector3{0., 0., 2.});
+  inputs.emplace_back(RangeXD<3, double>{{1., 1., 1.}, {2., 2., 10.}},
+                      Vector3{2., 0., 0.});
+
+  const MultiRangeBField bfield(std::move(inputs));
+
+  auto bcache = bfield.makeCache(mfContext);
+
+  // Test a point inside the first volume.
+  {
+    Result<Vector3> r = bfield.getField({0.5, 0.5, 0.5}, bcache);
+    BOOST_CHECK(r.ok());
+    BOOST_CHECK_CLOSE((*r)[0], 0., 0.05);
+    BOOST_CHECK_CLOSE((*r)[1], 0., 0.05);
+    BOOST_CHECK_CLOSE((*r)[2], 2., 0.05);
+  }
+
+  // Test a point inside the second volume, not overlapping the first.
+  {
+    Result<Vector3> r = bfield.getField({1.5, 1.5, 5.}, bcache);
+    BOOST_CHECK(r.ok());
+    BOOST_CHECK_CLOSE((*r)[0], 2., 0.05);
+    BOOST_CHECK_CLOSE((*r)[1], 0., 0.05);
+    BOOST_CHECK_CLOSE((*r)[2], 0., 0.05);
+  }
+
+  // Test a point inside the first volume again.
+  {
+    Result<Vector3> r = bfield.getField({2.5, 2.5, 2.5}, bcache);
+    BOOST_CHECK(r.ok());
+    BOOST_CHECK_CLOSE((*r)[0], 0., 0.05);
+    BOOST_CHECK_CLOSE((*r)[1], 0., 0.05);
+    BOOST_CHECK_CLOSE((*r)[2], 2., 0.05);
+  }
+
+  // Test a point inside the second volume in the overlap region.
+  {
+    Result<Vector3> r = bfield.getField({1.5, 1.5, 1.5}, bcache);
+    BOOST_CHECK(r.ok());
+    BOOST_CHECK_CLOSE((*r)[0], 2., 0.05);
+    BOOST_CHECK_CLOSE((*r)[1], 0., 0.05);
+    BOOST_CHECK_CLOSE((*r)[2], 0., 0.05);
+  }
+
+  // Test a point outside all volumes.
+  {
+    Result<Vector3> r = bfield.getField({-1., -1., -1}, bcache);
+    BOOST_CHECK(!r.ok());
+  }
+}
+}  // namespace Acts::Test

docs/core/magnetic_field.md

@@ -235,6 +235,28 @@ analytical implementation and is much faster to lookup:
 :::{doxygenfunction} Acts::solenoidFieldMap
 :::
 
+### Multi-range constant field
+
+The multi-range constant field allows modelling cases where a magnetic field
+can be described as multiple (potentially overlapping) regions, each of which
+has its own constant magnetic field. This provides more flexibility than the
+{class}`Acts::ConstantBField` while providing higher performance than
+{class}`Acts::InterpolatedBFieldMap`.
+
+This magnetic field provider is configured using a list of pairs, where each
+pair defines a region in three-dimensional space as well as a field vector.
+Magnetic field lookup then proceeds by finding the _last_ region in the
+user-provided list that contains the requested coordinate and returning the
+corresponding field vector.
+
+The implementation uses a simple caching mechanism to store the last matched
+region, providing improved performance for consecutive lookups within the same
+region. This is thread-safe when each thread uses its own cache instance. The
+field configuration itself is immutable after construction.
+
+:::{doxygenclass} Acts::MultiRangeBField
+:::
+
 ## Full provider interface
 
 :::{doxygenclass} Acts::MagneticFieldProvider