Core/MagneticField/InterpolatedBFieldMapTests.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/MagneticField/InterpolatedBFieldMap.hpp"
0013 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0014 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0015 #include "Acts/Utilities/Axis.hpp"
0016 #include "Acts/Utilities/AxisDefinitions.hpp"
0017 #include "Acts/Utilities/Grid.hpp"
0018 #include "Acts/Utilities/Result.hpp"
0019 #include "Acts/Utilities/VectorHelpers.hpp"
0020
0021 #include <array>
0022 #include <cstddef>
0023 #include <optional>
0024 #include <utility>
0025
0026 using Acts::VectorHelpers::perp;
0027
0028 namespace Acts::Test {
0029
0030 // Create a test context
0031 MagneticFieldContext mfContext = MagneticFieldContext();
0032
0033 BOOST_AUTO_TEST_CASE(InterpolatedBFieldMap_rz) {
0034   // definition of dummy BField
0035   struct BField {
0036     static Vector3 value(const std::array<double, 2>& rz) {
0037       double r = rz.at(0);
0038       double z = rz.at(1);
0039       // linear in r and z so interpolation should be exact
0040       return Vector3(r * z, 3 * r, -2 * z);
0041     }
0042   };
0043
0044   // map (x,y,z) -> (r,z)
0045   auto transformPos = [](const Vector3& pos) {
0046     return Vector2(perp(pos), pos.z());
0047   };
0048
0049   // map (Bx,By,Bz) -> (Bx,By,Bz)
0050   auto transformBField = [](const Vector3& field, const Vector3&) {
0051     return field;
0052   };
0053
0054   // magnetic field known on grid in (r,z)
0055   Axis r(0.0, 4.0, 4u);
0056   Axis z(-5, 7, 6u);
0057
0058   Grid g(Type<Vector3>, std::move(r), std::move(z));
0059
0060   using Grid_t = decltype(g);
0061   using BField_t = InterpolatedBFieldMap<Grid_t>;
0062
0063   // set grid values
0064   for (std::size_t i = 1; i <= g.numLocalBins().at(0) + 1; ++i) {
0065     for (std::size_t j = 1; j <= g.numLocalBins().at(1) + 1; ++j) {
0066       Grid_t::index_t indices = {{i, j}};
0067       const auto& llCorner = g.lowerLeftBinEdge(indices);
0068       g.atLocalBins(indices) = BField::value(llCorner);
0069     }
0070   }
0071
0072   // create BField service
0073   BField_t b{{transformPos, transformBField, std::move(g)}};
0074
0075   auto bCacheAny = b.makeCache(mfContext);
0076   BField_t::Cache& bCache = bCacheAny.as<BField_t::Cache>();
0077
0078   auto check = [&](double i) {
0079     BOOST_CHECK(b.isInside({0, 0, i * 4.9}));
0080     BOOST_CHECK(!b.isInside({0, 0, i * 5.1}));
0081
0082     BOOST_CHECK(b.isInside({i * 2.9, 0, 0}));
0083     BOOST_CHECK(!b.isInside({i * 3.1, 0, 0}));
0084
0085     BOOST_CHECK(b.isInside({0, i * 2.9, 0}));
0086     BOOST_CHECK(!b.isInside({0, i * 3.1, 0}));
0087
0088     BOOST_CHECK(b.isInside({2, 2.2, 0}));
0089     BOOST_CHECK(!b.isInside({2, 3, 0}));
0090
0091     BOOST_CHECK(b.isInside({i * 2, 2.2, 0}));
0092     BOOST_CHECK(!b.isInside({i * 2, 3, 0}));
0093
0094     BOOST_CHECK(b.isInside({2, i * 2.2, 0}));
0095     BOOST_CHECK(!b.isInside({2, i * 3, 0}));
0096
0097     BOOST_CHECK(b.isInside({i * 2, i * 2.2, 0}));
0098     BOOST_CHECK(!b.isInside({i * 2, i * 3, 0}));
0099   };
0100
0101   check(1);
0102   check(-1);
0103
0104   Vector3 pos;
0105   pos << -3, 2.5,
0106       1.7;  // this was previously checked, but is actually out of bounds
0107   BOOST_CHECK(!b.isInside(pos));
0108   BOOST_CHECK(!b.getField(pos, bCacheAny).ok());
0109
0110   pos << -1.6, 2.5, 1.7;
0111   BOOST_CHECK(b.isInside(pos));
0112   CHECK_CLOSE_REL(b.getField(pos, bCacheAny).value(),
0113                   BField::value({{perp(pos), pos.z()}}), 1e-6);
0114
0115   auto& c = *bCache.fieldCell;
0116   BOOST_CHECK(c.isInside(transformPos(pos)));
0117   CHECK_CLOSE_REL(c.getField(transformPos(pos)),
0118                   BField::value({{perp(pos), pos.z()}}), 1e-6);
0119
0120   SquareMatrix3 deriv;
0121
0122   pos << 1, 1, -5.5;  // this position is outside the grid
0123   BOOST_CHECK(!b.isInside(pos));
0124   BOOST_CHECK(!b.getField(pos, bCacheAny).ok());
0125
0126   pos << 1, 6, -1.7;  // this position is outside the grid
0127   BOOST_CHECK(!b.isInside(pos));
0128   BOOST_CHECK(!b.getField(pos, bCacheAny).ok());
0129
0130   pos << 0, 1.5, -2.5;
0131   BOOST_CHECK(b.isInside(pos));
0132   bCacheAny = b.makeCache(mfContext);
0133   BField_t::Cache& bCache2 = bCacheAny.as<BField_t::Cache>();
0134   CHECK_CLOSE_REL(b.getField(pos, bCacheAny).value(),
0135                   BField::value({{perp(pos), pos.z()}}), 1e-6);
0136   c = *bCache2.fieldCell;
0137   BOOST_CHECK(c.isInside(transformPos(pos)));
0138   CHECK_CLOSE_REL(c.getField(transformPos(pos)),
0139                   BField::value({{perp(pos), pos.z()}}), 1e-6);
0140
0141   pos << 2, 3, -4;
0142   BOOST_CHECK(!b.isInside(pos));
0143
0144   pos << 2, 2.2, -4;
0145   BOOST_CHECK(b.isInside(pos));
0146   bCacheAny = b.makeCache(mfContext);
0147   BField_t::Cache& bCache3 = bCacheAny.as<BField_t::Cache>();
0148   CHECK_CLOSE_REL(b.getField(pos, bCacheAny).value(),
0149                   BField::value({{perp(pos), pos.z()}}), 1e-6);
0150   c = *bCache3.fieldCell;
0151   BOOST_CHECK(c.isInside(transformPos(pos)));
0152   CHECK_CLOSE_REL(c.getField(transformPos(pos)),
0153                   BField::value({{perp(pos), pos.z()}}), 1e-6);
0154
0155   // some field cell tests
0156   BOOST_CHECK(!c.isInside(transformPos((pos << 3, 2, -3.7).finished())));
0157   BOOST_CHECK(!c.isInside(transformPos((pos << -2, 3, -4.7).finished())));
0158   BOOST_CHECK(!c.isInside(transformPos((pos << -2, 3, 4.7).finished())));
0159   BOOST_CHECK(c.isInside(transformPos((pos << 0, 2, -4.7).finished())));
0160   BOOST_CHECK(!c.isInside(transformPos((pos << 5, 2, 14.).finished())));
0161 }
0162 }  // namespace Acts::Test