Belle II Software light-2406-ragdoll
BeamlineFieldMapInterpolation Class Reference

The BeamlineFieldMapInterpolation class. More...

Collaboration diagram for BeamlineFieldMapInterpolation:

Public Member Functions

const TriangularInterpolationgetTriangularInterpolation () const
 Expose the triangular interpolation to outside.
 
 BeamlineFieldMapInterpolation ()=default
 Default constructor.
 
 ~BeamlineFieldMapInterpolation ()=default
 Default destructor.
 
void init (const string &fieldmapname, const string &interpolname, double validRadius)
 Initializes the magnetic field component.
 
int zIndexAndWeight (double z, double &w1) const
 For a given Z coordinate calculate the index of Z slice and corresponding weight.
 
bool inRange (const ROOT::Math::XYZVector &v) const
 Check the space point if the interpolation exists.
 
ROOT::Math::XYZVector interpolateField (const ROOT::Math::XYZVector &v) const
 Interpolate the magnetic field vector at the specified space point.
 

Protected Attributes

vector< ROOT::Math::XYZVector > m_B
 Buffer for the magnetic field map.
 
TriangularInterpolation m_triInterpol
 Object to locate point in a triangular mesh.
 
int m_nxy {0}
 Number of field points in XY plane.
 
int m_nz {0}
 Number of field slices in Z direction.
 
int m_nz1 {0}
 Start Z slice number for the finer Z grid.
 
int m_nz2 {0}
 End Z slice number for the finer Z grid.
 
int m_nr {0}
 Number of grid points in R direction.
 
int m_nphi {0}
 Number of grid points in Phi direction.
 
double m_rj {0}
 Separation radius between triangular and cylindrical meshes.
 
double m_rj2 {0}
 Square of the separation radius between triangular and cylindrical meshes.
 
double m_zj {0}
 Z border of finer Z grid.
 
double m_dz0 {0}
 Coarse Z grid pitch.
 
double m_dz1 {0}
 Finer Z grid pitch.
 
double m_idz0 {0}
 Inverse of coarse Z grid pitch.
 
double m_idz1 {0}
 Inverse of finer Z grid pitch.
 
double m_idphi {0}
 Repciprocal of Phi grid.
 
double m_idr {0}
 Repciprocal of R grid.
 
double m_rmax {0}
 Maximal radius where interpolation is still valid.
 
double m_zmax {0}
 Maximal Z where interpolation is still valid.
 

Detailed Description

The BeamlineFieldMapInterpolation class.

This class interpolates a magnetic field map around beamline. The magnetic field map is stored as a grid in cylindrical coordinates for outer radiuses and triangular mesh for inner part It is defined by a maximum radius and a maximum z value, a pitch size in both, r and z, and the number of grid points.

Definition at line 285 of file BFieldComponentBeamline.cc.

Member Function Documentation

◆ getTriangularInterpolation()

const TriangularInterpolation & getTriangularInterpolation ( ) const
inline

Expose the triangular interpolation to outside.

Returns
constant reference to the interpolation

Definition at line 292 of file BFieldComponentBeamline.cc.

292{return m_triInterpol;}
TriangularInterpolation m_triInterpol
Object to locate point in a triangular mesh.

◆ init()

void init ( const string &  fieldmapname,
const string &  interpolname,
double  validRadius 
)
inline

Initializes the magnetic field component.

This method opens the magnetic field map file and triangular mesh.

Parameters
fieldmapnameFile name containing the field map
interpolnameFile name containing triangular mesh
validRadiusMaximum radius up to which interpolation is valid

Definition at line 310 of file BFieldComponentBeamline.cc.

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 }
double m_rj
Separation radius between triangular and cylindrical meshes.
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.
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.
double m_rmax
Maximal radius where interpolation is still valid.
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
void init(vector< xy_t > &points, vector< triangle_t > &triangles, double d)
Calculate extents of a triangular mesh and build spatial index.

◆ inRange()

bool inRange ( const ROOT::Math::XYZVector &  v) const
inline

Check the space point if the interpolation exists.

Parameters
vThe space point in Cartesian coordinates (x,y,z) in [cm] at which the interpolation exists.
Returns
The magnetic field vector at the given space point in [T]. Returns false if the space point lies outside the valid region.

Definition at line 527 of file BFieldComponentBeamline.cc.

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 }

◆ interpolateField()

