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