Core/Utilities/BFieldMapUtilsTests.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2016-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/data/test_case.hpp>
0010 #include <boost/test/unit_test.hpp>
0011
0012 #include "Acts/Definitions/Algebra.hpp"
0013 #include "Acts/MagneticField/BFieldMapUtils.hpp"
0014 #include "Acts/MagneticField/InterpolatedBFieldMap.hpp"
0015 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0016 #include "Acts/Utilities/Result.hpp"
0017 #include "Acts/Utilities/VectorHelpers.hpp"
0018
0019 #include <array>
0020 #include <cstddef>
0021 #include <random>
0022 #include <utility>
0023 #include <vector>
0024
0025 namespace bdata = boost::unit_test::data;
0026
0027 using Acts::VectorHelpers::perp;
0028
0029 namespace Acts {
0030
0031 using namespace detail;
0032
0033 namespace Test {
0034
0035 BOOST_AUTO_TEST_CASE(bfield_creation) {
0036   // create grid values
0037   std::vector<double> rPos = {0., 1., 2., 3.};
0038   std::vector<double> xPos = {0., 1., 2., 3.};
0039   std::vector<double> yPos = {0., 1., 2., 3.};
0040   std::vector<double> zPos = {0., 1., 2., 3.};
0041
0042   // create b field in rz
0043   std::vector<Acts::Vector2> bField_rz;
0044   for (int i = 0; i < 27; i++) {
0045     bField_rz.push_back(Acts::Vector2(i, i));
0046   }
0047
0048   auto localToGlobalBin_rz = [](std::array<std::size_t, 2> binsRZ,
0049                                 std::array<std::size_t, 2> nBinsRZ) {
0050     return (binsRZ.at(1) * nBinsRZ.at(0) + binsRZ.at(0));
0051   };
0052   // create b field map in rz
0053   auto map_rz =
0054       Acts::fieldMapRZ(localToGlobalBin_rz, rPos, zPos, bField_rz, 1, 1, false);
0055   // check number of bins, minima & maxima
0056   std::vector<std::size_t> nBins_rz = {rPos.size(), zPos.size()};
0057   std::vector<double> minima_rz = {0., 0.};
0058   std::vector<double> maxima_rz = {3., 3.};
0059   BOOST_CHECK(map_rz.getNBins() == nBins_rz);
0060   // check minimum (should be first value because bin values are always
0061   // assigned to the left boundary)
0062   BOOST_CHECK(map_rz.getMin() == minima_rz);
0063   // check maximum (should be last value + 1 step because bin values are
0064   // always assigned to the left boundary)
0065   BOOST_CHECK(map_rz.getMax() == maxima_rz);
0066
0067   auto localToGlobalBin_xyz = [](std::array<std::size_t, 3> binsXYZ,
0068                                  std::array<std::size_t, 3> nBinsXYZ) {
0069     return (binsXYZ.at(0) * (nBinsXYZ.at(1) * nBinsXYZ.at(2)) +
0070             binsXYZ.at(1) * nBinsXYZ.at(2) + binsXYZ.at(2));
0071   };
0072
0073   rPos = {0., 1., 2., 3.};
0074   xPos = {0., 1., 2., 3.};
0075   yPos = {0., 1., 2., 3.};
0076   zPos = {0., 1., 2., 3.};
0077
0078   // create b field in xyz
0079   std::vector<Acts::Vector3> bField_xyz;
0080   for (int i = 0; i < 81; i++) {
0081     bField_xyz.push_back(Acts::Vector3(i, i, i));
0082   }
0083
0084   // create b field map in xyz
0085   auto map_xyz = Acts::fieldMapXYZ(localToGlobalBin_xyz, xPos, yPos, zPos,
0086                                    bField_xyz, 1, 1, false);
0087   // check number of bins, minima & maxima
0088   std::vector<std::size_t> nBins_xyz = {xPos.size(), yPos.size(), zPos.size()};
0089   std::vector<double> minima_xyz = {0., 0., 0.};
0090   std::vector<double> maxima_xyz = {3., 3., 3.};
0091   BOOST_CHECK(map_xyz.getNBins() == nBins_xyz);
0092   // check minimum (should be first value because bin values are always
0093   // assigned to the left boundary)
0094   BOOST_CHECK(map_xyz.getMin() == minima_xyz);
0095   // check maximum (should be last value + 1 step because bin values are
0096   // always assigned to the left boundary)
0097   BOOST_CHECK(map_xyz.getMax() == maxima_xyz);
0098
0099   // check if filled value is expected value in rz
0100   Acts::Vector3 pos0_rz(0., 0., 0.);
0101   Acts::Vector3 pos1_rz(1., 0., 1.);
0102   Acts::Vector3 pos2_rz(0., 2., 2.);
0103   auto value0_rz = map_rz.getField(pos0_rz).value();
0104   auto value1_rz = map_rz.getField(pos1_rz).value();
0105   auto value2_rz = map_rz.getField(pos2_rz).value();
0106   // calculate what the value should be at this point
0107   auto b0_rz =
0108       bField_rz.at(localToGlobalBin_rz({{0, 0}}, {{rPos.size(), zPos.size()}}));
0109   auto b1_rz =
0110       bField_rz.at(localToGlobalBin_rz({{1, 1}}, {{rPos.size(), zPos.size()}}));
0111   auto b2_rz =
0112       bField_rz.at(localToGlobalBin_rz({{2, 2}}, {{rPos.size(), zPos.size()}}));
0113   Acts::Vector3 bField0_rz(b0_rz.x(), 0., b0_rz.y());
0114   Acts::Vector3 bField1_rz(b1_rz.x(), 0., b1_rz.y());
0115   Acts::Vector3 bField2_rz(0., b2_rz.x(), b2_rz.y());
0116   // check the value
0117   // in rz case field is phi symmetric (check radius)
0118   CHECK_CLOSE_ABS(perp(value0_rz), perp(bField0_rz), 1e-9);
0119   CHECK_CLOSE_ABS(value0_rz.z(), bField0_rz.z(), 1e-9);
0120   CHECK_CLOSE_ABS(perp(value1_rz), perp(bField1_rz), 1e-9);
0121   CHECK_CLOSE_ABS(value1_rz.z(), bField1_rz.z(), 1e-9);
0122   CHECK_CLOSE_ABS(perp(value2_rz), perp(bField2_rz), 1e-9);
0123   CHECK_CLOSE_ABS(value2_rz.z(), bField2_rz.z(), 1e-9);
0124
0125   // check if filled value is expected value in rz
0126   Acts::Vector3 pos0_xyz(0., 0., 0.);
0127   Acts::Vector3 pos1_xyz(1., 1., 1.);
0128   Acts::Vector3 pos2_xyz(2., 2., 2.);
0129   auto value0_xyz = map_xyz.getField(pos0_xyz).value();
0130   auto value1_xyz = map_xyz.getField(pos1_xyz).value();
0131   auto value2_xyz = map_xyz.getField(pos2_xyz).value();
0132   // calculate what the value should be at this point
0133   auto b0_xyz = bField_xyz.at(localToGlobalBin_xyz(
0134       {{0, 0, 0}}, {{xPos.size(), yPos.size(), zPos.size()}}));
0135   auto b1_xyz = bField_xyz.at(localToGlobalBin_xyz(
0136       {{1, 1, 1}}, {{xPos.size(), yPos.size(), zPos.size()}}));
0137   auto b2_xyz = bField_xyz.at(localToGlobalBin_xyz(
0138       {{2, 2, 2}}, {{xPos.size(), yPos.size(), zPos.size()}}));
0139   // check the value
0140   BOOST_CHECK_EQUAL(value0_xyz, b0_xyz);
0141   BOOST_CHECK_EQUAL(value1_xyz, b1_xyz);
0142   BOOST_CHECK_EQUAL(value2_xyz, b2_xyz);
0143
0144   // checkx-,y-,z-symmetry - need to check close (because of interpolation
0145   // there can be small differences)
0146   CHECK_CLOSE_ABS(value0_xyz, b0_xyz, 1e-9);
0147   CHECK_CLOSE_ABS(value1_xyz, b1_xyz, 1e-9);
0148   CHECK_CLOSE_ABS(value2_xyz, b2_xyz, 1e-9);
0149 }
0150
0151 BOOST_AUTO_TEST_CASE(bfield_symmetry) {
0152   // create grid values
0153   std::vector<double> rPos = {0., 1., 2.};
0154   std::vector<double> xPos = {0., 1., 2.};
0155   std::vector<double> yPos = {0., 1., 2.};
0156   std::vector<double> zPos = {0., 1., 2.};
0157   // the bfield values in rz
0158   std::vector<Acts::Vector2> bField_rz;
0159   for (int i = 0; i < 9; i++) {
0160     bField_rz.push_back(Acts::Vector2(i, i));
0161   }
0162   // the field map in rz
0163   auto map_rz = Acts::fieldMapRZ(
0164       [](std::array<std::size_t, 2> binsRZ,
0165          std::array<std::size_t, 2> nBinsRZ) {
0166         return (binsRZ.at(1) * nBinsRZ.at(0) + binsRZ.at(0));
0167       },
0168       rPos, zPos, bField_rz, 1, 1, true);
0169
0170   // check number of bins, minima & maxima
0171   std::vector<std::size_t> nBins_rz = {rPos.size(), 2 * zPos.size() - 1};
0172   std::vector<double> minima_rz = {0., -2.};
0173   std::vector<double> maxima_rz = {2., 2.};
0174   BOOST_CHECK(map_rz.getNBins() == nBins_rz);
0175   auto vec = map_rz.getNBins();
0176   auto vec0 = map_rz.getMin();
0177   // check minimum (should be first value because bin values are always
0178   // assigned to the left boundary)
0179   BOOST_CHECK(map_rz.getMin() == minima_rz);
0180   // check maximum (should be last value + 1 step because bin values are
0181   // always assigned to the left boundary)
0182   BOOST_CHECK(map_rz.getMax() == maxima_rz);
0183
0184   // the bfield values in xyz
0185   std::vector<Acts::Vector3> bField_xyz;
0186   for (int i = 0; i < 27; i++) {
0187     bField_xyz.push_back(Acts::Vector3(i, i, i));
0188   }
0189   // the field map in xyz
0190   auto map_xyz = Acts::fieldMapXYZ(
0191       [](std::array<std::size_t, 3> binsXYZ,
0192          std::array<std::size_t, 3> nBinsXYZ) {
0193         return (binsXYZ.at(0) * (nBinsXYZ.at(1) * nBinsXYZ.at(2)) +
0194                 binsXYZ.at(1) * nBinsXYZ.at(2) + binsXYZ.at(2));
0195       },
0196       xPos, yPos, zPos, bField_xyz, 1, 1, true);
0197
0198   // check number of bins, minima & maxima
0199   std::vector<std::size_t> nBins_xyz = {
0200       2 * xPos.size() - 1, 2 * yPos.size() - 1, 2 * zPos.size() - 1};
0201   std::vector<double> minima_xyz = {-2., -2., -2.};
0202   std::vector<double> maxima_xyz = {2., 2., 2.};
0203   BOOST_CHECK(map_xyz.getNBins() == nBins_xyz);
0204   // check minimum (should be first value because bin values are always
0205   // assigned to the left boundary)
0206   BOOST_CHECK(map_xyz.getMin() == minima_xyz);
0207   // check maximum (should be last value + 1 step because bin values are
0208   // always assigned to the left boundary)
0209   BOOST_CHECK(map_xyz.getMax() == maxima_xyz);
0210
0211   Acts::Vector3 pos0(1.2, 1.3, 1.4);
0212   Acts::Vector3 pos1(1.2, 1.3, -1.4);
0213   Acts::Vector3 pos2(-1.2, 1.3, 1.4);
0214   Acts::Vector3 pos3(1.2, -1.3, 1.4);
0215   Acts::Vector3 pos4(-1.2, -1.3, 1.4);
0216
0217   auto value0_rz = map_rz.getField(pos0).value();
0218   auto value1_rz = map_rz.getField(pos1).value();
0219   auto value2_rz = map_rz.getField(pos2).value();
0220   auto value3_rz = map_rz.getField(pos3).value();
0221   auto value4_rz = map_rz.getField(pos4).value();
0222
0223   auto value0_xyz = map_xyz.getField(pos0).value();
0224   auto value1_xyz = map_xyz.getField(pos1).value();
0225   auto value2_xyz = map_xyz.getField(pos2).value();
0226   auto value3_xyz = map_xyz.getField(pos3).value();
0227   auto value4_xyz = map_xyz.getField(pos4).value();
0228
0229   // check z- and phi-symmetry
0230   CHECK_CLOSE_REL(perp(value0_rz), perp(value1_rz), 1e-10);
0231   CHECK_CLOSE_REL(value0_rz.z(), value1_rz.z(), 1e-10);
0232   CHECK_CLOSE_REL(perp(value0_rz), perp(value2_rz), 1e-10);
0233   CHECK_CLOSE_REL(value0_rz.z(), value2_rz.z(), 1e-10);
0234   CHECK_CLOSE_REL(perp(value0_rz), perp(value3_rz), 1e-10);
0235   CHECK_CLOSE_REL(value0_rz.z(), value3_rz.z(), 1e-10);
0236   CHECK_CLOSE_REL(perp(value0_rz), perp(value4_rz), 1e-10);
0237   CHECK_CLOSE_REL(value0_rz.z(), value4_rz.z(), 1e-10);
0238
0239   // checkx-,y-,z-symmetry - need to check close (because of interpolation
0240   // there can be small differences)
0241   CHECK_CLOSE_REL(value0_xyz, value1_xyz, 1e-10);
0242   CHECK_CLOSE_REL(value0_xyz, value2_xyz, 1e-10);
0243   CHECK_CLOSE_REL(value0_xyz, value3_xyz, 1e-10);
0244   CHECK_CLOSE_REL(value0_xyz, value4_xyz, 1e-10);
0245 }
0246
0247 /// Unit test for symmetric data
0248 BOOST_DATA_TEST_CASE(
0249     bfield_symmetry_random,
0250     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 0,
0251                    bdata::distribution =
0252                        std::uniform_real_distribution<double>(-10., 10.))) ^
0253         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 0,
0254                        bdata::distribution =
0255                            std::uniform_real_distribution<double>(-10., 10.))) ^
0256         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 0,
0257                        bdata::distribution =
0258                            std::uniform_real_distribution<double>(-20., 20.))) ^
0259         bdata::xrange(10),
0260     x, y, z, index) {
0261   (void)index;
0262
0263   std::vector<double> rPos;
0264   std::vector<double> xPos;
0265   std::vector<double> yPos;
0266   std::vector<double> zPos;
0267   double maxR = 20.;
0268   double maxZ = 30.;
0269   double maxBr = 10.;
0270   double maxBz = 20.;
0271   std::size_t nBins = 10;
0272   double stepR = maxR / nBins;
0273   double stepZ = maxZ / nBins;
0274   double bStepR = maxBr / nBins;
0275   double bStepZ = maxBz / nBins;
0276
0277   for (std::size_t i = 0; i < nBins; i++) {
0278     rPos.push_back(i * stepR);
0279     xPos.push_back(i * stepR);
0280     yPos.push_back(i * stepR);
0281     zPos.push_back(i * stepZ);
0282   }
0283   // bfield in rz
0284   std::vector<Acts::Vector2> bField_rz;
0285   for (std::size_t i = 0; i < nBins * nBins; i++) {
0286     bField_rz.push_back(Acts::Vector2(i * bStepR, i * bStepZ));
0287   }
0288   // the map in rz
0289   auto map_rz = Acts::fieldMapRZ(
0290       [](std::array<std::size_t, 2> binsRZ,
0291          std::array<std::size_t, 2> nBinsRZ) {
0292         return (binsRZ.at(1) * nBinsRZ.at(0) + binsRZ.at(0));
0293       },
0294       rPos, zPos, bField_rz, 1, 1, true);
0295
0296   // check number of bins, minima & maxima
0297   std::vector<std::size_t> nBins_rz = {rPos.size(), 2 * zPos.size() - 1};
0298   std::vector<double> minima_rz = {0., -((nBins - 1) * stepZ)};
0299   std::vector<double> maxima_rz = {(nBins - 1) * stepR, (nBins - 1) * stepZ};
0300   BOOST_CHECK(map_rz.getNBins() == nBins_rz);
0301   // check minimum (should be first value because bin values are always
0302   // assigned to the left boundary)
0303   CHECK_CLOSE_ABS(map_rz.getMin(), minima_rz, 1e-10);
0304   // check maximum (should be last value + 1 step because bin values are
0305   // always assigned to the left boundary)
0306   CHECK_CLOSE_ABS(map_rz.getMax(), maxima_rz, 1e-10);
0307
0308   // bfield in xyz
0309   std::vector<Acts::Vector3> bField_xyz;
0310   for (std::size_t i = 0; i < nBins * nBins * nBins; i++) {
0311     bField_xyz.push_back(Acts::Vector3(i * bStepR, i * bStepR, i * bStepZ));
0312   }
0313   // the map in xyz
0314   auto map_xyz = Acts::fieldMapXYZ(
0315       [](std::array<std::size_t, 3> binsXYZ,
0316          std::array<std::size_t, 3> nBinsXYZ) {
0317         return (binsXYZ.at(0) * (nBinsXYZ.at(1) * nBinsXYZ.at(2)) +
0318                 binsXYZ.at(1) * nBinsXYZ.at(2) + binsXYZ.at(2));
0319       },
0320       xPos, yPos, zPos, bField_xyz, 1, 1, true);
0321   // check number of bins, minima & maxima
0322   std::vector<std::size_t> nBins_xyz = {
0323       2 * xPos.size() - 1, 2 * yPos.size() - 1, 2 * zPos.size() - 1};
0324   std::vector<double> minima_xyz = {
0325       -((nBins - 1) * stepR), -((nBins - 1) * stepR), -((nBins - 1) * stepZ)};
0326   std::vector<double> maxima_xyz = {(nBins - 1) * stepR, (nBins - 1) * stepR,
0327                                     (nBins - 1) * stepZ};
0328   BOOST_CHECK(map_xyz.getNBins() == nBins_xyz);
0329   // check minimum (should be first value because bin values are always
0330   // assigned to the left boundary)
0331   CHECK_CLOSE_REL(map_xyz.getMin(), minima_xyz, 1e-10);
0332   // check maximum (should be last value + 1 step because bin values are
0333   // always assigned to the left boundary)
0334   CHECK_CLOSE_REL(map_xyz.getMax(), maxima_xyz, 1e-10);
0335
0336   Acts::Vector3 pos0(x, y, z);
0337   Acts::Vector3 pos1(x, y, -z);
0338   Acts::Vector3 pos2(-x, y, z);
0339   Acts::Vector3 pos3(x, -y, z);
0340   Acts::Vector3 pos4(-x, -y, z);
0341
0342   auto value0_rz = map_rz.getField(pos0).value();
0343   auto value1_rz = map_rz.getField(pos1).value();
0344   auto value2_rz = map_rz.getField(pos2).value();
0345   auto value3_rz = map_rz.getField(pos3).value();
0346   auto value4_rz = map_rz.getField(pos4).value();
0347
0348   // check z- and phi-symmetry
0349   CHECK_CLOSE_REL(perp(value0_rz), perp(value1_rz), 1e-10);
0350   CHECK_CLOSE_REL(value0_rz.z(), value1_rz.z(), 1e-10);
0351   CHECK_CLOSE_REL(perp(value0_rz), perp(value2_rz), 1e-10);
0352   CHECK_CLOSE_REL(value0_rz.z(), value2_rz.z(), 1e-10);
0353   CHECK_CLOSE_REL(perp(value0_rz), perp(value3_rz), 1e-10);
0354   CHECK_CLOSE_REL(value0_rz.z(), value3_rz.z(), 1e-10);
0355   CHECK_CLOSE_REL(perp(value0_rz), perp(value4_rz), 1e-10);
0356   CHECK_CLOSE_REL(value0_rz.z(), value4_rz.z(), 1e-10);
0357
0358   auto value0_xyz = map_xyz.getField(pos0).value();
0359   auto value1_xyz = map_xyz.getField(pos1).value();
0360   auto value2_xyz = map_xyz.getField(pos2).value();
0361   auto value3_xyz = map_xyz.getField(pos3).value();
0362   auto value4_xyz = map_xyz.getField(pos4).value();
0363
0364   // checkx-,y-,z-symmetry - need to check close (because of interpolation
0365   // there can be small differences)
0366   CHECK_CLOSE_REL(value0_xyz, value1_xyz, 1e-10);
0367   CHECK_CLOSE_REL(value0_xyz, value2_xyz, 1e-10);
0368   CHECK_CLOSE_REL(value0_xyz, value3_xyz, 1e-10);
0369   CHECK_CLOSE_REL(value0_xyz, value4_xyz, 1e-10);
0370 }
0371 }  // namespace Test
0372 }  // namespace Acts