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