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