From b3933b06cd9722f52f568cae3513123f343923e4 Mon Sep 17 00:00:00 2001 From: Stephen Nicholas Swatman Date: Thu, 5 Dec 2024 14:26:33 +0100 Subject: [PATCH] feat: Add a multi-region constant B-field provider (#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. ## 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. --- .../Acts/MagneticField/MultiRangeBField.hpp | 67 ++++++++++++++++ Core/src/MagneticField/CMakeLists.txt | 6 +- Core/src/MagneticField/MultiRangeBField.cpp | 78 +++++++++++++++++++ Examples/Python/src/Geometry.cpp | 15 ++++ Examples/Python/src/MagneticField.cpp | 6 ++ Examples/Python/tests/test_magnetic_field.py | 38 +++++++++ .../Core/MagneticField/CMakeLists.txt | 1 + .../MagneticField/MultiRangeBFieldTests.cpp | 74 ++++++++++++++++++ docs/core/magnetic_field.md | 22 ++++++ 9 files changed, 306 insertions(+), 1 deletion(-) create mode 100644 Core/include/Acts/MagneticField/MultiRangeBField.hpp create mode 100644 Core/src/MagneticField/MultiRangeBField.cpp create mode 100644 Tests/UnitTests/Core/MagneticField/MultiRangeBFieldTests.cpp diff --git a/Core/include/Acts/MagneticField/MultiRangeBField.hpp b/Core/include/Acts/MagneticField/MultiRangeBField.hpp new file mode 100644 index 00000000000..60b90097161 --- /dev/null +++ b/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 index = {}; + }; + + using BFieldRange = std::pair, Vector3>; + + // The different ranges and their corresponding field vectors. Note that + // regions positioned _later_ in this vector take priority over earlier + // regions. + std::vector 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& ranges); + + explicit MultiRangeBField(std::vector&& 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 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 getFieldGradient( + const Vector3& position, ActsMatrix<3, 3>& /*unused*/, + MagneticFieldProvider::Cache& cache) const override; +}; +} // namespace Acts diff --git a/Core/src/MagneticField/CMakeLists.txt b/Core/src/MagneticField/CMakeLists.txt index 4b9ebb0be0f..3b164d39505 100644 --- a/Core/src/MagneticField/CMakeLists.txt +++ b/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 ) diff --git a/Core/src/MagneticField/MultiRangeBField.cpp b/Core/src/MagneticField/MultiRangeBField.cpp new file mode 100644 index 00000000000..8899b50d802 --- /dev/null +++ b/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& ranges) + : fieldRanges(ranges) {} + +MultiRangeBField::MultiRangeBField( + std::vector&& ranges) + : fieldRanges(std::move(ranges)) {} + +MagneticFieldProvider::Cache MultiRangeBField::makeCache( + const MagneticFieldContext& mctx) const { + return MagneticFieldProvider::Cache(std::in_place_type, mctx); +} + +Result 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 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(); + 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().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::success(std::get<1>(fieldRanges[*foundRange])); + } else { + return Result::failure(MagneticFieldError::OutOfBounds); + } +} + +Result MultiRangeBField::getFieldGradient( + const Vector3& position, ActsMatrix<3, 3>& /*unused*/, + MagneticFieldProvider::Cache& cache) const { + return getField(position, cache); +} +} // namespace Acts diff --git a/Examples/Python/src/Geometry.cpp b/Examples/Python/src/Geometry.cpp index d69b35a86ba..4c15e0dcf94 100644 --- a/Examples/Python/src/Geometry.cpp +++ b/Examples/Python/src/Geometry.cpp @@ -492,6 +492,21 @@ void addExperimentalGeometry(Context& ctx) { .def(py::init, const Extent&>()); } + { + using RangeXDDim3 = Acts::RangeXD<3u, double>; + + py::class_(m, "RangeXDDim3") + .def(py::init([](const std::array& range0, + const std::array& range1, + const std::array& 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_>(m, "NullBField") .def(py::init<>()); + py::class_>(m, "MultiRangeBField") + .def(py::init< + std::vector, Acts::Vector3>>>()); + { using Config = Acts::SolenoidBField::Config; diff --git a/Examples/Python/tests/test_magnetic_field.py b/Examples/Python/tests/test_magnetic_field.py index a51ebd55ab6..d1336e3ca8a 100644 --- a/Examples/Python/tests/test_magnetic_field.py +++ b/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) diff --git a/Tests/UnitTests/Core/MagneticField/CMakeLists.txt b/Tests/UnitTests/Core/MagneticField/CMakeLists.txt index 60a829bf548..4600cc8352e 100644 --- a/Tests/UnitTests/Core/MagneticField/CMakeLists.txt +++ b/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) diff --git a/Tests/UnitTests/Core/MagneticField/MultiRangeBFieldTests.cpp b/Tests/UnitTests/Core/MagneticField/MultiRangeBFieldTests.cpp new file mode 100644 index 00000000000..3644db058b4 --- /dev/null +++ b/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 +#include + +#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, 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 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 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 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 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 r = bfield.getField({-1., -1., -1}, bcache); + BOOST_CHECK(!r.ok()); + } +} +} // namespace Acts::Test diff --git a/docs/core/magnetic_field.md b/docs/core/magnetic_field.md index 1924234bde0..488d5890203 100644 --- a/docs/core/magnetic_field.md +++ b/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