Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BFieldMapUtils.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file BFieldMapUtils.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017-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 "Acts/MagneticField/BFieldMapUtils.hpp"
10 
11 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
12 #include "Acts/MagneticField/SolenoidBField.hpp"
13 #include "Acts/Utilities/Result.hpp"
14 #include "Acts/Utilities/VectorHelpers.hpp"
15 #include "Acts/Utilities/detail/Axis.hpp"
16 #include "Acts/Utilities/detail/Grid.hpp"
17 #include "Acts/Utilities/detail/grid_helper.hpp"
18 
19 #include <algorithm>
20 #include <cmath>
21 #include <cstdlib>
22 #include <initializer_list>
23 #include <limits>
24 #include <set>
25 #include <tuple>
26 
27 using Acts::VectorHelpers::perp;
28 using Acts::VectorHelpers::phi;
29 
30 Acts::InterpolatedBFieldMap<
31  Acts::detail::Grid<Acts::Vector2, Acts::detail::EquidistantAxis,
32  Acts::detail::EquidistantAxis>>
33 Acts::fieldMapRZ(const std::function<size_t(std::array<size_t, 2> binsRZ,
34  std::array<size_t, 2> nBinsRZ)>&
35  localToGlobalBin,
36  std::vector<double> rPos, std::vector<double> zPos,
37  std::vector<Acts::Vector2> bField, double lengthUnit,
38  double BFieldUnit, bool firstQuadrant) {
39  // [1] Create Grid
40  // sort the values
41  std::sort(rPos.begin(), rPos.end());
42  std::sort(zPos.begin(), zPos.end());
43  // Get unique values
44  rPos.erase(std::unique(rPos.begin(), rPos.end()), rPos.end());
45  zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end());
46  rPos.shrink_to_fit();
47  zPos.shrink_to_fit();
48  // get the number of bins
49  size_t nBinsR = rPos.size();
50  size_t nBinsZ = zPos.size();
51 
52  // get the minimum and maximum. We just sorted the vectors, so these are just
53  // the first and last elements.
54  double rMin = rPos[0];
55  double zMin = zPos[0];
56  double rMax = rPos[nBinsR - 1];
57  double zMax = zPos[nBinsZ - 1];
58  // calculate maxima (add one last bin, because bin value always corresponds to
59  // left boundary)
60  double stepZ = std::fabs(zMax - zMin) / (nBinsZ - 1);
61  double stepR = std::fabs(rMax - rMin) / (nBinsR - 1);
62  rMax += stepR;
63  zMax += stepZ;
64  if (firstQuadrant) {
65  zMin = -zPos[nBinsZ - 1];
66  nBinsZ = static_cast<size_t>(2. * nBinsZ - 1);
67  }
68 
69  // Create the axis for the grid
70  Acts::detail::EquidistantAxis rAxis(rMin * lengthUnit, rMax * lengthUnit,
71  nBinsR);
72  Acts::detail::EquidistantAxis zAxis(zMin * lengthUnit, zMax * lengthUnit,
73  nBinsZ);
74 
75  // Create the grid
76  using Grid_t =
77  Acts::detail::Grid<Acts::Vector2, Acts::detail::EquidistantAxis,
78  Acts::detail::EquidistantAxis>;
79  Grid_t grid(std::make_tuple(std::move(rAxis), std::move(zAxis)));
80 
81  // [2] Set the bField values
82  for (size_t i = 1; i <= nBinsR; ++i) {
83  for (size_t j = 1; j <= nBinsZ; ++j) {
84  std::array<size_t, 2> nIndices = {{rPos.size(), zPos.size()}};
85  Grid_t::index_t indices = {{i, j}};
86  if (firstQuadrant) {
87  // std::vectors begin with 0 and we do not want the user needing to
88  // take underflow or overflow bins in account this is why we need to
89  // subtract by one
90  size_t n = std::abs(int(j) - int(zPos.size()));
91  Grid_t::index_t indicesFirstQuadrant = {{i - 1, n}};
92 
93  grid.atLocalBins(indices) =
94  bField.at(localToGlobalBin(indicesFirstQuadrant, nIndices)) *
95  BFieldUnit;
96  } else {
97  // std::vectors begin with 0 and we do not want the user needing to
98  // take underflow or overflow bins in account this is why we need to
99  // subtract by one
100  grid.atLocalBins(indices) =
101  bField.at(localToGlobalBin({{i - 1, j - 1}}, nIndices)) *
102  BFieldUnit;
103  }
104  }
105  }
106  grid.setExteriorBins(Acts::Vector2::Zero());
107 
108  // [3] Create the transformation for the position
109  // map (x,y,z) -> (r,z)
110  auto transformPos = [](const Acts::Vector3& pos) {
111  return Acts::Vector2(perp(pos), pos.z());
112  };
113 
114  // [4] Create the transformation for the bfield
115  // map (Br,Bz) -> (Bx,By,Bz)
116  auto transformBField = [](const Acts::Vector2& field,
117  const Acts::Vector3& pos) {
118  double r_sin_theta_2 = pos.x() * pos.x() + pos.y() * pos.y();
119  double cos_phi = 0, sin_phi = 0;
120  if (r_sin_theta_2 > std::numeric_limits<double>::min()) {
121  double inv_r_sin_theta = 1. / sqrt(r_sin_theta_2);
122  cos_phi = pos.x() * inv_r_sin_theta;
123  sin_phi = pos.y() * inv_r_sin_theta;
124  } else {
125  cos_phi = 1.;
126  sin_phi = 0.;
127  }
128  return Acts::Vector3(field.x() * cos_phi, field.x() * sin_phi, field.y());
129  };
130 
131  // [5] Create the mapper & BField Service
132  // create field mapping
133  return Acts::InterpolatedBFieldMap<Grid_t>(
134  {transformPos, transformBField, std::move(grid)});
135 }
136 
137 Acts::InterpolatedBFieldMap<Acts::detail::Grid<
138  Acts::Vector3, Acts::detail::EquidistantAxis, Acts::detail::EquidistantAxis,
139  Acts::detail::EquidistantAxis>>
140 Acts::fieldMapXYZ(const std::function<size_t(std::array<size_t, 3> binsXYZ,
141  std::array<size_t, 3> nBinsXYZ)>&
142  localToGlobalBin,
143  std::vector<double> xPos, std::vector<double> yPos,
144  std::vector<double> zPos, std::vector<Acts::Vector3> bField,
145  double lengthUnit, double BFieldUnit, bool firstOctant) {
146  // [1] Create Grid
147  // Sort the values
148  std::sort(xPos.begin(), xPos.end());
149  std::sort(yPos.begin(), yPos.end());
150  std::sort(zPos.begin(), zPos.end());
151  // Get unique values
152  xPos.erase(std::unique(xPos.begin(), xPos.end()), xPos.end());
153  yPos.erase(std::unique(yPos.begin(), yPos.end()), yPos.end());
154  zPos.erase(std::unique(zPos.begin(), zPos.end()), zPos.end());
155  xPos.shrink_to_fit();
156  yPos.shrink_to_fit();
157  zPos.shrink_to_fit();
158  // get the number of bins
159  size_t nBinsX = xPos.size();
160  size_t nBinsY = yPos.size();
161  size_t nBinsZ = zPos.size();
162 
163  // Create the axis for the grid
164  // get minima and maximia. We just sorted the vectors, so these are just the
165  // first and last elements.
166  double xMin = xPos[0];
167  double yMin = yPos[0];
168  double zMin = zPos[0];
169  // get maxima
170  double xMax = xPos[nBinsX - 1];
171  double yMax = yPos[nBinsY - 1];
172  double zMax = zPos[nBinsZ - 1];
173  // calculate maxima (add one last bin, because bin value always corresponds to
174  // left boundary)
175  double stepZ = std::fabs(zMax - zMin) / (nBinsZ - 1);
176  double stepY = std::fabs(yMax - yMin) / (nBinsY - 1);
177  double stepX = std::fabs(xMax - xMin) / (nBinsX - 1);
178  xMax += stepX;
179  yMax += stepY;
180  zMax += stepZ;
181 
182  // If only the first octant is given
183  if (firstOctant) {
184  xMin = -xPos[nBinsX - 1];
185  yMin = -yPos[nBinsY - 1];
186  zMin = -zPos[nBinsZ - 1];
187  nBinsX = 2 * nBinsX - 1;
188  nBinsY = 2 * nBinsY - 1;
189  nBinsZ = 2 * nBinsZ - 1;
190  }
191  Acts::detail::EquidistantAxis xAxis(xMin * lengthUnit, xMax * lengthUnit,
192  nBinsX);
193  Acts::detail::EquidistantAxis yAxis(yMin * lengthUnit, yMax * lengthUnit,
194  nBinsY);
195  Acts::detail::EquidistantAxis zAxis(zMin * lengthUnit, zMax * lengthUnit,
196  nBinsZ);
197  // Create the grid
198  using Grid_t =
199  Acts::detail::Grid<Acts::Vector3, Acts::detail::EquidistantAxis,
200  Acts::detail::EquidistantAxis,
201  Acts::detail::EquidistantAxis>;
202  Grid_t grid(
203  std::make_tuple(std::move(xAxis), std::move(yAxis), std::move(zAxis)));
204 
205  // [2] Set the bField values
206  for (size_t i = 1; i <= nBinsX; ++i) {
207  for (size_t j = 1; j <= nBinsY; ++j) {
208  for (size_t k = 1; k <= nBinsZ; ++k) {
209  Grid_t::index_t indices = {{i, j, k}};
210  std::array<size_t, 3> nIndices = {
211  {xPos.size(), yPos.size(), zPos.size()}};
212  if (firstOctant) {
213  // std::vectors begin with 0 and we do not want the user needing to
214  // take underflow or overflow bins in account this is why we need to
215  // subtract by one
216  size_t m = std::abs(int(i) - (int(xPos.size())));
217  size_t n = std::abs(int(j) - (int(yPos.size())));
218  size_t l = std::abs(int(k) - (int(zPos.size())));
219  Grid_t::index_t indicesFirstOctant = {{m, n, l}};
220 
221  grid.atLocalBins(indices) =
222  bField.at(localToGlobalBin(indicesFirstOctant, nIndices)) *
223  BFieldUnit;
224 
225  } else {
226  // std::vectors begin with 0 and we do not want the user needing to
227  // take underflow or overflow bins in account this is why we need to
228  // subtract by one
229  grid.atLocalBins(indices) =
230  bField.at(localToGlobalBin({{i - 1, j - 1, k - 1}}, nIndices)) *
231  BFieldUnit;
232  }
233  }
234  }
235  }
236  grid.setExteriorBins(Acts::Vector3::Zero());
237 
238  // [3] Create the transformation for the position
239  // map (x,y,z) -> (r,z)
240  auto transformPos = [](const Acts::Vector3& pos) { return pos; };
241 
242  // [4] Create the transformation for the bfield
243  // map (Bx,By,Bz) -> (Bx,By,Bz)
244  auto transformBField = [](const Acts::Vector3& field,
245  const Acts::Vector3& /*pos*/) { return field; };
246 
247  // [5] Create the mapper & BField Service
248  // create field mapping
249  return Acts::InterpolatedBFieldMap<Grid_t>(
250  {transformPos, transformBField, std::move(grid)});
251 }
252 
253 Acts::InterpolatedBFieldMap<
254  Acts::detail::Grid<Acts::Vector2, Acts::detail::EquidistantAxis,
255  Acts::detail::EquidistantAxis>>
256 Acts::solenoidFieldMap(std::pair<double, double> rlim,
257  std::pair<double, double> zlim,
258  std::pair<size_t, size_t> nbins,
259  const SolenoidBField& field) {
260  double rMin = 0, rMax = 0, zMin = 0, zMax = 0;
261  std::tie(rMin, rMax) = rlim;
262  std::tie(zMin, zMax) = zlim;
263 
264  size_t nBinsR = 0, nBinsZ = 0;
265  std::tie(nBinsR, nBinsZ) = nbins;
266 
267  double stepZ = std::abs(zMax - zMin) / (nBinsZ - 1);
268  double stepR = std::abs(rMax - rMin) / (nBinsR - 1);
269 
270  rMax += stepR;
271  zMax += stepZ;
272 
273  // Create the axis for the grid
274  Acts::detail::EquidistantAxis rAxis(rMin, rMax, nBinsR);
275  Acts::detail::EquidistantAxis zAxis(zMin, zMax, nBinsZ);
276 
277  // Create the grid
278  using Grid_t =
279  Acts::detail::Grid<Acts::Vector2, Acts::detail::EquidistantAxis,
280  Acts::detail::EquidistantAxis>;
281  Grid_t grid(std::make_tuple(std::move(rAxis), std::move(zAxis)));
282 
283  // Create the transformation for the position
284  // map (x,y,z) -> (r,z)
285  auto transformPos = [](const Acts::Vector3& pos) {
286  return Acts::Vector2(perp(pos), pos.z());
287  };
288 
289  // Create the transformation for the bfield
290  // map (Br,Bz) -> (Bx,By,Bz)
291  auto transformBField = [](const Acts::Vector2& bfield,
292  const Acts::Vector3& pos) {
293  double r_sin_theta_2 = pos.x() * pos.x() + pos.y() * pos.y();
294  double cos_phi = 0, sin_phi = 0;
295  if (r_sin_theta_2 > std::numeric_limits<double>::min()) {
296  double inv_r_sin_theta = 1. / sqrt(r_sin_theta_2);
297  cos_phi = pos.x() * inv_r_sin_theta;
298  sin_phi = pos.y() * inv_r_sin_theta;
299  } else {
300  cos_phi = 1.;
301  sin_phi = 0.;
302  }
303  return Acts::Vector3(bfield.x() * cos_phi, bfield.x() * sin_phi,
304  bfield.y());
305  };
306 
307  // iterate over all bins, set their value to the solenoid value
308  // at their lower left position
309  for (size_t i = 0; i <= nBinsR + 1; i++) {
310  for (size_t j = 0; j <= nBinsZ + 1; j++) {
311  Grid_t::index_t index({i, j});
312  if (i == 0 || j == 0 || i == nBinsR + 1 || j == nBinsZ + 1) {
313  // under or overflow bin, set zero
314  grid.atLocalBins(index) = Grid_t::value_type(0, 0);
315  } else {
316  // regular bin, get lower left boundary
317  Grid_t::point_t lowerLeft = grid.lowerLeftBinEdge(index);
318  // do lookup
319  Vector2 B = field.getField(Vector2(lowerLeft[0], lowerLeft[1]));
320  grid.atLocalBins(index) = B;
321  }
322  }
323  }
324 
325  // Create the mapper & BField Service
326  // create field mapping
327  Acts::InterpolatedBFieldMap<Grid_t> map(
328  {transformPos, transformBField, std::move(grid)});
329  return map;
330 }