Core/MagneticField/InterpolatedBFieldMapTests.cpp

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