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