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;
66 [[nodiscard]]
const vector<xy_t>&
getPoints()
const {
return m_points;}
68 [[nodiscard]]
const vector<triangle_t>&
getTriangles()
const {
return m_triangles;}
90 void init(vector<xy_t>& points, vector<triangle_t>& triangles,
double d)
92 std::swap(points, m_points);
93 std::swap(triangles, m_triangles);
94 const double inf = numeric_limits<double>::infinity();
95 double xmin = inf, xmax = -inf;
96 double ymin = inf, ymax = -inf;
97 auto limit = [&xmin, &xmax, &ymin, &ymax](
const xy_t & p) ->
void {
98 xmin = std::min(xmin, p.x); xmax = std::max(xmax, p.x);
99 ymin = std::min(ymin, p.y); ymax = std::max(ymax, p.y);
101 for (
const auto& t : m_triangles) {
102 limit(m_points[t.j0]);
103 limit(m_points[t.j1]);
104 limit(m_points[t.j2]);
106 double xw = xmax - xmin;
107 double yw = ymax - ymin;
109 xmin -= xw * (1. / 256), xmax += xw * (1. / 256);
110 ymin -= yw * (1. / 256), ymax += yw * (1. / 256);
111 int nx = lrint(xw * (1 + 1. / 128) / d), ny = lrint(yw * (1 + 1. / 128) / d);
112 nx = std::max(nx, 2);
113 ny = std::max(ny, 2);
114 makeIndex(nx, xmin, xmax, ny, ymin, ymax);
124 const xy_t& p0 = m_points[p.j0], &p1 = m_points[p.j1], &p2 = m_points[p.j2];
125 return {(p0.
x + p1.x + p2.x)* (1. / 3), (p0.
y + p1.y + p2.y)* (1. / 3)};
138 void makeIndex(
int nx,
double xmin,
double xmax,
int ny,
double ymin,
double ymax)
140 m_nx = nx, m_ny = ny, m_xmin = xmin, m_xmax = xmax, m_ymin = ymin, m_ymax = ymax;
141 m_ixnorm = (nx - 1) / (xmax - xmin), m_iynorm = (ny - 1) / (ymax - ymin);
142 double dx = (xmax - xmin) / (nx - 1), dy = (ymax - ymin) / (ny - 1);
144 m_spatialIndex.resize(nx * ny);
145 m_triangleCenters.resize(m_triangles.size());
146 for (
unsigned int i = 0; i < m_triangles.size(); i++) m_triangleCenters[i] = triangleCenter(m_triangles[i]);
148 m_triangleAreas.resize(m_triangles.size());
149 auto getTriangleArea = [
this](
const triangle_t& p) ->
double{
150 const xy_t& p0 = m_points[p.j0], &p1 = m_points[p.j1], &p2 = m_points[p.j2];
151 double d21x = p2.x - p1.x, d21y = p2.y - p1.y;
152 double d01x = p0.
x - p1.x, d01y = p0.
y - p1.y;
153 return 1 / (d21x * d01y - d21y * d01x);
155 for (
unsigned int i = 0; i < m_triangles.size(); i++) m_triangleAreas[i] = getTriangleArea(m_triangles[i]);
158 for (
int ix = 0; ix < nx; ix++) {
159 double x = xmin + ix * dx;
160 for (
int iy = 0; iy < ny; iy++) {
161 double y = ymin + iy * dy;
162 int imin = -1;
double dmin = 1e100;
163 for (
unsigned int i = 0; i < m_triangleCenters.size(); i++) {
164 const xy_t& p = m_triangleCenters[i];
165 double d = pow(p.x - x, 2) + pow(p.y - y, 2);
166 if (d < dmin) {imin = i; dmin = d;}
168 int k = iy + ny * ix;
169 m_spatialIndex[k] = imin;
183 [[nodiscard]]
short int sideCross(
short int prev,
short int curr,
const xy_t& r,
const xy_t& v0)
const
185 const double vx = r.x - v0.
x, vy = r.y - v0.
y;
186 auto isCrossed = [&vx, &vy, &v0](
const xy_t & p1,
const xy_t & p0) ->
bool {
187 double u0x = p0.x, u0y = p0.y;
188 double ux = p1.x - u0x, uy = p1.y - u0y;
189 double dx = u0x - v0.
x, dy = u0y - v0.
y;
190 double D = uy * vx - ux * vy;
191 double t = dx * vy - dy * vx;
192 double s = dx * uy - dy * ux;
193 return ((t < D) != (t < 0)) && ((s < D) != (s < 0));
197 const xy_t& p0 = m_points[p.j0], &p1 = m_points[p.j1], &p2 = m_points[p.j2];
198 if (p.n0 != prev && isCrossed(p1, p2))
return p.n0;
199 if (p.n1 != prev && isCrossed(p2, p0))
return p.n1;
200 if (p.n2 != prev && isCrossed(p0, p1))
return p.n2;
201 return m_triangles.size();
213 void weights(
short int i,
const xy_t& r,
double& w0,
double& w1,
double& w2)
const
216 const xy_t& p0 = m_points[p.j0], &p1 = m_points[p.j1], &p2 = m_points[p.j2];
217 double dx2 = p2.x - r.x, dy2 = p2.y - r.y;
218 double d21x = p2.x - p1.x, d21y = p2.y - p1.y;
219 double d02x = p0.
x - p2.x, d02y = p0.
y - p2.y;
220 w0 = (dx2 * d21y - dy2 * d21x) * m_triangleAreas[i];
221 w1 = (dx2 * d02y - dy2 * d02x) * m_triangleAreas[i];
236 unsigned int ix = lrint((r0.
x - m_xmin) * m_ixnorm), iy = lrint((r0.
y - m_ymin) * m_iynorm);
237 short int curr = (ix < m_nx && iy < m_ny) ? m_spatialIndex[iy + m_ny * ix] : 0;
238 xy_t r = m_triangleCenters[curr];
239 const short int end = m_triangles.size();
240 short int prev = end;
242 short int next = sideCross(prev, curr, r, r0);
243 if (next == end)
break;
312 void init(
const string& fieldmapname,
const string& interpolname,
double validRadius)
314 if (fieldmapname.empty()) {
315 B2ERROR(
"The filename for the beamline magnetic field component is empty !");
318 std::string l_fieldmapname = FileSystem::findFile(
"/data/" + fieldmapname);
320 if (interpolname.empty()) {
321 B2ERROR(
"The filename for the triangulation of the beamline magnetic field component is empty !");
324 std::string l_interpolname = FileSystem::findFile(
"/data/" + interpolname);
326 m_rmax = validRadius;
328 B2DEBUG(50,
"Delaunay triangulation of the beamline field map: " << l_interpolname);
329 B2DEBUG(50,
"Beamline field map: " << l_fieldmapname);
332 ifstream INd(l_interpolname);
334 B2DEBUG(50,
"Total number of triangles: " << nts);
335 vector<triangle_t> ts;
340 while (INd >> nts >> p.j0 >> p.j1 >> p.j2 >> p.n0 >> p.n1 >> p.n2) ts.push_back(p);
344 io::filtering_istream IN;
345 IN.push(io::gzip_decompressor());
346 IN.push(io::file_source(l_fieldmapname));
349 IN >> nrphi >> m_rj >> m_nr >> m_nphi;
350 IN >> m_nz >> m_zj >> m_dz0 >> m_dz1;
351 m_idphi = m_nphi / M_PI;
355 int nz0 = 2 * int(m_zj / m_dz0);
356 m_zmax = (m_nz - (nz0 + 1)) / 2 * m_dz1 + m_zj;
357 m_nz1 = (m_nz - (nz0 + 1)) / 2;
361 struct cs_t {
double c, s;};
362 vector<cs_t> cs(nrphi);
364 vector<B2Vector3D> tbc;
366 char cbuf[256]; IN.getline(cbuf, 256);
368 for (
int j = 0; j < nrphi; j++) {
369 IN.getline(cbuf, 256);
371 double r = strtod(next, &next);
372 double phi = strtod(next, &next);
374 double Br = strtod(next, &next);
375 double Bphi = strtod(next, &next);
376 double Bz = strtod(next,
nullptr);
378 rmax = std::max(r, rmax);
381 }
else if (phi == 180) {
385 cs[j] = {cos(phi), sin(phi)};
387 double x = r * cs[j].c, y = r * cs[j].s;
388 pc.push_back({x, y});
389 if (cs[j].s == 0) Bphi = 0;
390 double Bx = Br * cs[j].c - Bphi * cs[j].s;
391 double By = Br * cs[j].s + Bphi * cs[j].c;
392 tbc.emplace_back(Bx, By, Bz);
396 m_idr = m_nr / (rmax - m_rj);
397 m_rmax = std::min(m_rmax, rmax);
398 bool reduce = m_rj > m_rmax;
405 ip.resize(nrphi,
false);
406 vector<bool> it(ts.size(),
false);
407 auto inside = [
this](
const xy_t & p)->
bool{
408 return p.x * p.x + p.y * p.y < m_rmax * m_rmax;
411 for (
int i = 0, imax = ts.size(); i < imax; i++) {
413 const xy_t& p0 = pc[p.j0], &p1 = pc[p.j1], &p2 = pc[p.j2];
414 if (inside(p0) || inside(p1) || inside(p2)) {
422 vector<short int> pindx(nrphi, -1);
424 for (
int i = 0, imax = ip.size(); i < imax; i++) {
425 if (ip[i]) pindx[i] = rnp++;
429 for (
int i = 0, imax = pc.size(); i < imax; i++) {
430 if (ip[i]) rpc.push_back(pc[i]);
433 vector<short int> tindx(ts.size(), -1);
435 for (
int i = 0, imax = it.size(); i < imax; i++) {
436 if (it[i]) tindx[i] = rnt++;
438 vector<triangle_t> rts;
440 short int nt = ts.size();
441 auto newind = [&nt, &tindx, &rnt](
short int n) ->
short int {
return (n < nt) ? tindx[n] : rnt;};
442 for (
int i = 0, imax = ts.size(); i < imax; i++) {
445 rts.push_back({pindx[t.j0], pindx[t.j1], pindx[t.j2], newind(t.n0), newind(t.n1), newind(t.n2)});
449 B2DEBUG(50,
"Reduce map size to cover only region R<" << m_rmax <<
" cm: Ntriangles=" << rnt <<
" Nxypoints = " << rnp <<
450 " Nzslices=" << m_nz <<
" NBpoints = " << rnp * m_nz);
454 ip.resize(nrphi,
true);
458 m_triInterpol.init(pc, ts, 0.1);
460 vector<B2Vector3F> bc(m_nxy * m_nz);
461 unsigned int count = 0;
462 for (
int i = 0; i < nrphi; i++) {
466 for (
int i = 1; i < m_nz; ++i) {
467 for (
int j = 0; j < nrphi; j++) {
468 IN.getline(cbuf, 256);
469 if (!ip[j])
continue;
471 next = strchr(next,
' ');
472 next = strchr(next + 1,
' ');
473 next = strchr(next + 1,
' ');
474 double Br = strtod(next, &next);
475 double Bphi = strtod(next, &next);
476 double Bz = strtod(next,
nullptr);
477 if (cs[j].s == 0) Bphi = 0;
478 double Bx = Br * cs[j].c - Bphi * cs[j].s;
479 double By = Br * cs[j].s + Bphi * cs[j].c;
480 bc[count++].SetXYZ(Bx, By, Bz);
483 assert(count == bc.size());
498 if (std::abs(z) > m_zmax)
return -1;
502 fz = (z + m_zmax) * m_idz1;
503 }
else if (z < m_zj) {
504 fz = (z + m_zj) * m_idz0;
507 fz = (z - m_zj) * m_idz1;
510 int jz =
static_cast<int>(fz);
531 if (std::abs(v.z()) > m_zmax)
return false;
532 double R2 = v.x() * v.x() + v.y() * v.y();
533 if (R2 > m_rmax * m_rmax)
return false;
549 double R2 = v.x() * v.x() + v.y() * v.y();
550 if (R2 > m_rmax * m_rmax)
return res;
552 int iz = zIndexAndWeight(v.z(), wz1);
553 if (iz < 0)
return res;
554 double wz0 = 1 - wz1;
557 xy_t xy = {v.x(), std::abs(v.y())};
559 short int it = m_triInterpol.findTriangle(xy);
560 m_triInterpol.weights(it, xy, w0, w1, w2);
561 auto t = m_triInterpol.getTriangles().begin() + it;
562 int j0 = t->j0, j1 = t->j1, j2 = t->j2;
563 const B2Vector3F* B = m_B.data() + m_nxy * iz;
564 B2Vector3D b = (B[j0] * w0 + B[j1] * w1 + B[j2] * w2) * wz0;
566 b += (B[j0] * w0 + B[j1] * w1 + B[j2] * w2) * wz1;
569 double r = sqrt(R2), phi = atan2(std::abs(v.y()), v.x());
570 double fr = (r - m_rj) * m_idr;
571 double fphi = phi * m_idphi;
573 int ir =
static_cast<int>(fr);
574 int iphi =
static_cast<int>(fphi);
577 iphi -= (iphi == m_nphi);
579 double wr1 = fr - ir, wr0 = 1 - wr1;
580 double wphi1 = fphi - iphi, wphi0 = 1 - wphi1;
581 const int nr1 = m_nr + 1, nphi1 = m_nphi + 1;
582 int j00 = m_nxy - nr1 * (nphi1 - iphi) + ir;
587 double w00 = wr0 * wphi0, w01 = wphi0 * wr1, w10 = wphi1 * wr0, w11 = wphi1 * wr1;
588 const B2Vector3F* B = m_B.data() + m_nxy * iz;
589 B2Vector3D b = (B[j00] * w00 + B[j01] * w01 + B[j10] * w10 + B[j11] * w11) * wz0;
591 b += (B[j00] * w00 + B[j01] * w01 + B[j10] * w10 + B[j11] * w11) * wz1;
594 if (v.y() < 0) res.SetY(-res.y());
638 void BFieldComponentBeamline::initialize()
642 m_ler->
init(m_mapFilename_ler, m_interFilename_ler, m_mapRegionR[1]);
643 m_her->init(m_mapFilename_her, m_interFilename_her, m_mapRegionR[1]);
646 BFieldComponentBeamline::~BFieldComponentBeamline()
648 if (m_ler)
delete m_ler;
649 if (m_her)
delete m_her;
652 bool BFieldComponentBeamline::isInRange(
const B2Vector3D& p)
const
654 if (!m_ler || !m_her)
return false;
655 double s = m_sinBeamCrossAngle, c = m_cosBeamCrossAngle;
657 double xc = v.x() * c, zs = v.z() * s, zc = v.z() * c, xs = v.x() * s;
660 return m_ler->inRange(lv) || m_her->inRange(hv);
666 double s = m_sinBeamCrossAngle, c = m_cosBeamCrossAngle;
668 double xc = v.x() * c, zs = v.z() * s, zc = v.z() * c, xs = v.x() * s;
676 double mhb = std::abs(rhb.x()) + std::abs(rhb.y()) + std::abs(rhb.z());
677 double mlb = std::abs(rlb.x()) + std::abs(rlb.y()) + std::abs(rlb.z());
679 if (mhb < 1e-10) res = rlb;
680 else if (mlb < 1e-10) res = rhb;
682 res = 0.5 * (rlb + rhb);
687 void BFieldComponentBeamline::terminate()
702 if (*gInstance ==
nullptr) {
709 BFieldComponentBeamline::BFieldComponentBeamline()
712 if (*gInstance !=
nullptr) {
713 B2WARNING(
"BFieldComponentBeamline: object already instantiated");