Belle II Software light-2406-ragdoll
BFieldComponentBeamline.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/BFieldComponentBeamline.h>
10
11#include <framework/utilities/FileSystem.h>
12#include <framework/logging/Logger.h>
13
14#include <boost/iostreams/filtering_stream.hpp>
15#include <boost/iostreams/device/file.hpp>
16#include <boost/iostreams/filter/gzip.hpp>
17
18#include <cmath>
19#include <limits>
20
21using namespace std;
22namespace io = boost::iostreams;
23namespace Belle2 {
30 struct triangle_t {
32 short int j0{0};
34 short int j1{0};
36 short int j2{0};
38 short int n0{0};
40 short int n1{0};
42 short int n2{0};
43 };
44
46 struct xy_t {
47 double x{0};
48 double y{0};
49 };
50
62 public:
64 [[nodiscard]] const vector<xy_t>& getPoints() const { return m_points;}
66 [[nodiscard]] const vector<triangle_t>& getTriangles() const { return m_triangles;}
67
70
72 TriangularInterpolation(vector<xy_t>& pc, vector<triangle_t>& ts, double d)
73 {
74 init(pc, ts, d);
75 }
76
79
88 void init(vector<xy_t>& points, vector<triangle_t>& triangles, double d)
89 {
90 std::swap(points, m_points);
91 std::swap(triangles, m_triangles);
92 const double inf = numeric_limits<double>::infinity();
93 double xmin = inf, xmax = -inf;
94 double ymin = inf, ymax = -inf;
95 auto limit = [&xmin, &xmax, &ymin, &ymax](const xy_t & p) -> void {
96 xmin = std::min(xmin, p.x); xmax = std::max(xmax, p.x);
97 ymin = std::min(ymin, p.y); ymax = std::max(ymax, p.y);
98 };
99 for (const auto& t : m_triangles) {
100 limit(m_points[t.j0]);
101 limit(m_points[t.j1]);
102 limit(m_points[t.j2]);
103 }
104 double xw = xmax - xmin;
105 double yw = ymax - ymin;
106
107 xmin -= xw * (1. / 256), xmax += xw * (1. / 256);
108 ymin -= yw * (1. / 256), ymax += yw * (1. / 256);
109 int nx = lrint(xw * (1 + 1. / 128) / d), ny = lrint(yw * (1 + 1. / 128) / d);
110 nx = std::max(nx, 2);
111 ny = std::max(ny, 2);
112 makeIndex(nx, xmin, xmax, ny, ymin, ymax);
113 }
114
120 [[nodiscard]] xy_t triangleCenter(const triangle_t& p) const
121 {
122 const xy_t& p0 = m_points[p.j0], &p1 = m_points[p.j1], &p2 = m_points[p.j2];
123 return {(p0.x + p1.x + p2.x)* (1. / 3), (p0.y + p1.y + p2.y)* (1. / 3)};
124 }
125
136 void makeIndex(int nx, double xmin, double xmax, int ny, double ymin, double ymax)
137 {
138 m_nx = nx, m_ny = ny, m_xmin = xmin, m_xmax = xmax, m_ymin = ymin, m_ymax = ymax;
139 m_ixnorm = (nx - 1) / (xmax - xmin), m_iynorm = (ny - 1) / (ymax - ymin);
140 double dx = (xmax - xmin) / (nx - 1), dy = (ymax - ymin) / (ny - 1);
141 // cout<<nx<<" "<<xmin<<" "<<xmax<<" "<<ny<<" "<<ymin<<" "<<ymax<<endl;
142 m_spatialIndex.resize(nx * ny);
143 m_triangleCenters.resize(m_triangles.size());
144 for (unsigned int i = 0; i < m_triangles.size(); i++) m_triangleCenters[i] = triangleCenter(m_triangles[i]);
145
146 m_triangleAreas.resize(m_triangles.size());
147 auto getTriangleArea = [this](const triangle_t& p) ->double{
148 const xy_t& p0 = m_points[p.j0], &p1 = m_points[p.j1], &p2 = m_points[p.j2];
149 double d21x = p2.x - p1.x, d21y = p2.y - p1.y;
150 double d01x = p0.x - p1.x, d01y = p0.y - p1.y;
151 return 1 / (d21x * d01y - d21y * d01x);
152 };
153 for (unsigned int i = 0; i < m_triangles.size(); i++) m_triangleAreas[i] = getTriangleArea(m_triangles[i]);
154
155
156 for (int ix = 0; ix < nx; ix++) {
157 double x = xmin + ix * dx;
158 for (int iy = 0; iy < ny; iy++) {
159 double y = ymin + iy * dy;
160 int imin = -1; double dmin = 1e100;
161 for (unsigned int i = 0; i < m_triangleCenters.size(); i++) {
162 const xy_t& p = m_triangleCenters[i];
163 double d = pow(p.x - x, 2) + pow(p.y - y, 2);
164 if (d < dmin) {imin = i; dmin = d;}
165 }
166 int k = iy + ny * ix;
167 m_spatialIndex[k] = imin;
168 }
169 }
170 }
171
181 [[nodiscard]] short int sideCross(short int prev, short int curr, const xy_t& r, const xy_t& v0) const
182 {
183 const double vx = r.x - v0.x, vy = r.y - v0.y;
184 auto isCrossed = [&vx, &vy, &v0](const xy_t & p1, const xy_t & p0) -> bool {
185 double u0x = p0.x, u0y = p0.y;
186 double ux = p1.x - u0x, uy = p1.y - u0y;
187 double dx = u0x - v0.x, dy = u0y - v0.y;
188 double D = uy * vx - ux * vy;
189 double t = dx * vy - dy * vx;
190 double s = dx * uy - dy * ux;
191 return ((t < D) != (t < 0)) && ((s < D) != (s < 0));
192 };
193
194 const triangle_t& p = m_triangles[curr];
195 const xy_t& p0 = m_points[p.j0], &p1 = m_points[p.j1], &p2 = m_points[p.j2];
196 if (p.n0 != prev && isCrossed(p1, p2)) return p.n0;
197 if (p.n1 != prev && isCrossed(p2, p0)) return p.n1;
198 if (p.n2 != prev && isCrossed(p0, p1)) return p.n2;
199 return m_triangles.size();
200 }
201
211 void weights(short int i, const xy_t& r, double& w0, double& w1, double& w2) const
212 {
213 const triangle_t& p = m_triangles[i];
214 const xy_t& p0 = m_points[p.j0], &p1 = m_points[p.j1], &p2 = m_points[p.j2];
215 double dx2 = p2.x - r.x, dy2 = p2.y - r.y;
216 double d21x = p2.x - p1.x, d21y = p2.y - p1.y;
217 double d02x = p0.x - p2.x, d02y = p0.y - p2.y;
218 w0 = (dx2 * d21y - dy2 * d21x) * m_triangleAreas[i];
219 w1 = (dx2 * d02y - dy2 * d02x) * m_triangleAreas[i];
220 w2 = 1 - (w0 + w1);
221 }
222
232 [[nodiscard]] short int findTriangle(const xy_t& r0) const
233 {
234 unsigned int ix = lrint((r0.x - m_xmin) * m_ixnorm), iy = lrint((r0.y - m_ymin) * m_iynorm);
235 short int curr = (ix < m_nx && iy < m_ny) ? m_spatialIndex[iy + m_ny * ix] : 0;
236 xy_t r = m_triangleCenters[curr];
237 const short int end = m_triangles.size();
238 short int prev = end;
239 do {
240 short int next = sideCross(prev, curr, r, r0);
241 if (next == end) break;
242 prev = curr;
243 curr = next;
244 } while (true);
245 return curr;
246 }
247 protected:
249 vector<triangle_t> m_triangles;
251 vector<xy_t> m_points;
253 vector<xy_t> m_triangleCenters;
255 vector<double> m_triangleAreas;
257 vector<short int> m_spatialIndex;
259 double m_xmin;
261 double m_xmax;
263 double m_ymin;
265 double m_ymax;
267 unsigned int m_nx;
269 unsigned int m_ny;
271 double m_ixnorm{1};
273 double m_iynorm{1};
274 };
275
286 public:
301
310 void init(const string& fieldmapname, const string& interpolname, double validRadius)
311 {
312 if (fieldmapname.empty()) {
313 B2ERROR("The filename for the beamline magnetic field component is empty !");
314 return;
315 }
316 std::string l_fieldmapname = FileSystem::findFile("/data/" + fieldmapname);
317
318 if (interpolname.empty()) {
319 B2ERROR("The filename for the triangulation of the beamline magnetic field component is empty !");
320 return;
321 }
322 std::string l_interpolname = FileSystem::findFile("/data/" + interpolname);
323
324 m_rmax = validRadius;
325
326 B2DEBUG(50, "Delaunay triangulation of the beamline field map: " << l_interpolname);
327 B2DEBUG(50, "Beamline field map: " << l_fieldmapname);
328
329 // load interpolation triangles
330 ifstream INd(l_interpolname);
331 int nts; INd >> nts;
332 B2DEBUG(50, "Total number of triangles: " << nts);
333 vector<triangle_t> ts;
334 ts.reserve(nts);
335 // load triangle definitions from file
336 {
337 triangle_t p;
338 while (INd >> nts >> p.j0 >> p.j1 >> p.j2 >> p.n0 >> p.n1 >> p.n2) ts.push_back(p);
339 }
340
341 //Load the map file
342 io::filtering_istream IN;
343 IN.push(io::gzip_decompressor());
344 IN.push(io::file_source(l_fieldmapname));
345
346 int nrphi;
347 IN >> nrphi >> m_rj >> m_nr >> m_nphi;
348 IN >> m_nz >> m_zj >> m_dz0 >> m_dz1;
349 m_idphi = m_nphi / M_PI;
350 m_rj2 = m_rj * m_rj;
351 m_idz0 = 1 / m_dz0;
352 m_idz1 = 1 / m_dz1;
353 int nz0 = 2 * int(m_zj / m_dz0);
354 m_zmax = (m_nz - (nz0 + 1)) / 2 * m_dz1 + m_zj;
355 m_nz1 = (m_nz - (nz0 + 1)) / 2;
356 m_nz2 = m_nz1 + nz0;
357 // cout<<_zmax<<" "<<nz0<<" "<<m_nz1<<" "<<m_nz2<<endl;
358
359 struct cs_t {double c, s;};
360 vector<cs_t> cs(nrphi);
361 vector<xy_t> pc;
362 vector<ROOT::Math::XYZVector> tbc;
363 pc.reserve(nrphi);
364 char cbuf[256]; IN.getline(cbuf, 256);
365 double rmax = 0;
366 for (int j = 0; j < nrphi; j++) {
367 IN.getline(cbuf, 256);
368 char* next = cbuf;
369 double r = strtod(next, &next);
370 double phi = strtod(next, &next);
371 strtod(next, &next);
372 double Br = strtod(next, &next);
373 double Bphi = strtod(next, &next);
374 double Bz = strtod(next, nullptr);
375 r *= 100;
376 rmax = std::max(r, rmax);
377 if (phi == 0) {
378 cs[j] = { 1, 0};
379 } else if (phi == 180) {
380 cs[j] = { -1, 0};
381 } else {
382 phi *= M_PI / 180;
383 cs[j] = {cos(phi), sin(phi)};
384 }
385 double x = r * cs[j].c, y = r * cs[j].s;
386 pc.push_back({x, y});
387 if (cs[j].s == 0) Bphi = 0;
388 double Bx = Br * cs[j].c - Bphi * cs[j].s;
389 double By = Br * cs[j].s + Bphi * cs[j].c;
390 tbc.emplace_back(Bx, By, Bz);
391 }
392 // cout<<"Field map has data points up to R="<<rmax<<" cm."<<endl;
393
394 m_idr = m_nr / (rmax - m_rj);
395 m_rmax = std::min(m_rmax, rmax);
396 bool reduce = m_rj > m_rmax;
397
398 // if valid radius is within the triangular mesh try to reduce
399 // field map keeping only points which participate to the
400 // interpolation
401 vector<bool> ip;
402 if (reduce) {
403 ip.resize(nrphi, false);
404 vector<bool> it(ts.size(), false);
405 auto inside = [this](const xy_t & p)->bool{
406 return p.x * p.x + p.y * p.y < m_rmax * m_rmax;
407 };
408
409 for (int i = 0, imax = ts.size(); i < imax; i++) {
410 const triangle_t& p = ts[i];
411 const xy_t& p0 = pc[p.j0], &p1 = pc[p.j1], &p2 = pc[p.j2];
412 if (inside(p0) || inside(p1) || inside(p2)) {
413 it[i] = true;
414 ip[p.j0] = true;
415 ip[p.j1] = true;
416 ip[p.j2] = true;
417 }
418 }
419
420 vector<short int> pindx(nrphi, -1);
421 int rnp = 0;
422 for (int i = 0, imax = ip.size(); i < imax; i++) {
423 if (ip[i]) pindx[i] = rnp++;
424 }
425 vector<xy_t> rpc;
426 rpc.reserve(rnp);
427 for (int i = 0, imax = pc.size(); i < imax; i++) {
428 if (ip[i]) rpc.push_back(pc[i]);
429 }
430
431 vector<short int> tindx(ts.size(), -1);
432 short int rnt = 0;
433 for (int i = 0, imax = it.size(); i < imax; i++) {
434 if (it[i]) tindx[i] = rnt++;
435 }
436 vector<triangle_t> rts;
437 rts.reserve(rnt);
438 short int nt = ts.size();
439 auto newind = [&nt, &tindx, &rnt](short int n) -> short int {return (n < nt) ? tindx[n] : rnt;};
440 for (int i = 0, imax = ts.size(); i < imax; i++) {
441 if (it[i]) {
442 const triangle_t& t = ts[i];
443 rts.push_back({pindx[t.j0], pindx[t.j1], pindx[t.j2], newind(t.n0), newind(t.n1), newind(t.n2)});
444 }
445 }
446
447 B2DEBUG(50, "Reduce map size to cover only region R<" << m_rmax << " cm: Ntriangles=" << rnt << " Nxypoints = " << rnp <<
448 " Nzslices=" << m_nz << " NBpoints = " << rnp * m_nz);
449 std::swap(rpc, pc);
450 std::swap(rts, ts);
451 } else {
452 ip.resize(nrphi, true);
453 }
454 m_nxy = pc.size();
455
456 m_triInterpol.init(pc, ts, 0.1);
457
458 vector<ROOT::Math::XYZVector> bc(m_nxy * m_nz);
459 unsigned int count = 0;
460 for (int i = 0; i < nrphi; i++) {
461 if (ip[i]) bc[count++] = ROOT::Math::XYZVector(tbc[i]);
462 }
463
464 for (int i = 1; i < m_nz; ++i) {
465 for (int j = 0; j < nrphi; j++) {
466 IN.getline(cbuf, 256);
467 if (!ip[j]) continue;
468 char* next = cbuf;
469 next = strchr(next, ' ');
470 next = strchr(next + 1, ' ');
471 next = strchr(next + 1, ' ');
472 double Br = strtod(next, &next);
473 double Bphi = strtod(next, &next);
474 double Bz = strtod(next, nullptr);
475 if (cs[j].s == 0) Bphi = 0;
476 double Bx = Br * cs[j].c - Bphi * cs[j].s;
477 double By = Br * cs[j].s + Bphi * cs[j].c;
478 bc[count++].SetXYZ(Bx, By, Bz);
479 }
480 }
481 assert(count == bc.size());
482 swap(bc, m_B);
483 }
484
494 int zIndexAndWeight(double z, double& w1) const
495 {
496 if (std::abs(z) > m_zmax) return -1;
497 double fz;
498 int iz = 0;
499 if (z < - m_zj) {
500 fz = (z + m_zmax) * m_idz1;
501 } else if (z < m_zj) {
502 fz = (z + m_zj) * m_idz0;
503 iz = m_nz1;
504 } else {
505 fz = (z - m_zj) * m_idz1;
506 iz = m_nz2;
507 }
508 int jz = static_cast<int>(fz);
509 w1 = fz - jz;
510 iz += jz;
511 if (iz == m_nz) {
512 --iz;
513 w1 = 1;
514 }
515 return iz;
516 }
517
527 [[nodiscard]] bool inRange(const ROOT::Math::XYZVector& v) const
528 {
529 if (std::abs(v.Z()) > m_zmax) return false;
530 double R2 = v.X() * v.X() + v.Y() * v.Y();
531 if (R2 > m_rmax * m_rmax) return false;
532 return true;
533 }
534
544 [[nodiscard]] ROOT::Math::XYZVector interpolateField(const ROOT::Math::XYZVector& v) const
545 {
546 ROOT::Math::XYZVector res = {0, 0, 0};
547 double R2 = v.X() * v.X() + v.Y() * v.Y();
548 if (R2 > m_rmax * m_rmax) return res;
549 double wz1;
550 int iz = zIndexAndWeight(v.Z(), wz1);
551 if (iz < 0) return res;
552 double wz0 = 1 - wz1;
553
554 if (R2 < m_rj2) { // triangular interpolation
555 xy_t xy = {v.X(), std::abs(v.Y())};
556 double w0, w1, w2;
557 short int it = m_triInterpol.findTriangle(xy);
558 m_triInterpol.weights(it, xy, w0, w1, w2);
559 auto t = m_triInterpol.getTriangles().begin() + it;
560 int j0 = t->j0, j1 = t->j1, j2 = t->j2;
561 const ROOT::Math::XYZVector* B = m_B.data() + m_nxy * iz;
562 ROOT::Math::XYZVector b = (B[j0] * w0 + B[j1] * w1 + B[j2] * w2) * wz0;
563 B += m_nxy; // next z-slice
564 b += (B[j0] * w0 + B[j1] * w1 + B[j2] * w2) * wz1;
565 res = b;
566 } else {// r-phi grid
567 double r = sqrt(R2), phi = atan2(std::abs(v.Y()), v.X());
568 double fr = (r - m_rj) * m_idr;
569 double fphi = phi * m_idphi;
570
571 int ir = static_cast<int>(fr);
572 int iphi = static_cast<int>(fphi);
573
574 ir -= (ir == m_nr);
575 iphi -= (iphi == m_nphi);
576
577 double wr1 = fr - ir, wr0 = 1 - wr1;
578 double wphi1 = fphi - iphi, wphi0 = 1 - wphi1;
579 const int nr1 = m_nr + 1, nphi1 = m_nphi + 1;
580 int j00 = m_nxy - nr1 * (nphi1 - iphi) + ir;
581 int j01 = j00 + 1;
582 int j10 = j00 + nr1;
583 int j11 = j01 + nr1;
584
585 double w00 = wr0 * wphi0, w01 = wphi0 * wr1, w10 = wphi1 * wr0, w11 = wphi1 * wr1;
586 const ROOT::Math::XYZVector* B = m_B.data() + m_nxy * iz;
587 ROOT::Math::XYZVector b = (B[j00] * w00 + B[j01] * w01 + B[j10] * w10 + B[j11] * w11) * wz0;
588 B += m_nxy; // next z-slice
589 b += (B[j00] * w00 + B[j01] * w01 + B[j10] * w10 + B[j11] * w11) * wz1;
590 res = b;
591 }
592 if (v.Y() < 0) res.SetY(-res.Y());
593 return res;
594 }
595 protected:
597 vector<ROOT::Math::XYZVector> m_B;
601 int m_nxy{0};
603 int m_nz{0};
605 int m_nz1{0};
607 int m_nz2{0};
609 int m_nr{0};
611 int m_nphi{0};
613 double m_rj{0};
615 double m_rj2{0};
617 double m_zj{0};
619 double m_dz0{0};
621 double m_dz1{0};
623 double m_idz0{0};
625 double m_idz1{0};
627 double m_idphi{0};
629 double m_idr{0};
631 double m_rmax{0};
633 double m_zmax{0};
634 };
635
637 {
642 }
643
645 {
646 if (m_ler) delete m_ler;
647 if (m_her) delete m_her;
648 }
649
650 bool BFieldComponentBeamline::isInRange(const ROOT::Math::XYZVector& p) const
651 {
652 if (!m_ler || !m_her) return false;
654 ROOT::Math::XYZVector v = -p; // invert coordinates to match ANSYS one
655 double xc = v.X() * c, zs = v.Z() * s, zc = v.Z() * c, xs = v.X() * s;
656 ROOT::Math::XYZVector hv{xc - zs, v.Y(), zc + xs};
657 ROOT::Math::XYZVector lv{xc + zs, v.Y(), zc - xs};
658 return m_ler->inRange(lv) || m_her->inRange(hv);
659 }
660
661 ROOT::Math::XYZVector BFieldComponentBeamline::calculate(const ROOT::Math::XYZVector& p) const
662 {
663 ROOT::Math::XYZVector res;
665 ROOT::Math::XYZVector v = -p; // invert coordinates to match ANSYS one
666 double xc = v.X() * c, zs = v.Z() * s, zc = v.Z() * c, xs = v.X() * s;
667 ROOT::Math::XYZVector hv{xc - zs, v.Y(), zc + xs};
668 ROOT::Math::XYZVector lv{xc + zs, v.Y(), zc - xs};
669 ROOT::Math::XYZVector hb = m_her->interpolateField(hv);
670 ROOT::Math::XYZVector lb = m_ler->interpolateField(lv);
671 ROOT::Math::XYZVector rhb{hb.X()* c + hb.Z()* s, hb.Y(), hb.Z()* c - hb.X()* s};
672 ROOT::Math::XYZVector rlb{lb.X()* c - lb.Z()* s, lb.Y(), lb.Z()* c + lb.X()* s};
673
674 double mhb = std::abs(rhb.X()) + std::abs(rhb.Y()) + std::abs(rhb.Z());
675 double mlb = std::abs(rlb.X()) + std::abs(rlb.Y()) + std::abs(rlb.Z());
676
677 if (mhb < 1e-10) res = rlb;
678 else if (mlb < 1e-10) res = rhb;
679 else {
680 res = 0.5 * (rlb + rhb);
681 }
682 return res;
683 }
684
686 {
687 }
688
691 {
692 static BFieldComponentBeamline* gInstance = nullptr;
693 return &gInstance;
694 }
695
698 {
699 auto gInstance = GetInstancePtr();
700 if (*gInstance == nullptr) {
701 // Constructor creates a new instance, inits gInstance.
703 }
704 return **gInstance;
705 }
706
708 {
709 auto gInstance = GetInstancePtr();
710 if (*gInstance != nullptr) {
711 B2WARNING("BFieldComponentBeamline: object already instantiated");
712 } else {
713 *gInstance = this;
714 }
715 }
716
718}
The BFieldComponentBeamline class.
std::string m_mapFilename_ler
The filename of the magnetic field map.
BeamlineFieldMapInterpolation * m_ler
Actual magnetic field interpolation object for LER.
BeamlineFieldMapInterpolation * m_her
Actual magnetic field interpolation object for HER.
double m_sinBeamCrossAngle
The sin of the crossing angle of the beams
double m_cosBeamCrossAngle
The cos of the crossing angle of the beams
std::string m_interFilename_ler
The filename of the map for interpolation.
std::string m_mapFilename_her
Parameter to set Angle of the beam.
std::string m_interFilename_her
The filename of the map for interpolation.
double m_mapRegionR[2]
The min and max boundaries of the map region in r.
The BeamlineFieldMapInterpolation class.
const TriangularInterpolation & getTriangularInterpolation() const
Expose the triangular interpolation to outside.
double m_rj
Separation radius between triangular and cylindrical meshes.
bool inRange(const ROOT::Math::XYZVector &v) const
Check the space point if the interpolation exists.
int zIndexAndWeight(double z, double &w1) const
For a given Z coordinate calculate the index of Z slice and corresponding weight.
double m_idz1
Inverse of finer Z grid pitch.
int m_nphi
Number of grid points in Phi direction.
int m_nr
Number of grid points in R direction.
double m_zmax
Maximal Z where interpolation is still valid.
double m_rj2
Square of the separation radius between triangular and cylindrical meshes.
int m_nz
Number of field slices in Z direction.
ROOT::Math::XYZVector interpolateField(const ROOT::Math::XYZVector &v) const
Interpolate the magnetic field vector at the specified space point.
int m_nz2
End Z slice number for the finer Z grid.
vector< ROOT::Math::XYZVector > m_B
Buffer for the magnetic field map.
int m_nxy
Number of field points in XY plane.
double m_idz0
Inverse of coarse Z grid pitch.
TriangularInterpolation m_triInterpol
Object to locate point in a triangular mesh.
void init(const string &fieldmapname, const string &interpolname, double validRadius)
Initializes the magnetic field component.
double m_rmax
Maximal radius where interpolation is still valid.
BeamlineFieldMapInterpolation()=default
Default constructor.
~BeamlineFieldMapInterpolation()=default
Default destructor.
int m_nz1
Start Z slice number for the finer Z grid.
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
The TriangularInterpolation class.
vector< short int > m_spatialIndex
Spatial index.
double m_xmin
Border of the region where the spatial index is constructed.
unsigned int m_ny
Spatial index grid size.
double m_ixnorm
Reciprocals to speedup the index calculation.
~TriangularInterpolation()=default
Destructor.
TriangularInterpolation(vector< xy_t > &pc, vector< triangle_t > &ts, double d)
More complex constructor.
vector< triangle_t > m_triangles
Triangle list.
void init(vector< xy_t > &points, vector< triangle_t > &triangles, double d)
Calculate extents of a triangular mesh and build spatial index.
xy_t triangleCenter(const triangle_t &p) const
Calculate triangle center.
double m_ymin
Border of the region where the spatial index is constructed.
vector< xy_t > m_triangleCenters
Triangle centers.
short int sideCross(short int prev, short int curr, const xy_t &r, const xy_t &v0) const
Determine which triangle side is crossed by a line segment defined by r and v0 points.
const vector< triangle_t > & getTriangles() const
returns list of triangles
double m_xmax
Border of the region where the spatial index is constructed.
void makeIndex(int nx, double xmin, double xmax, int ny, double ymin, double ymax)
Make spatial index.
const vector< xy_t > & getPoints() const
returns list of verticies
vector< double > m_triangleAreas
Triangle areas.
unsigned int m_nx
Spatial index grid size.
short int findTriangle(const xy_t &r0) const
Find the triangle which contain the point.
TriangularInterpolation()=default
Default constructor.
void weights(short int i, const xy_t &r, double &w0, double &w1, double &w2) const
Calculate barycentric coordinates of a point inside triangle.
double m_ymax
Border of the region where the spatial index is constructed.
double m_iynorm
Reciprocals to speedup the index calculation.
virtual void initialize() override
Initializes the magnetic field component.
BFieldComponentBeamline ** GetInstancePtr()
Static function holding the instance.
bool isInRange(const ROOT::Math::XYZVector &point) const
Check presence of beamline field at the specific space point in the detector coordinate frame.
virtual void terminate() override
Terminates the magnetic field component.
virtual ROOT::Math::XYZVector calculate(const ROOT::Math::XYZVector &point) const override
Calculates the magnetic field vector at the specified space point.
virtual ~BFieldComponentBeamline()
The BFieldComponentBeamline destructor.
static BFieldComponentBeamline & Instance()
BFieldComponentBeamline instance.
BFieldComponentBeamline()
The BFieldComponentBeamline constructor.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.
short int n1
2nd adjacent triangle in a list of triangles
short int j1
2nd vertex index in a list of xy-points
short int n2
3rd adjacent triangle in a list of triangles
short int j0
1st vertex index in a list of xy-points
short int n0
1st adjacent triangle in a list of triangles
short int j2
3rd vertex index in a list of xy-points
A simple 2d vector stucture.
double y
y component
double x
x component