EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BFieldMapUtilsTests.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BFieldMapUtilsTests.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2016-2018 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include <boost/test/data/test_case.hpp>
10 #include <boost/test/unit_test.hpp>
11 
12 #include "Acts/MagneticField/BFieldMapUtils.hpp"
13 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
14 #include "Acts/Utilities/Helpers.hpp"
15 #include "Acts/Utilities/detail/Axis.hpp"
16 
17 namespace bdata = boost::unit_test::data;
18 
19 using Acts::VectorHelpers::perp;
20 
21 namespace Acts {
22 
23 using namespace detail;
24 
25 namespace Test {
26 
27 BOOST_AUTO_TEST_CASE(bfield_creation) {
28  // create grid values
29  std::vector<double> rPos = {0., 1., 2.};
30  std::vector<double> xPos = {0., 1., 2.};
31  std::vector<double> yPos = {0., 1., 2.};
32  std::vector<double> zPos = {0., 1., 2.};
33 
34  // create b field in rz
35  std::vector<Acts::Vector2D> bField_rz;
36  for (int i = 0; i < 9; i++) {
37  bField_rz.push_back(Acts::Vector2D(i, i));
38  }
39 
40  auto localToGlobalBin_rz = [](std::array<size_t, 2> binsRZ,
41  std::array<size_t, 2> nBinsRZ) {
42  return (binsRZ.at(1) * nBinsRZ.at(0) + binsRZ.at(0));
43  };
44  // create b field mapper in rz
45  auto mapper_rz = Acts::fieldMapperRZ(localToGlobalBin_rz, rPos, zPos,
46  bField_rz, 1, 1, false);
47  // check number of bins, minima & maxima
48  std::vector<size_t> nBins_rz = {rPos.size(), zPos.size()};
49  std::vector<double> minima_rz = {0., 0.};
50  std::vector<double> maxima_rz = {3., 3.};
51  BOOST_CHECK(mapper_rz.getNBins() == nBins_rz);
52  // check minimum (should be first value because bin values are always
53  // assigned to the left boundary)
54  BOOST_CHECK(mapper_rz.getMin() == minima_rz);
55  // check maximum (should be last value + 1 step because bin values are
56  // always assigned to the left boundary)
57  BOOST_CHECK(mapper_rz.getMax() == maxima_rz);
58 
59  // create b field in xyz
60  std::vector<Acts::Vector3D> bField_xyz;
61  for (int i = 0; i < 27; i++) {
62  bField_xyz.push_back(Acts::Vector3D(i, i, i));
63  }
64 
65  auto localToGlobalBin_xyz = [](std::array<size_t, 3> binsXYZ,
66  std::array<size_t, 3> nBinsXYZ) {
67  return (binsXYZ.at(0) * (nBinsXYZ.at(1) * nBinsXYZ.at(2)) +
68  binsXYZ.at(1) * nBinsXYZ.at(2) + binsXYZ.at(2));
69  };
70 
71  // create b field mapper in xyz
72  auto mapper_xyz = Acts::fieldMapperXYZ(localToGlobalBin_xyz, xPos, yPos, zPos,
73  bField_xyz, 1, 1, false);
74  // check number of bins, minima & maxima
75  std::vector<size_t> nBins_xyz = {xPos.size(), yPos.size(), zPos.size()};
76  std::vector<double> minima_xyz = {0., 0., 0.};
77  std::vector<double> maxima_xyz = {3., 3., 3.};
78  BOOST_CHECK(mapper_xyz.getNBins() == nBins_xyz);
79  // check minimum (should be first value because bin values are always
80  // assigned to the left boundary)
81  BOOST_CHECK(mapper_xyz.getMin() == minima_xyz);
82  // check maximum (should be last value + 1 step because bin values are
83  // always assigned to the left boundary)
84  BOOST_CHECK(mapper_xyz.getMax() == maxima_xyz);
85 
86  // check if filled value is expected value in rz
87  Acts::Vector3D pos0_rz(0., 0., 0.);
88  Acts::Vector3D pos1_rz(1., 0., 1.);
89  Acts::Vector3D pos2_rz(0., 2., 2.);
90  auto value0_rz = mapper_rz.getField(pos0_rz);
91  auto value1_rz = mapper_rz.getField(pos1_rz);
92  auto value2_rz = mapper_rz.getField(pos2_rz);
93  // calculate what the value should be at this point
94  auto b0_rz =
95  bField_rz.at(localToGlobalBin_rz({{0, 0}}, {{rPos.size(), zPos.size()}}));
96  auto b1_rz =
97  bField_rz.at(localToGlobalBin_rz({{1, 1}}, {{rPos.size(), zPos.size()}}));
98  auto b2_rz =
99  bField_rz.at(localToGlobalBin_rz({{2, 2}}, {{rPos.size(), zPos.size()}}));
100  Acts::Vector3D bField0_rz(b0_rz.x(), 0., b0_rz.y());
101  Acts::Vector3D bField1_rz(b1_rz.x(), 0., b1_rz.y());
102  Acts::Vector3D bField2_rz(0., b2_rz.x(), b2_rz.y());
103  // check the value
104  // in rz case field is phi symmetric (check radius)
105  CHECK_CLOSE_ABS(perp(value0_rz), perp(bField0_rz), 1e-9);
106  CHECK_CLOSE_ABS(value0_rz.z(), bField0_rz.z(), 1e-9);
107  CHECK_CLOSE_ABS(perp(value1_rz), perp(bField1_rz), 1e-9);
108  CHECK_CLOSE_ABS(value1_rz.z(), bField1_rz.z(), 1e-9);
109  CHECK_CLOSE_ABS(perp(value2_rz), perp(bField2_rz), 1e-9);
110  CHECK_CLOSE_ABS(value2_rz.z(), bField2_rz.z(), 1e-9);
111 
112  // check if filled value is expected value in rz
113  Acts::Vector3D pos0_xyz(0., 0., 0.);
114  Acts::Vector3D pos1_xyz(1., 1., 1.);
115  Acts::Vector3D pos2_xyz(2., 2., 2.);
116  auto value0_xyz = mapper_xyz.getField(pos0_xyz);
117  auto value1_xyz = mapper_xyz.getField(pos1_xyz);
118  auto value2_xyz = mapper_xyz.getField(pos2_xyz);
119  // calculate what the value should be at this point
120  auto b0_xyz = bField_xyz.at(localToGlobalBin_xyz(
121  {{0, 0, 0}}, {{xPos.size(), yPos.size(), zPos.size()}}));
122  auto b1_xyz = bField_xyz.at(localToGlobalBin_xyz(
123  {{1, 1, 1}}, {{xPos.size(), yPos.size(), zPos.size()}}));
124  auto b2_xyz = bField_xyz.at(localToGlobalBin_xyz(
125  {{2, 2, 2}}, {{xPos.size(), yPos.size(), zPos.size()}}));
126  // check the value
127  BOOST_CHECK_EQUAL(value0_xyz, b0_xyz);
128  BOOST_CHECK_EQUAL(value1_xyz, b1_xyz);
129  BOOST_CHECK_EQUAL(value2_xyz, b2_xyz);
130 
131  // checkx-,y-,z-symmetry - need to check close (because of interpolation
132  // there can be small differences)
133  CHECK_CLOSE_ABS(value0_xyz, b0_xyz, 1e-9);
134  CHECK_CLOSE_ABS(value1_xyz, b1_xyz, 1e-9);
135  CHECK_CLOSE_ABS(value2_xyz, b2_xyz, 1e-9);
136 }
137 
138 BOOST_AUTO_TEST_CASE(bfield_symmetry) {
139  // create grid values
140  std::vector<double> rPos = {0., 1., 2.};
141  std::vector<double> xPos = {0., 1., 2.};
142  std::vector<double> yPos = {0., 1., 2.};
143  std::vector<double> zPos = {0., 1., 2.};
144  // the bfield values in rz
145  std::vector<Acts::Vector2D> bField_rz;
146  for (int i = 0; i < 9; i++) {
147  bField_rz.push_back(Acts::Vector2D(i, i));
148  }
149  // the field mapper in rz
150  auto mapper_rz = Acts::fieldMapperRZ(
151  [](std::array<size_t, 2> binsRZ, std::array<size_t, 2> nBinsRZ) {
152  return (binsRZ.at(1) * nBinsRZ.at(0) + binsRZ.at(0));
153  },
154  rPos, zPos, bField_rz, 1, 1, true);
155 
156  // check number of bins, minima & maxima
157  std::vector<size_t> nBins_rz = {rPos.size(), 2 * zPos.size() - 1};
158  std::vector<double> minima_rz = {0., -2.};
159  std::vector<double> maxima_rz = {3., 3.};
160  BOOST_CHECK(mapper_rz.getNBins() == nBins_rz);
161  auto vec = mapper_rz.getNBins();
162  auto vec0 = mapper_rz.getMin();
163  // check minimum (should be first value because bin values are always
164  // assigned to the left boundary)
165  BOOST_CHECK(mapper_rz.getMin() == minima_rz);
166  // check maximum (should be last value + 1 step because bin values are
167  // always assigned to the left boundary)
168  BOOST_CHECK(mapper_rz.getMax() == maxima_rz);
169 
170  // the bfield values in xyz
171  std::vector<Acts::Vector3D> bField_xyz;
172  for (int i = 0; i < 27; i++) {
173  bField_xyz.push_back(Acts::Vector3D(i, i, i));
174  }
175  // the field mapper in xyz
176  auto mapper_xyz = Acts::fieldMapperXYZ(
177  [](std::array<size_t, 3> binsXYZ, std::array<size_t, 3> nBinsXYZ) {
178  return (binsXYZ.at(0) * (nBinsXYZ.at(1) * nBinsXYZ.at(2)) +
179  binsXYZ.at(1) * nBinsXYZ.at(2) + binsXYZ.at(2));
180  },
181  xPos, yPos, zPos, bField_xyz, 1, 1, true);
182 
183  // check number of bins, minima & maxima
184  std::vector<size_t> nBins_xyz = {2 * xPos.size() - 1, 2 * yPos.size() - 1,
185  2 * zPos.size() - 1};
186  std::vector<double> minima_xyz = {-2., -2., -2.};
187  std::vector<double> maxima_xyz = {3., 3., 3.};
188  BOOST_CHECK(mapper_xyz.getNBins() == nBins_xyz);
189  // check minimum (should be first value because bin values are always
190  // assigned to the left boundary)
191  BOOST_CHECK(mapper_xyz.getMin() == minima_xyz);
192  // check maximum (should be last value + 1 step because bin values are
193  // always assigned to the left boundary)
194  BOOST_CHECK(mapper_xyz.getMax() == maxima_xyz);
195 
196  Acts::Vector3D pos0(1.2, 1.3, 1.4);
197  Acts::Vector3D pos1(1.2, 1.3, -1.4);
198  Acts::Vector3D pos2(-1.2, 1.3, 1.4);
199  Acts::Vector3D pos3(1.2, -1.3, 1.4);
200  Acts::Vector3D pos4(-1.2, -1.3, 1.4);
201 
202  auto value0_rz = mapper_rz.getField(pos0);
203  auto value1_rz = mapper_rz.getField(pos1);
204  auto value2_rz = mapper_rz.getField(pos2);
205  auto value3_rz = mapper_rz.getField(pos3);
206  auto value4_rz = mapper_rz.getField(pos4);
207 
208  auto value0_xyz = mapper_xyz.getField(pos0);
209  auto value1_xyz = mapper_xyz.getField(pos1);
210  auto value2_xyz = mapper_xyz.getField(pos2);
211  auto value3_xyz = mapper_xyz.getField(pos3);
212  auto value4_xyz = mapper_xyz.getField(pos4);
213 
214  // check z- and phi-symmetry
215  CHECK_CLOSE_REL(perp(value0_rz), perp(value1_rz), 1e-10);
216  CHECK_CLOSE_REL(value0_rz.z(), value1_rz.z(), 1e-10);
217  CHECK_CLOSE_REL(perp(value0_rz), perp(value2_rz), 1e-10);
218  CHECK_CLOSE_REL(value0_rz.z(), value2_rz.z(), 1e-10);
219  CHECK_CLOSE_REL(perp(value0_rz), perp(value3_rz), 1e-10);
220  CHECK_CLOSE_REL(value0_rz.z(), value3_rz.z(), 1e-10);
221  CHECK_CLOSE_REL(perp(value0_rz), perp(value4_rz), 1e-10);
222  CHECK_CLOSE_REL(value0_rz.z(), value4_rz.z(), 1e-10);
223 
224  // checkx-,y-,z-symmetry - need to check close (because of interpolation
225  // there can be small differences)
226  CHECK_CLOSE_REL(value0_xyz, value1_xyz, 1e-10);
227  CHECK_CLOSE_REL(value0_xyz, value2_xyz, 1e-10);
228  CHECK_CLOSE_REL(value0_xyz, value3_xyz, 1e-10);
229  CHECK_CLOSE_REL(value0_xyz, value4_xyz, 1e-10);
230 }
231 
233 BOOST_DATA_TEST_CASE(
234  bfield_symmetry_random,
235  bdata::random(
236  (bdata::seed = 0,
237  bdata::distribution = std::uniform_real_distribution<>(-10., 10.))) ^
238  bdata::random((bdata::seed = 0,
239  bdata::distribution =
240  std::uniform_real_distribution<>(-10., 10.))) ^
241  bdata::random((bdata::seed = 0,
242  bdata::distribution =
243  std::uniform_real_distribution<>(-20., 20.))) ^
244  bdata::xrange(10),
245  x, y, z, index) {
246  (void)index;
247 
248  std::vector<double> rPos;
249  std::vector<double> xPos;
250  std::vector<double> yPos;
251  std::vector<double> zPos;
252  double maxR = 20.;
253  double maxZ = 30.;
254  double maxBr = 10.;
255  double maxBz = 20.;
256  size_t nBins = 10;
257  double stepR = maxR / nBins;
258  double stepZ = maxZ / nBins;
259  double bStepR = maxBr / nBins;
260  double bStepZ = maxBz / nBins;
261 
262  for (size_t i = 0; i < nBins; i++) {
263  rPos.push_back(i * stepR);
264  xPos.push_back(i * stepR);
265  yPos.push_back(i * stepR);
266  zPos.push_back(i * stepZ);
267  }
268  // bfield in rz
269  std::vector<Acts::Vector2D> bField_rz;
270  for (size_t i = 0; i < nBins * nBins; i++) {
271  bField_rz.push_back(Acts::Vector2D(i * bStepR, i * bStepZ));
272  }
273  // the mapper in rz
274  auto mapper_rz = Acts::fieldMapperRZ(
275  [](std::array<size_t, 2> binsRZ, std::array<size_t, 2> nBinsRZ) {
276  return (binsRZ.at(1) * nBinsRZ.at(0) + binsRZ.at(0));
277  },
278  rPos, zPos, bField_rz, 1, 1, true);
279 
280  // check number of bins, minima & maxima
281  std::vector<size_t> nBins_rz = {rPos.size(), 2 * zPos.size() - 1};
282  std::vector<double> minima_rz = {0., -((nBins - 1) * stepZ)};
283  std::vector<double> maxima_rz = {nBins * stepR, nBins * stepZ};
284  BOOST_CHECK(mapper_rz.getNBins() == nBins_rz);
285  // check minimum (should be first value because bin values are always
286  // assigned to the left boundary)
287  CHECK_CLOSE_ABS(mapper_rz.getMin(), minima_rz, 1e-10);
288  // check maximum (should be last value + 1 step because bin values are
289  // always assigned to the left boundary)
290  CHECK_CLOSE_ABS(mapper_rz.getMax(), maxima_rz, 1e-10);
291 
292  // bfield in xyz
293  std::vector<Acts::Vector3D> bField_xyz;
294  for (size_t i = 0; i < nBins * nBins * nBins; i++) {
295  bField_xyz.push_back(Acts::Vector3D(i * bStepR, i * bStepR, i * bStepZ));
296  }
297  // the mapper in xyz
298  auto mapper_xyz = Acts::fieldMapperXYZ(
299  [](std::array<size_t, 3> binsXYZ, std::array<size_t, 3> nBinsXYZ) {
300  return (binsXYZ.at(0) * (nBinsXYZ.at(1) * nBinsXYZ.at(2)) +
301  binsXYZ.at(1) * nBinsXYZ.at(2) + binsXYZ.at(2));
302  },
303  xPos, yPos, zPos, bField_xyz, 1, 1, true);
304  // check number of bins, minima & maxima
305  std::vector<size_t> nBins_xyz = {2 * xPos.size() - 1, 2 * yPos.size() - 1,
306  2 * zPos.size() - 1};
307  std::vector<double> minima_xyz = {
308  -((nBins - 1) * stepR), -((nBins - 1) * stepR), -((nBins - 1) * stepZ)};
309  std::vector<double> maxima_xyz = {nBins * stepR, nBins * stepR,
310  nBins * stepZ};
311  BOOST_CHECK(mapper_xyz.getNBins() == nBins_xyz);
312  // check minimum (should be first value because bin values are always
313  // assigned to the left boundary)
314  CHECK_CLOSE_REL(mapper_xyz.getMin(), minima_xyz, 1e-10);
315  // check maximum (should be last value + 1 step because bin values are
316  // always assigned to the left boundary)
317  CHECK_CLOSE_REL(mapper_xyz.getMax(), maxima_xyz, 1e-10);
318 
319  Acts::Vector3D pos0(x, y, z);
320  Acts::Vector3D pos1(x, y, -z);
321  Acts::Vector3D pos2(-x, y, z);
322  Acts::Vector3D pos3(x, -y, z);
323  Acts::Vector3D pos4(-x, -y, z);
324 
325  auto value0_rz = mapper_rz.getField(pos0);
326  auto value1_rz = mapper_rz.getField(pos1);
327  auto value2_rz = mapper_rz.getField(pos2);
328  auto value3_rz = mapper_rz.getField(pos3);
329  auto value4_rz = mapper_rz.getField(pos4);
330 
331  // check z- and phi-symmetry
332  CHECK_CLOSE_REL(perp(value0_rz), perp(value1_rz), 1e-10);
333  CHECK_CLOSE_REL(value0_rz.z(), value1_rz.z(), 1e-10);
334  CHECK_CLOSE_REL(perp(value0_rz), perp(value2_rz), 1e-10);
335  CHECK_CLOSE_REL(value0_rz.z(), value2_rz.z(), 1e-10);
336  CHECK_CLOSE_REL(perp(value0_rz), perp(value3_rz), 1e-10);
337  CHECK_CLOSE_REL(value0_rz.z(), value3_rz.z(), 1e-10);
338  CHECK_CLOSE_REL(perp(value0_rz), perp(value4_rz), 1e-10);
339  CHECK_CLOSE_REL(value0_rz.z(), value4_rz.z(), 1e-10);
340 
341  auto value0_xyz = mapper_xyz.getField(pos0);
342  auto value1_xyz = mapper_xyz.getField(pos1);
343  auto value2_xyz = mapper_xyz.getField(pos2);
344  auto value3_xyz = mapper_xyz.getField(pos3);
345  auto value4_xyz = mapper_xyz.getField(pos4);
346 
347  // checkx-,y-,z-symmetry - need to check close (because of interpolation
348  // there can be small differences)
349  CHECK_CLOSE_REL(value0_xyz, value1_xyz, 1e-10);
350  CHECK_CLOSE_REL(value0_xyz, value2_xyz, 1e-10);
351  CHECK_CLOSE_REL(value0_xyz, value3_xyz, 1e-10);
352  CHECK_CLOSE_REL(value0_xyz, value4_xyz, 1e-10);
353 }
354 } // namespace Test
355 } // namespace Acts