Belle II Software  release-05-01-25
BFieldComponentBeamline.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010-2011 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Andreas Moll, Kazutaka Sumisawa, Alexei Sibidanov *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
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 #include <limits>
22 
23 using namespace std;
24 namespace io = boost::iostreams;
25 namespace Belle2 {
32  struct triangle_t {
34  short int j0{0};
36  short int j1{0};
38  short int j2{0};
40  short int n0{0};
42  short int n1{0};
44  short int n2{0};
45  };
46 
48  struct xy_t {
49  double x{0};
50  double y{0};
51  };
52 
64  public:
66  [[nodiscard]] const vector<xy_t>& getPoints() const { return m_points;}
68  [[nodiscard]] const vector<triangle_t>& getTriangles() const { return m_triangles;}
69 
71  TriangularInterpolation() = default;
72 
74  TriangularInterpolation(vector<xy_t>& pc, vector<triangle_t>& ts, double d)
75  {
76  init(pc, ts, d);
77  }
78 
80  ~TriangularInterpolation() = default;
81 
90  void init(vector<xy_t>& points, vector<triangle_t>& triangles, double d)
91  {
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);
100  };
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]);
105  }
106  double xw = xmax - xmin;
107  double yw = ymax - ymin;
108 
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);
115  }
116 
122  [[nodiscard]] xy_t triangleCenter(const triangle_t& p) const
123  {
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)};
126  }
127 
138  void makeIndex(int nx, double xmin, double xmax, int ny, double ymin, double ymax)
139  {
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);
143  // cout<<nx<<" "<<xmin<<" "<<xmax<<" "<<ny<<" "<<ymin<<" "<<ymax<<endl;
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]);
147 
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);
154  };
155  for (unsigned int i = 0; i < m_triangles.size(); i++) m_triangleAreas[i] = getTriangleArea(m_triangles[i]);
156 
157 
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;}
167  }
168  int k = iy + ny * ix;
169  m_spatialIndex[k] = imin;
170  }
171  }
172  }
173 
183  [[nodiscard]] short int sideCross(short int prev, short int curr, const xy_t& r, const xy_t& v0) const
184  {
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));
194  };
195 
196  const triangle_t& p = m_triangles[curr];
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();
202  }
203 
213  void weights(short int i, const xy_t& r, double& w0, double& w1, double& w2) const
214  {
215  const triangle_t& p = m_triangles[i];
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];
222  w2 = 1 - (w0 + w1);
223  }
224 
234  [[nodiscard]] short int findTriangle(const xy_t& r0) const
235  {
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;
241  do {
242  short int next = sideCross(prev, curr, r, r0);
243  if (next == end) break;
244  prev = curr;
245  curr = next;
246  } while (true);
247  return curr;
248  }
249  protected:
251  vector<triangle_t> m_triangles;
253  vector<xy_t> m_points;
255  vector<xy_t> m_triangleCenters;
257  vector<double> m_triangleAreas;
259  vector<short int> m_spatialIndex;
261  double m_xmin;
263  double m_xmax;
265  double m_ymin;
267  double m_ymax;
269  unsigned int m_nx;
271  unsigned int m_ny;
273  double m_ixnorm{1};
275  double m_iynorm{1};
276  };
277 
288  public:
294  [[nodiscard]] const TriangularInterpolation& getTriangularInterpolation() const {return m_triInterpol;}
298  BeamlineFieldMapInterpolation() = default;
302  ~BeamlineFieldMapInterpolation() = default;
303 
312  void init(const string& fieldmapname, const string& interpolname, double validRadius)
313  {
314  if (fieldmapname.empty()) {
315  B2ERROR("The filename for the beamline magnetic field component is empty !");
316  return;
317  }
318  std::string l_fieldmapname = FileSystem::findFile("/data/" + fieldmapname);
319 
320  if (interpolname.empty()) {
321  B2ERROR("The filename for the triangulation of the beamline magnetic field component is empty !");
322  return;
323  }
324  std::string l_interpolname = FileSystem::findFile("/data/" + interpolname);
325 
326  m_rmax = validRadius;
327 
328  B2DEBUG(50, "Delaunay triangulation of the beamline field map: " << l_interpolname);
329  B2DEBUG(50, "Beamline field map: " << l_fieldmapname);
330 
331  // load interpolation triangles
332  ifstream INd(l_interpolname);
333  int nts; INd >> nts;
334  B2DEBUG(50, "Total number of triangles: " << nts);
335  vector<triangle_t> ts;
336  ts.reserve(nts);
337  // load triangle definitions from file
338  {
339  triangle_t p;
340  while (INd >> nts >> p.j0 >> p.j1 >> p.j2 >> p.n0 >> p.n1 >> p.n2) ts.push_back(p);
341  }
342 
343  //Load the map file
344  io::filtering_istream IN;
345  IN.push(io::gzip_decompressor());
346  IN.push(io::file_source(l_fieldmapname));
347 
348  int nrphi;
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;
352  m_rj2 = m_rj * m_rj;
353  m_idz0 = 1 / m_dz0;
354  m_idz1 = 1 / m_dz1;
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;
358  m_nz2 = m_nz1 + nz0;
359  // cout<<_zmax<<" "<<nz0<<" "<<m_nz1<<" "<<m_nz2<<endl;
360 
361  struct cs_t {double c, s;};
362  vector<cs_t> cs(nrphi);
363  vector<xy_t> pc;
364  vector<B2Vector3D> tbc;
365  pc.reserve(nrphi);
366  char cbuf[256]; IN.getline(cbuf, 256);
367  double rmax = 0;
368  for (int j = 0; j < nrphi; j++) {
369  IN.getline(cbuf, 256);
370  char* next = cbuf;
371  double r = strtod(next, &next);
372  double phi = strtod(next, &next);
373  strtod(next, &next);
374  double Br = strtod(next, &next);
375  double Bphi = strtod(next, &next);
376  double Bz = strtod(next, nullptr);
377  r *= 100;
378  rmax = std::max(r, rmax);
379  if (phi == 0) {
380  cs[j] = { 1, 0};
381  } else if (phi == 180) {
382  cs[j] = { -1, 0};
383  } else {
384  phi *= M_PI / 180;
385  cs[j] = {cos(phi), sin(phi)};
386  }
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);
393  }
394  // cout<<"Field map has data points up to R="<<rmax<<" cm."<<endl;
395 
396  m_idr = m_nr / (rmax - m_rj);
397  m_rmax = std::min(m_rmax, rmax);
398  bool reduce = m_rj > m_rmax;
399 
400  // if valid radius is within the triangular mesh try to reduce
401  // field map keeping only points which participate to the
402  // interpolation
403  vector<bool> ip;
404  if (reduce) {
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;
409  };
410 
411  for (int i = 0, imax = ts.size(); i < imax; i++) {
412  const triangle_t& p = ts[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)) {
415  it[i] = true;
416  ip[p.j0] = true;
417  ip[p.j1] = true;
418  ip[p.j2] = true;
419  }
420  }
421 
422  vector<short int> pindx(nrphi, -1);
423  int rnp = 0;
424  for (int i = 0, imax = ip.size(); i < imax; i++) {
425  if (ip[i]) pindx[i] = rnp++;
426  }
427  vector<xy_t> rpc;
428  rpc.reserve(rnp);
429  for (int i = 0, imax = pc.size(); i < imax; i++) {
430  if (ip[i]) rpc.push_back(pc[i]);
431  }
432 
433  vector<short int> tindx(ts.size(), -1);
434  short int rnt = 0;
435  for (int i = 0, imax = it.size(); i < imax; i++) {
436  if (it[i]) tindx[i] = rnt++;
437  }
438  vector<triangle_t> rts;
439  rts.reserve(rnt);
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++) {
443  if (it[i]) {
444  const triangle_t& t = ts[i];
445  rts.push_back({pindx[t.j0], pindx[t.j1], pindx[t.j2], newind(t.n0), newind(t.n1), newind(t.n2)});
446  }
447  }
448 
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);
451  std::swap(rpc, pc);
452  std::swap(rts, ts);
453  } else {
454  ip.resize(nrphi, true);
455  }
456  m_nxy = pc.size();
457 
458  m_triInterpol.init(pc, ts, 0.1);
459 
460  vector<B2Vector3F> bc(m_nxy * m_nz);
461  unsigned int count = 0;
462  for (int i = 0; i < nrphi; i++) {
463  if (ip[i]) bc[count++] = B2Vector3F(tbc[i]);
464  }
465 
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;
470  char* next = cbuf;
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);
481  }
482  }
483  assert(count == bc.size());
484  swap(bc, m_B);
485  }
486 
496  int zIndexAndWeight(double z, double& w1) const
497  {
498  if (std::abs(z) > m_zmax) return -1;
499  double fz;
500  int iz = 0;
501  if (z < - m_zj) {
502  fz = (z + m_zmax) * m_idz1;
503  } else if (z < m_zj) {
504  fz = (z + m_zj) * m_idz0;
505  iz = m_nz1;
506  } else {
507  fz = (z - m_zj) * m_idz1;
508  iz = m_nz2;
509  }
510  int jz = static_cast<int>(fz);
511  w1 = fz - jz;
512  iz += jz;
513  if (iz == m_nz) {
514  --iz;
515  w1 = 1;
516  }
517  return iz;
518  }
519 
529  [[nodiscard]] bool inRange(const B2Vector3D& v) const
530  {
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;
534  return true;
535  }
536 
546  [[nodiscard]] B2Vector3D interpolateField(const B2Vector3D& v) const
547  {
548  B2Vector3D res = {0, 0, 0};
549  double R2 = v.x() * v.x() + v.y() * v.y();
550  if (R2 > m_rmax * m_rmax) return res;
551  double wz1;
552  int iz = zIndexAndWeight(v.z(), wz1);
553  if (iz < 0) return res;
554  double wz0 = 1 - wz1;
555 
556  if (R2 < m_rj2) { // triangular interpolation
557  xy_t xy = {v.x(), std::abs(v.y())};
558  double w0, w1, w2;
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;
565  B += m_nxy; // next z-slice
566  b += (B[j0] * w0 + B[j1] * w1 + B[j2] * w2) * wz1;
567  res = b;
568  } else {// r-phi grid
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;
572 
573  int ir = static_cast<int>(fr);
574  int iphi = static_cast<int>(fphi);
575 
576  ir -= (ir == m_nr);
577  iphi -= (iphi == m_nphi);
578 
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;
583  int j01 = j00 + 1;
584  int j10 = j00 + nr1;
585  int j11 = j01 + nr1;
586 
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;
590  B += m_nxy; // next z-slice
591  b += (B[j00] * w00 + B[j01] * w01 + B[j10] * w10 + B[j11] * w11) * wz1;
592  res = b;
593  }
594  if (v.y() < 0) res.SetY(-res.y());
595  return res;
596  }
597  protected:
599  vector<B2Vector3F> m_B;
603  int m_nxy{0};
605  int m_nz{0};
607  int m_nz1{0};
609  int m_nz2{0};
611  int m_nr{0};
613  int m_nphi{0};
615  double m_rj{0};
617  double m_rj2{0};
619  double m_zj{0};
621  double m_dz0{0};
623  double m_dz1{0};
625  double m_idz0{0};
627  double m_idz1{0};
629  double m_idphi{0};
631  double m_idr{0};
633  double m_rmax{0};
635  double m_zmax{0};
636  };
637 
638  void BFieldComponentBeamline::initialize()
639  {
640  if (!m_ler) m_ler = new BeamlineFieldMapInterpolation;
641  if (!m_her) m_her = new BeamlineFieldMapInterpolation;
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]);
644  }
645 
646  BFieldComponentBeamline::~BFieldComponentBeamline()
647  {
648  if (m_ler) delete m_ler;
649  if (m_her) delete m_her;
650  }
651 
652  bool BFieldComponentBeamline::isInRange(const B2Vector3D& p) const
653  {
654  if (!m_ler || !m_her) return false;
655  double s = m_sinBeamCrossAngle, c = m_cosBeamCrossAngle;
656  B2Vector3D v = -p; // invert coordinates to match ANSYS one
657  double xc = v.x() * c, zs = v.z() * s, zc = v.z() * c, xs = v.x() * s;
658  B2Vector3D hv{xc - zs, v.y(), zc + xs};
659  B2Vector3D lv{xc + zs, v.y(), zc - xs};
660  return m_ler->inRange(lv) || m_her->inRange(hv);
661  }
662 
663  B2Vector3D BFieldComponentBeamline::calculate(const B2Vector3D& p) const
664  {
665  B2Vector3D res;
666  double s = m_sinBeamCrossAngle, c = m_cosBeamCrossAngle;
667  B2Vector3D v = -p; // invert coordinates to match ANSYS one
668  double xc = v.x() * c, zs = v.z() * s, zc = v.z() * c, xs = v.x() * s;
669  B2Vector3D hv{xc - zs, v.y(), zc + xs};
670  B2Vector3D lv{xc + zs, v.y(), zc - xs};
671  B2Vector3D hb = m_her->interpolateField(hv);
672  B2Vector3D lb = m_ler->interpolateField(lv);
673  B2Vector3D rhb{hb.x()* c + hb.z()* s, hb.y(), hb.z()* c - hb.x()* s};
674  B2Vector3D rlb{lb.x()* c - lb.z()* s, lb.y(), lb.z()* c + lb.x()* s};
675 
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());
678 
679  if (mhb < 1e-10) res = rlb;
680  else if (mlb < 1e-10) res = rhb;
681  else {
682  res = 0.5 * (rlb + rhb);
683  }
684  return res;
685  }
686 
687  void BFieldComponentBeamline::terminate()
688  {
689  }
690 
693  {
694  static BFieldComponentBeamline* gInstance = nullptr;
695  return &gInstance;
696  }
697 
699  BFieldComponentBeamline& BFieldComponentBeamline::Instance()
700  {
701  auto gInstance = GetInstancePtr();
702  if (*gInstance == nullptr) {
703  // Constructor creates a new instance, inits gInstance.
705  }
706  return **gInstance;
707  }
708 
709  BFieldComponentBeamline::BFieldComponentBeamline()
710  {
711  auto gInstance = GetInstancePtr();
712  if (*gInstance != nullptr) {
713  B2WARNING("BFieldComponentBeamline: object already instantiated");
714  } else {
715  *gInstance = this;
716  }
717  }
718 
720 }
Belle2::BeamlineFieldMapInterpolation::m_B
vector< B2Vector3F > m_B
Buffer for the magnetic field map.
Definition: BFieldComponentBeamline.cc:599
Belle2::TriangularInterpolation::triangleCenter
xy_t triangleCenter(const triangle_t &p) const
Calculate triangle center.
Definition: BFieldComponentBeamline.cc:122
Belle2::TriangularInterpolation::getPoints
const vector< xy_t > & getPoints() const
returns list of verticies
Definition: BFieldComponentBeamline.cc:66
Belle2::xy_t::x
double x
x component
Definition: BFieldComponentBeamline.cc:49
Belle2::TriangularInterpolation::m_triangleAreas
vector< double > m_triangleAreas
Triangle areas.
Definition: BFieldComponentBeamline.cc:257
Belle2::BeamlineFieldMapInterpolation::init
void init(const string &fieldmapname, const string &interpolname, double validRadius)
Initializes the magnetic field component.
Definition: BFieldComponentBeamline.cc:312
Belle2::xy_t
A simple 2d vector stucture.
Definition: BFieldComponentBeamline.cc:48
Belle2::BeamlineFieldMapInterpolation::zIndexAndWeight
int zIndexAndWeight(double z, double &w1) const
For a given Z coordinate calculate the index of Z slice and corresponding weight.
Definition: BFieldComponentBeamline.cc:496
Belle2::TriangularInterpolation::m_xmax
double m_xmax
Border of the region where the spatial index is constructed.
Definition: BFieldComponentBeamline.cc:263
Belle2::B2Vector3::x
DataType x() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:424
Belle2::BeamlineFieldMapInterpolation::interpolateField
B2Vector3D interpolateField(const B2Vector3D &v) const
Interpolate the magnetic field vector at the specified space point.
Definition: BFieldComponentBeamline.cc:546
Belle2::BeamlineFieldMapInterpolation::inRange
bool inRange(const B2Vector3D &v) const
Check the space point if the interpolation exists.
Definition: BFieldComponentBeamline.cc:529
Belle2::TriangularInterpolation::init
void init(vector< xy_t > &points, vector< triangle_t > &triangles, double d)
Calculate extents of a triangular mesh and build spatial index.
Definition: BFieldComponentBeamline.cc:90
Belle2::TriangularInterpolation::findTriangle
short int findTriangle(const xy_t &r0) const
Find the triangle which contain the point.
Definition: BFieldComponentBeamline.cc:234
Belle2::TriangularInterpolation::sideCross
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.
Definition: BFieldComponentBeamline.cc:183
Belle2::TriangularInterpolation::m_points
vector< xy_t > m_points
Vertex list.
Definition: BFieldComponentBeamline.cc:253
Belle2::TriangularInterpolation::TriangularInterpolation
TriangularInterpolation(vector< xy_t > &pc, vector< triangle_t > &ts, double d)
More complex constructor.
Definition: BFieldComponentBeamline.cc:74
Belle2::GetInstancePtr
BFieldComponentBeamline ** GetInstancePtr()
Static function holding the instance.
Definition: BFieldComponentBeamline.cc:692
Belle2::B2Vector3< double >
Belle2::BeamlineFieldMapInterpolation::m_triInterpol
TriangularInterpolation m_triInterpol
Object to locate point in a triangular mesh.
Definition: BFieldComponentBeamline.cc:601
Belle2::TriangularInterpolation::m_ymax
double m_ymax
Border of the region where the spatial index is constructed.
Definition: BFieldComponentBeamline.cc:267
Belle2::BeamlineFieldMapInterpolation::getTriangularInterpolation
const TriangularInterpolation & getTriangularInterpolation() const
Expose the triangular interpolation to outside.
Definition: BFieldComponentBeamline.cc:294
Belle2::BeamlineFieldMapInterpolation
The BeamlineFieldMapInterpolation class.
Definition: BFieldComponentBeamline.cc:287
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TriangularInterpolation::m_ymin
double m_ymin
Border of the region where the spatial index is constructed.
Definition: BFieldComponentBeamline.cc:265
Belle2::TriangularInterpolation::getTriangles
const vector< triangle_t > & getTriangles() const
returns list of triangles
Definition: BFieldComponentBeamline.cc:68
Belle2::TriangularInterpolation::m_xmin
double m_xmin
Border of the region where the spatial index is constructed.
Definition: BFieldComponentBeamline.cc:261
Belle2::TriangularInterpolation
The TriangularInterpolation class.
Definition: BFieldComponentBeamline.cc:63
Belle2::B2Vector3::y
DataType y() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:426
Belle2::TriangularInterpolation::makeIndex
void makeIndex(int nx, double xmin, double xmax, int ny, double ymin, double ymax)
Make spatial index.
Definition: BFieldComponentBeamline.cc:138
Belle2::ECL::triangle_t
struct for a triangle
Definition: BelleLathe.h:53
Belle2::TriangularInterpolation::m_nx
unsigned int m_nx
Spatial index grid size.
Definition: BFieldComponentBeamline.cc:269
Belle2::TriangularInterpolation::m_triangles
vector< triangle_t > m_triangles
Triangle list.
Definition: BFieldComponentBeamline.cc:251
Belle2::B2Vector3::z
DataType z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:428
Belle2::TriangularInterpolation::m_spatialIndex
vector< short int > m_spatialIndex
Spatial index.
Definition: BFieldComponentBeamline.cc:259
Belle2::TriangularInterpolation::weights
void weights(short int i, const xy_t &r, double &w0, double &w1, double &w2) const
Calculate barycentric coordinates of a point inside triangle.
Definition: BFieldComponentBeamline.cc:213
Belle2::TriangularInterpolation::m_ny
unsigned int m_ny
Spatial index grid size.
Definition: BFieldComponentBeamline.cc:271
Belle2::xy_t::y
double y
y component
Definition: BFieldComponentBeamline.cc:50
Belle2::B2Vector3F
B2Vector3< float > B2Vector3F
typedef for common usage with float
Definition: B2Vector3.h:510
Belle2::TriangularInterpolation::m_triangleCenters
vector< xy_t > m_triangleCenters
Triangle centers.
Definition: BFieldComponentBeamline.cc:255
Belle2::BFieldComponentBeamline
The BFieldComponentBeamline class.
Definition: BFieldComponentBeamline.h:42