ATLAS Offline Software
BFieldMesh.icc
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "CxxUtils/vec.h"
6 // constructor
7 template<class T>
8 BFieldMesh<T>::BFieldMesh(double zmin,
9  double zmax,
10  double rmin,
11  double rmax,
12  double phimin,
13  double phimax,
14  double bscale)
15  : m_scale(bscale)
16  , m_nomScale(bscale)
17 {
18  m_min = { zmin, rmin, phimin };
19  m_max = { zmax, rmax, phimax };
20 }
21 // set ranges
22 template<class T>
23 void
24 BFieldMesh<T>::setRange(double zmin,
25  double zmax,
26  double rmin,
27  double rmax,
28  double phimin,
29  double phimax)
30 {
31  m_min = { zmin, rmin, phimin };
32  m_max = { zmax, rmax, phimax };
33 }
34 
35 // set bscale
36 template<class T>
37 void
38 BFieldMesh<T>::setBscale(double bscale)
39 {
40  m_scale = m_nomScale = bscale;
41 }
42 
43 // scale bscale by a factor
44 template<class T>
45 void
46 BFieldMesh<T>::scaleBscale(double factor)
47 {
48  m_scale = factor * m_nomScale;
49 }
50 
51 //
52 // Reserve space in the vectors to avoid unnecessary memory re-allocations.
53 //
54 template<class T>
55 void
56 BFieldMesh<T>::reserve(int nz, int nr, int nphi, int nfield)
57 {
58  m_mesh[0].reserve(nz);
59  m_mesh[1].reserve(nr);
60  m_mesh[2].reserve(nphi);
61  m_field.reserve(nfield);
62 }
63 
64 template<class T>
65 void
66 BFieldMesh<T>::reserve(int nz, int nr, int nphi)
67 {
68  reserve(nz, nr, nphi, nz * nr * nphi);
69 }
70 
71 // add elements to vectors
72 template<class T>
73 void
74 BFieldMesh<T>::appendMesh(int axis, double value)
75 {
76  m_mesh[axis].push_back(value);
77 }
78 
79 template<class T>
80 void
81 BFieldMesh<T>::appendField(const BFieldVector<T>& field)
82 {
83  m_field.push_back(field);
84 }
85 
86 //
87 // Construct the look-up table to accelerate bin-finding.
88 //
89 template<class T>
90 void
91 BFieldMesh<T>::buildLUT()
92 {
93  for (int j = 0; j < 3; ++j) { // z, r, phi
94  // align the m_mesh edges to m_min/m_max
95  m_mesh[j].front() = m_min[j];
96  m_mesh[j].back() = m_max[j];
97  // determine the unit size, q, to be used in the LUTs
98  const double width = m_mesh[j].back() - m_mesh[j].front();
99  double q(width);
100  for (unsigned i = 0; i < m_mesh[j].size() - 1; ++i) {
101  q = std::min(q, m_mesh[j][i + 1] - m_mesh[j][i]);
102  }
103  // find the number of units in the LUT
104  int n = int(width / q) + 1;
105  q = width / (n + 0.5);
106  m_invUnit[j] = 1.0 / q; // new unit size
107  ++n;
108  int m = 0; // mesh number
109  for (int i = 0; i < n; ++i) { // LUT index
110  if (i * q + m_mesh[j].front() > m_mesh[j][m + 1]) {
111  m++;
112  }
113  m_LUT[j].push_back(m);
114  }
115  }
116  m_roff = m_mesh[2].size(); // index offset for incrementing r by 1
117  m_zoff = m_roff * m_mesh[1].size(); // index offset for incrementing z by 1
118 }
119 
120 template<class T>
121 int
122 BFieldMesh<T>::memSize() const
123 {
124  int size = 0;
125  size += sizeof(double) * 10;
126  size += sizeof(int) * 2;
127  for (int i = 0; i < 3; ++i) {
128  size += sizeof(double) * m_mesh[i].capacity();
129  size += sizeof(int) * m_LUT[i].capacity();
130  }
131  size += sizeof(BFieldVector<T>) * m_field.capacity();
132  return size;
133 }
134 
135 //
136 // Test if a point (z,r,phi) is inside this mesh region.
137 //
138 template<class T>
139 bool
140 BFieldMesh<T>::inside(double z, double r, double phi) const
141 {
142  // assume phi is in [-pi,pi].
143  // phimax() is in [0,2pi], but phimin() may be < 0 if the range crosses phi =
144  // 0. we have to test twice to get all possible cases.
145  if (phi < phimin()) {
146  phi += 2.0 * M_PI;
147  }
148  return (phi >= phimin() && phi <= phimax() && z >= zmin() && z <= zmax() &&
149  r >= rmin() && r <= rmax());
150 }
151 
152 //
153 // Find and return the cache of the bin containing (z,r,phi)
154 //
155 template<class T>
156 void
157 BFieldMesh<T>::getCache(double z,
158  double r,
159  double phi,
160  BFieldCache& cache,
161  double scaleFactor) const
162 {
163  // make sure phi is inside this zone
164  if (phi < phimin()) {
165  phi += 2.0 * M_PI;
166  }
167  // find the mesh, and relative location in the mesh
168  // z
169  const std::vector<double>& mz(m_mesh[0]);
170  int iz = static_cast<int>((z - zmin()) * m_invUnit[0]); // index to LUT
171  iz = m_LUT[0][iz]; // tentative mesh index from LUT
172  iz += (z > mz[iz + 1]);
173  // r
174  const std::vector<double>& mr(m_mesh[1]);
175  int ir = static_cast<int>((r - rmin()) * m_invUnit[1]); // index to LUT
176  ir = m_LUT[1][ir]; // tentative mesh index from LUT
177  ir += (r > mr[ir + 1]);
178  // phi
179  const std::vector<double>& mphi(m_mesh[2]);
180  int iphi = static_cast<int>((phi - phimin()) * m_invUnit[2]); // index to LUT
181  iphi = m_LUT[2][iphi]; // tentative mesh index from LUT
182  iphi += (phi > mphi[iphi + 1]);
183  // store the bin edges
184  cache.setRange(
185  mz[iz], mz[iz + 1], mr[ir], mr[ir + 1], mphi[iphi], mphi[iphi + 1]);
186  // store the B field at the 8 corners of the bin in the cache
187  const int im0 = iz * m_zoff + ir * m_roff + iphi; // index of the first corner
188  const double sf = scaleFactor;
189  // store the B scale
190  cache.setBscale(m_scale);
191  const BFieldVector<T>& fvec0 = m_field[im0];
192  const BFieldVector<T>& fvec1 = m_field[im0 + m_roff];
193  const BFieldVector<T>& fvec2 = m_field[im0 + m_zoff];
194  const BFieldVector<T>& fvec3 = m_field[im0 + m_zoff + m_roff];
195  const BFieldVector<T>& fvec4 = m_field[im0 + 1];
196  const BFieldVector<T>& fvec5 = m_field[im0 + m_roff + 1];
197  const BFieldVector<T>& fvec6 = m_field[im0 + m_zoff + 1];
198  const BFieldVector<T>& fvec7 = m_field[im0 + m_zoff + m_roff + 1];
199 
200  // In the usual case the m_field is of type short int
201  // else (special case) can be double
202  if constexpr (std::is_same<T, short>::value) {
203  CxxUtils::vec<int, 8> field1I = {fvec0[0], fvec1[0], fvec2[0], fvec3[0],
204  fvec4[0], fvec5[0], fvec6[0], fvec7[0]};
205 
206 
207  CxxUtils::vec<int, 8> field2I = {fvec0[1], fvec1[1], fvec2[1], fvec3[1],
208  fvec4[1], fvec5[1], fvec6[1], fvec7[1]};
209 
210 
211  CxxUtils::vec<int, 8> field3I = {fvec0[2], fvec1[2], fvec2[2], fvec3[2],
212  fvec4[2], fvec5[2], fvec6[2], fvec7[2]};
213 
214  CxxUtils::vec<double, 8> field1;
215  CxxUtils::vec<double, 8> field2;
216  CxxUtils::vec<double, 8> field3;
217  CxxUtils::vconvert(field1, field1I);
218  CxxUtils::vconvert(field2, field2I);
219  CxxUtils::vconvert(field3, field3I);
220 
221  cache.setField(sf * field1, sf * field2, sf * field3);
222  } else {
223  CxxUtils::vec<double, 8> field1 = {fvec0[0], fvec1[0], fvec2[0], fvec3[0],
224  fvec4[0], fvec5[0], fvec6[0], fvec7[0]};
225 
226  CxxUtils::vec<double, 8> field2 = {fvec0[1], fvec1[1], fvec2[1], fvec3[1],
227  fvec4[1], fvec5[1], fvec6[1], fvec7[1]};
228 
229  CxxUtils::vec<double, 8> field3 = {fvec0[2], fvec1[2], fvec2[2], fvec3[2],
230  fvec4[2], fvec5[2], fvec6[2], fvec7[2]};
231 
232  cache.setField(sf * field1, sf * field2, sf * field3);
233  }
234 }
235 
236 //
237 // Compute the magnetic field (and the derivatives) without caching
238 //
239 template<class T>
240 void
241 BFieldMesh<T>::getB(const double* ATH_RESTRICT xyz,
242  double* ATH_RESTRICT B,
243  double* ATH_RESTRICT deriv) const
244 {
245  // cylindrical coordinates
246  const double x = xyz[0];
247  const double y = xyz[1];
248  const double z = xyz[2];
249  const double r = sqrt(x * x + y * y);
250  double phi = std::atan2(y, x);
251  if (phi < phimin()) {
252  phi += 2.0 * M_PI;
253  }
254  // is it inside this map?
255  if (!inside(z, r, phi)) { // no
256  B[0] = B[1] = B[2] = 0.0;
257  if (deriv) {
258  for (int i = 0; i < 9; ++i) {
259  deriv[i] = 0.0;
260  }
261  }
262  return;
263  }
264  BFieldCache cache;
265  this->getCache(z, r, phi, cache, 1.0);
266  cache.getB(xyz, r, phi, B, deriv);
267 }
268 // accessors
269 template<class T>
270 double
271 BFieldMesh<T>::min(size_t axis) const
272 {
273  return m_min[axis];
274 }
275 
276 template<class T>
277 double
278 BFieldMesh<T>::max(size_t axis) const
279 {
280  return m_max[axis];
281 }
282 
283 template<class T>
284 double
285 BFieldMesh<T>::zmin() const
286 {
287  return m_min[0];
288 }
289 
290 template<class T>
291 double
292 BFieldMesh<T>::zmax() const
293 {
294  return m_max[0];
295 }
296 
297 template<class T>
298 double
299 BFieldMesh<T>::rmin() const
300 {
301  return m_min[1];
302 }
303 
304 template<class T>
305 double
306 BFieldMesh<T>::rmax() const
307 {
308  return m_max[1];
309 }
310 
311 template<class T>
312 double
313 BFieldMesh<T>::phimin() const
314 {
315  return m_min[2];
316 }
317 
318 template<class T>
319 double
320 BFieldMesh<T>::phimax() const
321 {
322  return m_max[2];
323 }
324 
325 template<class T>
326 unsigned
327 BFieldMesh<T>::nmesh(size_t axis) const
328 {
329  return m_mesh[axis].size();
330 }
331 
332 template<class T>
333 double
334 BFieldMesh<T>::mesh(size_t axis, size_t index) const
335 {
336  return m_mesh[axis][index];
337 }
338 
339 template<class T>
340 unsigned
341 BFieldMesh<T>::nfield() const
342 {
343  return m_field.size();
344 }
345 
346 template<class T>
347 const BFieldVector<T>&
348 BFieldMesh<T>::field(size_t index) const
349 {
350  return m_field[index];
351 }
352 
353 template<class T>
354 double
355 BFieldMesh<T>::bscale() const
356 {
357  return m_scale;
358 }
359 
Generated on Mon Sep 29 2025 21:07:03 for ATLAS Offline Software by   doxygen 1.8.18