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