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