ROOT::Math::XYZVector interpolateField ( const ROOT::Math::XYZVector &  v) const
inline

Interpolate the magnetic field vector at the specified space point.

Parameters
vThe space point in Cartesian coordinates (x,y,z) in [cm] at which the magnetic field vector should be calculated.
Returns
The magnetic field vector at the given space point in [T]. Returns a zero vector (0,0,0) if the space point lies outside the region described by the component.

Definition at line 544 of file BFieldComponentBeamline.cc.

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 }
int zIndexAndWeight(double z, double &w1) const
For a given Z coordinate calculate the index of Z slice and corresponding weight.
const vector< triangle_t > & getTriangles() const
returns list of triangles
short int findTriangle(const xy_t &r0) const
Find the triangle which contain the point.
void weights(short int i, const xy_t &r, double &w0, double &w1, double &w2) const
Calculate barycentric coordinates of a point inside triangle.

◆ zIndexAndWeight()

int zIndexAndWeight ( double  z,
double &  w1 
) const
inline

For a given Z coordinate calculate the index of Z slice and corresponding weight.

Parameters
zZ coordinate
w1weight of the slice
Returns
the index of Z slice. Returns -1 if Z coordinate is outside valid region region.

Definition at line 494 of file BFieldComponentBeamline.cc.

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 }

Member Data Documentation

◆ m_B

vector<ROOT::Math::XYZVector> m_B
protected

Buffer for the magnetic field map.

Definition at line 597 of file BFieldComponentBeamline.cc.

◆ m_dz0

double m_dz0 {0}
protected

Coarse Z grid pitch.

Definition at line 619 of file BFieldComponentBeamline.cc.

◆ m_dz1

double m_dz1 {0}
protected

Finer Z grid pitch.

Definition at line 621 of file BFieldComponentBeamline.cc.

◆ m_idphi

double m_idphi {0}
protected

Repciprocal of Phi grid.

Definition at line 627 of file BFieldComponentBeamline.cc.

◆ m_idr

double m_idr {0}
protected

Repciprocal of R grid.

Definition at line 629 of file BFieldComponentBeamline.cc.

◆ m_idz0

double m_idz0 {0}
protected

Inverse of coarse Z grid pitch.

Definition at line 623 of file BFieldComponentBeamline.cc.

◆ m_idz1

double m_idz1 {0}
protected

Inverse of finer Z grid pitch.

Definition at line 625 of file BFieldComponentBeamline.cc.

◆ m_nphi

int m_nphi {0}
protected

Number of grid points in Phi direction.

Definition at line 611 of file BFieldComponentBeamline.cc.

◆ m_nr

int m_nr {0}
protected

Number of grid points in R direction.

Definition at line 609 of file BFieldComponentBeamline.cc.

◆ m_nxy

int m_nxy {0}
protected

Number of field points in XY plane.

Definition at line 601 of file BFieldComponentBeamline.cc.

◆ m_nz

int m_nz {0}
protected

Number of field slices in Z direction.

Definition at line 603 of file BFieldComponentBeamline.cc.

◆ m_nz1

int m_nz1 {0}
protected

Start Z slice number for the finer Z grid.

Definition at line 605 of file BFieldComponentBeamline.cc.

◆ m_nz2

int m_nz2 {0}
protected

End Z slice number for the finer Z grid.

Definition at line 607 of file BFieldComponentBeamline.cc.

◆ m_rj

double m_rj {0}
protected

Separation radius between triangular and cylindrical meshes.

Definition at line 613 of file BFieldComponentBeamline.cc.

◆ m_rj2

double m_rj2 {0}
protected

Square of the separation radius between triangular and cylindrical meshes.

Definition at line 615 of file BFieldComponentBeamline.cc.

◆ m_rmax

double m_rmax {0}
protected

Maximal radius where interpolation is still valid.

Definition at line 631 of file BFieldComponentBeamline.cc.

◆ m_triInterpol

TriangularInterpolation m_triInterpol
protected

Object to locate point in a triangular mesh.

Definition at line 599 of file BFieldComponentBeamline.cc.

◆ m_zj

double m_zj {0}
protected

Z border of finer Z grid.

Definition at line 617 of file BFieldComponentBeamline.cc.

◆ m_zmax

double m_zmax {0}
protected

Maximal Z where interpolation is still valid.

Definition at line 633 of file BFieldComponentBeamline.cc.


The documentation for this class was generated from the following file: