Belle II Software  release-05-01-25
BFieldComponent3d.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010-20XX Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * *
7  * This software is provided "as is" without any warranty. *
8  **************************************************************************/
9 
10 #include <geometry/bfieldmap/BFieldComponent3d.h>
11 #include <geometry/bfieldmap/BFieldComponentBeamline.h>
12 
13 #include <framework/utilities/FileSystem.h>
14 #include <framework/logging/Logger.h>
15 
16 #include <boost/iostreams/filtering_stream.hpp>
17 #include <boost/iostreams/device/file.hpp>
18 #include <boost/iostreams/filter/gzip.hpp>
19 
20 #include <cmath>
21 
22 using namespace std;
23 using namespace Belle2;
24 namespace io = boost::iostreams;
25 
26 void BFieldComponent3d::initialize()
27 {
28 
29  // Input field map
30  if (m_mapFilename.empty()) {
31  B2ERROR("The filename for the 3d magnetic field component is empty !");
32  return;
33  }
34  string fullPath = FileSystem::findFile("/data/" + m_mapFilename);
35  if (!FileSystem::fileExists(fullPath)) {
36  B2ERROR("The 3d magnetic field map file '" << m_mapFilename << "' could not be found !");
37  return;
38  }
39 
40  // Options to reduce 3D to 2D
41  if (m_mapEnable != "rphiz" && m_mapEnable != "rphi" && m_mapEnable != "phiz" && m_mapEnable != "rz") {
42  B2ERROR("BField3d:: enabled coordinates must be \"rphiz\", \"rphi\", \"phiz\" or \"rz\"");
43  return;
44  }
45 
46  // Excluded region
47  m_exRegion = true;
48  if ((m_exRegionR[0] == m_exRegionR[1]) || (m_exRegionZ[0] == m_exRegionZ[1])) m_exRegion = false;
49 
50  // Print initial input parameters
51  B2DEBUG(100, "BField3d:: initial input parameters");
52  B2DEBUG(100, Form(" map filename: %s", m_mapFilename.c_str()));
53  B2DEBUG(100, Form(" map dimension: %s", m_mapEnable.c_str()));
54  if (m_interpolate) { B2DEBUG(100, Form(" map interpolation: on")); }
55  else { B2DEBUG(100, Form(" map interpolation: off")); }
56  B2DEBUG(100, Form(" map r pitch & range: %.2e [%.2e, %.2e] cm", m_gridPitch[0], m_mapRegionR[0], m_mapRegionR[1]));
57  B2DEBUG(100, Form(" map phi pitch: %.2e deg", m_gridPitch[1] * 180 / M_PI));
58  B2DEBUG(100, Form(" map z pitch & range: %.2e [%.2e, %.2e] cm", m_gridPitch[2], m_mapRegionZ[0], m_mapRegionZ[1]));
59  if (m_exRegion) {
60  B2DEBUG(100, Form(" map r excluded region: [%.2e, %.2e] cm", m_exRegionR[0], m_exRegionR[1]));
61  B2DEBUG(100, Form(" map z excluded region: [%.2e, %.2e] cm", m_exRegionZ[0], m_exRegionZ[1]));
62  }
63 
64  m_bmap.reserve(m_mapSize[0]*m_mapSize[1]*m_mapSize[2]);
65  // Load B-field map file
66  io::filtering_istream fieldMapFile;
67  fieldMapFile.push(io::gzip_decompressor());
68  fieldMapFile.push(io::file_source(fullPath));
69 
70  char tmp[256];
71  for (int k = 0; k < m_mapSize[2]; k++) { // z
72  for (int i = 0; i < m_mapSize[0]; i++) { // r
73  for (int j = 0; j < m_mapSize[1]; j++) { // phi
74  double Br, Bz, Bphi;
75  //r[m] phi[deg] z[m] Br[T] Bphi[T] Bz[T]
76  fieldMapFile.getline(tmp, 256);
77  // sscanf(tmp+33,"%lf %lf %lf",&Br,&Bphi,&Bz);
78  char* next;
79  Br = strtod(tmp + 33, &next);
80  Bphi = strtod(next, &next);
81  Bz = strtod(next, nullptr);
82  m_bmap.emplace_back(-Br, -Bphi, -Bz);
83  }
84  }
85  }
86 
87  // Introduce error on B field
88  if (m_errRegionR[0] != m_errRegionR[1]) {
89  auto it = m_bmap.begin();
90  for (int k = 0; k < m_mapSize[2]; k++) { // z
91  double r = m_mapRegionR[0];
92  for (int i = 0; i < m_mapSize[0]; i++, r += m_gridPitch[0]) { // r
93  if (!(r >= m_errRegionR[0] && r < m_errRegionR[1])) { it += m_mapSize[1]; continue;}
94  for (int j = 0; j < m_mapSize[1]; j++) { // phi
95  B2Vector3F& B = *it;
96  B.SetX(B.x() * m_errB[0]);
97  B.SetY(B.y() * m_errB[1]);
98  B.SetZ(B.z() * m_errB[2]);
99  }
100  }
101  }
102  }
103 
104  m_igridPitch[0] = 1 / m_gridPitch[0];
105  m_igridPitch[1] = 1 / m_gridPitch[1];
106  m_igridPitch[2] = 1 / m_gridPitch[2];
107 
108  B2DEBUG(100, Form("BField3d:: final map region & pitch: r [%.2e,%.2e] %.2e, phi %.2e, z [%.2e,%.2e] %.2e",
109  m_mapRegionR[0], m_mapRegionR[1], m_gridPitch[0], m_gridPitch[1],
110  m_mapRegionZ[0], m_mapRegionZ[1], m_gridPitch[2]));
111  B2DEBUG(100, "Memory consumption: " << m_bmap.size()*sizeof(B2Vector3F) / (1024 * 1024.) << " Mb");
112 }
113 
114 namespace {
130  template<int ORDER = 4> constexpr double fast_atan2_minimax(double y, double x)
131  {
132  static_assert(ORDER >= 1 && ORDER <= 6, "fast_atan2_minimax: Only orders 1-6 are supported");
133  constexpr double pi4 = M_PI / 4, pi2 = M_PI / 2;
134  const double ax = std::abs(x), ay = std::abs(y);
135  const double z = (ax - ay) / (ay + ax);
136  const double v = std::abs(z), v2 = v * v;
137  double p0{0}, p1{0};
138  if (ORDER == 2) {
139  p0 = 0.273;
140  } else if (ORDER == 3) {
141  p0 = 0.2447;
142  p1 = 0.0663;
143  } else if (ORDER == 4) {
144  // 4th order polynomial minimax approximation
145  // max error <=2.57373e-05 rad (0.00147464 deg)
146  constexpr double c4[] = {
147  +0.2132675884368258,
148  +0.23481662556227,
149  -0.2121597649928347,
150  +0.0485854027042442
151  };
152  p0 = c4[0] + v2 * c4[2];
153  p1 = c4[1] + v2 * c4[3];
154  } else if (ORDER == 5) {
155  // 5th order polynomial minimax approximation
156  // max error <=7.57429e-06 rad (0.000433975 deg)
157  constexpr double c5[] = {
158  +0.2141439315495107,
159  +0.2227491783659812,
160  -0.1628994411740733,
161  -0.02778537455524869,
162  +0.03962954416153075
163  };
164  p0 = c5[0] + v2 * (c5[2] + v2 * c5[4]);
165  p1 = c5[1] + v2 * (c5[3]);
166  } else if (ORDER == 6) {
167  // 6th order polynomial minimax approximation
168  // max error <=4.65134e-07 rad (2.66502e-05 deg)
169  constexpr double c6[] = {
170  +0.2145843590230225,
171  +0.2146820003230985,
172  -0.116250549812964,
173  -0.1428509550362758,
174  +0.1660612278047719,
175  -0.05086851503449636
176  };
177  p0 = c6[0] + v2 * (c6[2] + v2 * (c6[4]));
178  p1 = c6[1] + v2 * (c6[3] + v2 * (c6[5]));
179  }
180  double phi = (z + 1) * pi4 - (z * (v - 1)) * (p0 + v * p1);
181  phi = pi2 - copysign(phi, x);
182  return copysign(phi, y);
183  }
184 }
185 
186 B2Vector3D BFieldComponent3d::calculate(const B2Vector3D& point) const
187 {
188  auto getPhiIndexWeight = [this](double y, double x, double & wphi) -> unsigned int {
189  double phi = fast_atan2_minimax<4>(y, x);
190  wphi = phi * m_igridPitch[1];
191  auto iphi = static_cast<unsigned int>(wphi);
192  iphi = min(iphi, static_cast<unsigned int>(m_mapSize[1] - 2));
193  wphi -= iphi;
194  return iphi;
195  };
196 
197  B2Vector3D B(0, 0, 0);
198 
199  // If both '3d' and 'Beamline' components are defined in xml file,
200  // '3d' component returns zero field where 'Beamline' component is defined.
201  // If no 'Beamline' component is defined in xml file, the following function will never be called.
202  if (BFieldComponentBeamline::Instance().isInRange(point)) {
203  return B;
204  }
205 
206  double z = point.z();
207  // Check if the point lies inside the magnetic field boundaries
208  if (z < m_mapRegionZ[0] || z > m_mapRegionZ[1]) return B;
209 
210  double r2 = point.Perp2();
211  // Check if the point lies inside the magnetic field boundaries
212  if (r2 < m_mapRegionR[0]*m_mapRegionR[0] || r2 >= m_mapRegionR[1]*m_mapRegionR[1]) return B;
213  // Check if the point lies in the exclude region
214  if (m_exRegion && (z >= m_exRegionZ[0]) && (z < m_exRegionZ[1]) &&
215  (r2 >= m_exRegionR[0]*m_exRegionR[0]) && (r2 < m_exRegionR[1]*m_exRegionR[1])) return B;
216 
217  // Calculate the lower index of the point in the Z grid
218  // Note that z coordinate is inverted to match ANSYS frame
219  double wz = (m_mapRegionZ[1] - z) * m_igridPitch[2];
220  unsigned int iz = static_cast<int>(wz);
221  iz = min(iz, static_cast<unsigned int>(m_mapSize[2] - 2));
222  wz -= iz;
223 
224  if (r2 > 0) {
225  double r = sqrt(r2);
226 
227  // Calculate the lower index of the point in the R grid
228  double wr = (r - m_mapRegionR[0]) * m_igridPitch[0];
229  auto ir = static_cast<unsigned int>(wr);
230  ir = min(ir, static_cast<unsigned int>(m_mapSize[0] - 2));
231  wr -= ir;
232 
233  // Calculate the lower index of the point in the Phi grid
234  double ay = std::abs(point.y());
235  double wphi;
236  unsigned int iphi = getPhiIndexWeight(ay, point.x(), wphi);
237 
238  // Get B-field values from map
239  B2Vector3D b = interpolate(ir, iphi, iz, wr, wphi, wz); // in cylindrical system
240  double norm = 1 / r;
241  double s = ay * norm, c = point.x() * norm;
242  // Flip sign of By if y<0
243  const double sgny = (point.y() >= 0) - (point.y() < 0);
244  // in cartesian system
245  B.SetXYZ(-(b.x() * c - b.y() * s), -sgny * (b.x() * s + b.y() * c), b.z());
246  } else {
247  // Get B-field values from map in cartesian system assuming phi=0, so Bx = Br and By = Bphi
248  B = interpolate(0, 0, iz, 0, 0, wz);
249  }
250 
251  return B;
252 }
253 
254 void BFieldComponent3d::terminate()
255 {
256 }
257 
258 B2Vector3D BFieldComponent3d::interpolate(unsigned int ir, unsigned int iphi, unsigned int iz,
259  double wr1, double wphi1, double wz1) const
260 {
261  const unsigned int strideZ = m_mapSize[0] * m_mapSize[1];
262  const unsigned int strideR = m_mapSize[1];
263 
264  const double wz0 = 1 - wz1, wr0 = 1 - wr1, wphi0 = 1 - wphi1;
265  const unsigned int j000 = iz * strideZ + ir * strideR + iphi;
266  const unsigned int j001 = j000 + 1;
267  const unsigned int j010 = j000 + strideR;
268  const unsigned int j011 = j001 + strideR;
269  const unsigned int j100 = j000 + strideZ;
270  const unsigned int j101 = j001 + strideZ;
271  const unsigned int j110 = j010 + strideZ;
272  const unsigned int j111 = j011 + strideZ;
273  const double w00 = wphi0 * wr0;
274  const double w10 = wphi0 * wr1;
275  const double w01 = wphi1 * wr0;
276  const double w11 = wphi1 * wr1;
277  const vector<B2Vector3F>& B = m_bmap;
278  return
279  (B[j000] * w00 + B[j001] * w01 + B[j010] * w10 + B[j011] * w11) * wz0 +
280  (B[j100] * w00 + B[j101] * w01 + B[j110] * w10 + B[j111] * w11) * wz1;
281 }
VXDTFFilterTest::v2
const std::vector< double > v2
MATLAB generated random vector.
Definition: decorrelationMatrix.cc:44
Belle2::B2Vector3::x
DataType x() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:424
Belle2::B2Vector3< float >
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::B2Vector3::y
DataType y() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:426
Belle2::B2Vector3::Perp2
DataType Perp2() const
The transverse component squared (R^2 in cylindrical coordinate system).
Definition: B2Vector3.h:195
Belle2::B2Vector3::z
DataType z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:428