10 #include <geometry/bfieldmap/BFieldComponent3d.h>
11 #include <geometry/bfieldmap/BFieldComponentBeamline.h>
13 #include <framework/utilities/FileSystem.h>
14 #include <framework/logging/Logger.h>
16 #include <boost/iostreams/filtering_stream.hpp>
17 #include <boost/iostreams/device/file.hpp>
18 #include <boost/iostreams/filter/gzip.hpp>
24 namespace io = boost::iostreams;
26 void BFieldComponent3d::initialize()
30 if (m_mapFilename.empty()) {
31 B2ERROR(
"The filename for the 3d magnetic field component is empty !");
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 !");
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\"");
48 if ((m_exRegionR[0] == m_exRegionR[1]) || (m_exRegionZ[0] == m_exRegionZ[1])) m_exRegion =
false;
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]));
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]));
64 m_bmap.reserve(m_mapSize[0]*m_mapSize[1]*m_mapSize[2]);
66 io::filtering_istream fieldMapFile;
67 fieldMapFile.push(io::gzip_decompressor());
68 fieldMapFile.push(io::file_source(fullPath));
71 for (
int k = 0; k < m_mapSize[2]; k++) {
72 for (
int i = 0; i < m_mapSize[0]; i++) {
73 for (
int j = 0; j < m_mapSize[1]; j++) {
76 fieldMapFile.getline(tmp, 256);
79 Br = strtod(tmp + 33, &next);
80 Bphi = strtod(next, &next);
81 Bz = strtod(next,
nullptr);
82 m_bmap.emplace_back(-Br, -Bphi, -Bz);
88 if (m_errRegionR[0] != m_errRegionR[1]) {
89 auto it = m_bmap.begin();
90 for (
int k = 0; k < m_mapSize[2]; k++) {
91 double r = m_mapRegionR[0];
92 for (
int i = 0; i < m_mapSize[0]; i++, r += m_gridPitch[0]) {
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++) {
96 B.SetX(B.x() * m_errB[0]);
97 B.SetY(B.y() * m_errB[1]);
98 B.SetZ(B.z() * m_errB[2]);
104 m_igridPitch[0] = 1 / m_gridPitch[0];
105 m_igridPitch[1] = 1 / m_gridPitch[1];
106 m_igridPitch[2] = 1 / m_gridPitch[2];
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");
130 template<
int ORDER = 4> constexpr
double fast_atan2_minimax(
double y,
double x)
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;
140 }
else if (ORDER == 3) {
143 }
else if (ORDER == 4) {
146 constexpr
double c4[] = {
152 p0 = c4[0] +
v2 * c4[2];
153 p1 = c4[1] +
v2 * c4[3];
154 }
else if (ORDER == 5) {
157 constexpr
double c5[] = {
161 -0.02778537455524869,
164 p0 = c5[0] +
v2 * (c5[2] +
v2 * c5[4]);
165 p1 = c5[1] +
v2 * (c5[3]);
166 }
else if (ORDER == 6) {
169 constexpr
double c6[] = {
177 p0 = c6[0] +
v2 * (c6[2] +
v2 * (c6[4]));
178 p1 = c6[1] +
v2 * (c6[3] +
v2 * (c6[5]));
180 double phi = (z + 1) * pi4 - (z * (v - 1)) * (p0 + v * p1);
181 phi = pi2 - copysign(phi, x);
182 return copysign(phi, y);
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));
202 if (BFieldComponentBeamline::Instance().isInRange(point)) {
206 double z = point.
z();
208 if (z < m_mapRegionZ[0] || z > m_mapRegionZ[1])
return B;
210 double r2 = point.
Perp2();
212 if (r2 < m_mapRegionR[0]*m_mapRegionR[0] || r2 >= m_mapRegionR[1]*m_mapRegionR[1])
return B;
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;
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));
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));
234 double ay = std::abs(point.
y());
236 unsigned int iphi = getPhiIndexWeight(ay, point.
x(), wphi);
239 B2Vector3D b = interpolate(ir, iphi, iz, wr, wphi, wz);
241 double s = ay * norm, c = point.
x() * norm;
243 const double sgny = (point.
y() >= 0) - (point.
y() < 0);
245 B.SetXYZ(-(b.x() * c - b.y() * s), -sgny * (b.x() * s + b.y() * c), b.z());
248 B = interpolate(0, 0, iz, 0, 0, wz);
254 void BFieldComponent3d::terminate()
258 B2Vector3D BFieldComponent3d::interpolate(
unsigned int ir,
unsigned int iphi,
unsigned int iz,
259 double wr1,
double wphi1,
double wz1)
const
261 const unsigned int strideZ = m_mapSize[0] * m_mapSize[1];
262 const unsigned int strideR = m_mapSize[1];
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;
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;