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   SquareMatrix3 deriv;
0125
0126   pos << 1, 1, -5.5;  // this position is outside the grid
0127   BOOST_CHECK(!b.isInside(pos));
0128   BOOST_CHECK(!b.getField(pos, bCacheAny).ok());
0129
0130   pos << 1, 6, -1.7;  // this position is outside the grid
0131   BOOST_CHECK(!b.isInside(pos));
0132   BOOST_CHECK(!b.getField(pos, bCacheAny).ok());
0133
0134   pos << 0, 1.5, -2.5;
0135   BOOST_CHECK(b.isInside(pos));
0136   bCacheAny = b.makeCache(mfContext);
0137   BField_t::Cache& bCache2 = bCacheAny.as<BField_t::Cache>();
0138   CHECK_CLOSE_REL(b.getField(pos, bCacheAny).value(),
0139                   BField::value({{perp(pos), pos.z()}}), 1e-6);
0140   c = *bCache2.fieldCell;
0141   BOOST_CHECK(c.isInside(transformPos(pos)));
0142   CHECK_CLOSE_REL(c.getField(transformPos(pos)),
0143                   BField::value({{perp(pos), pos.z()}}), 1e-6);
0144
0145   pos << 2, 3, -4;
0146   BOOST_CHECK(!b.isInside(pos));
0147
0148   pos << 2, 2.2, -4;
0149   BOOST_CHECK(b.isInside(pos));
0150   bCacheAny = b.makeCache(mfContext);
0151   BField_t::Cache& bCache3 = bCacheAny.as<BField_t::Cache>();
0152   CHECK_CLOSE_REL(b.getField(pos, bCacheAny).value(),
0153                   BField::value({{perp(pos), pos.z()}}), 1e-6);
0154   c = *bCache3.fieldCell;
0155   BOOST_CHECK(c.isInside(transformPos(pos)));
0156   CHECK_CLOSE_REL(c.getField(transformPos(pos)),
0157                   BField::value({{perp(pos), pos.z()}}), 1e-6);
0158
0159   // some field cell tests
0160   BOOST_CHECK(!c.isInside(transformPos((pos << 3, 2, -3.7).finished())));
0161   BOOST_CHECK(!c.isInside(transformPos((pos << -2, 3, -4.7).finished())));
0162   BOOST_CHECK(!c.isInside(transformPos((pos << -2, 3, 4.7).finished())));
0163   BOOST_CHECK(c.isInside(transformPos((pos << 0, 2, -4.7).finished())));
0164   BOOST_CHECK(!c.isInside(transformPos((pos << 5, 2, 14.).finished())));
0165 }
0166
0167 BOOST_AUTO_TEST_SUITE_END()
0168
0169 }  // namespace ActsTests