Tests/IntegrationTests/InterpolatedSolenoidBFieldTest.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/data/test_case.hpp>
0010 #include <boost/test/unit_test.hpp>
0011
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/MagneticField/BFieldMapUtils.hpp"
0014 #include "Acts/MagneticField/InterpolatedBFieldMap.hpp"
0015 #include "Acts/MagneticField/SolenoidBField.hpp"
0016 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0017 #include "Acts/Utilities/Grid.hpp"
0018 #include "Acts/Utilities/detail/Axis.hpp"
0019
0020 #include <cmath>
0021 #include <fstream>
0022 #include <iostream>
0023 #include <random>
0024
0025 using namespace Acts::UnitLiterals;
0026
0027 namespace bdata = boost::unit_test::data;
0028
0029 namespace Acts::IntegrationTest {
0030
0031 const double L = 5.8_m;
0032 const double R = (2.56 + 2.46) * 0.5 * 0.5_m;
0033 const std::size_t nCoils = 1154;
0034 const double bMagCenter = 2_T;
0035 const std::size_t nBinsR = 150;
0036 const std::size_t nBinsZ = 200;
0037
0038 auto makeFieldMap(const SolenoidBField& field) {
0039   std::ofstream ostr("solenoidmap.csv");
0040   ostr << "i;j;r;z;B_r;B_z" << std::endl;
0041
0042   double rMin = 0;
0043   double rMax = R * 2.;
0044   double zMin = 2 * (-L / 2.);
0045   double zMax = 2 * (L / 2.);
0046
0047   std::cout << "rMin = " << rMin << std::endl;
0048   std::cout << "rMax = " << rMax << std::endl;
0049   std::cout << "zMin = " << zMin << std::endl;
0050   std::cout << "zMax = " << zMax << std::endl;
0051
0052   auto map =
0053       solenoidFieldMap({rMin, rMax}, {zMin, zMax}, {nBinsR, nBinsZ}, field);
0054   // I know this is the correct grid type
0055   using Grid_t = Acts::Grid<Acts::Vector2, Acts::detail::EquidistantAxis,
0056                             Acts::detail::EquidistantAxis>;
0057   const Grid_t& grid = map.getGrid();
0058   using index_t = Grid_t::index_t;
0059   using point_t = Grid_t::point_t;
0060
0061   for (std::size_t i = 0; i <= nBinsR + 1; i++) {
0062     for (std::size_t j = 0; j <= nBinsZ + 1; j++) {
0063       // std::cout << "(i,j) = " << i << "," << j << std::endl;
0064       index_t index({i, j});
0065       if (i == 0 || j == 0 || i == nBinsR + 1 || j == nBinsZ + 1) {
0066         // under or overflow bin
0067       } else {
0068         point_t lowerLeft = grid.lowerLeftBinEdge(index);
0069         Vector2 B = grid.atLocalBins(index);
0070         ostr << i << ";" << j << ";" << lowerLeft[0] << ";" << lowerLeft[1];
0071         ostr << ";" << B[0] << ";" << B[1] << std::endl;
0072       }
0073     }
0074   }
0075
0076   return map;
0077 }
0078
0079 Acts::SolenoidBField bSolenoidField({R, L, nCoils, bMagCenter});
0080 auto bFieldMap = makeFieldMap(bSolenoidField);
0081 auto bCache = bFieldMap.makeCache(Acts::MagneticFieldContext{});
0082
0083 struct StreamWrapper {
0084   StreamWrapper(std::ofstream ofstr) : m_ofstr(std::move(ofstr)) {
0085     m_ofstr << "x;y;z;B_x;B_y;B_z;Bm_x;Bm_y;Bm_z" << std::endl;
0086   }
0087
0088   std::ofstream m_ofstr;
0089 };
0090
0091 StreamWrapper valid(std::ofstream("magfield_lookup.csv"));
0092
0093 const int ntests = 10000;
0094 BOOST_DATA_TEST_CASE(
0095     solenoid_interpolated_bfield_comparison,
0096     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 1,
0097                    bdata::distribution = std::uniform_real_distribution<double>(
0098                        1.5 * (-L / 2.), 1.5 * L / 2.))) ^
0099         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 2,
0100                        bdata::distribution =
0101                            std::uniform_real_distribution<double>(0,
0102                                                                   R * 1.5))) ^
0103         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 3,
0104                        bdata::distribution =
0105                            std::uniform_real_distribution<double>(-M_PI,
0106                                                                   M_PI))) ^
0107         bdata::xrange(ntests),
0108     z, r, phi, index) {
0109   if (index % 1000 == 0) {
0110     std::cout << index << std::endl;
0111   }
0112
0113   Vector3 pos(r * std::cos(phi), r * std::sin(phi), z);
0114   Vector3 B = bSolenoidField.getField(pos) / Acts::UnitConstants::T;
0115   Vector3 Bm = bFieldMap.getField(pos, bCache).value() / Acts::UnitConstants::T;
0116
0117   // test less than 5% deviation
0118   if (std::abs(r - R) > 10 && (std::abs(z) < L / 3. || r > 20)) {
0119     // only if more than 10mm away from coil for all z
0120     // only if not close to r=0 for large z
0121     CHECK_CLOSE_REL(Bm.norm(), B.norm(), 0.05);
0122   }
0123
0124   std::ofstream& ofstr = valid.m_ofstr;
0125   ofstr << pos.x() << ";" << pos.y() << ";" << pos.z() << ";";
0126   ofstr << B.x() << ";" << B.y() << ";" << B.z() << ";";
0127   ofstr << Bm.x() << ";" << Bm.y() << ";" << Bm.z() << std::endl;
0128 }
0129
0130 }  // namespace Acts::IntegrationTest