The BeamlineFieldMapInterpolation class. More...

Public Member Functions
const TriangularInterpolation &	getTriangularInterpolation () 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;}

Belle2::BeamlineFieldMapInterpolation::m_triInterpol

TriangularInterpolation m_triInterpol

Object to locate point in a triangular mesh.

Definition: BFieldComponentBeamline.cc:599

◆ 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

fieldmapname	File name containing the field map
interpolname	File name containing triangular mesh
validRadius	Maximum radius up to which interpolation is valid

Definition at line 310 of file BFieldComponentBeamline.cc.

    {
      if (fieldmapname.empty()) {
        B2ERROR("The filename for the beamline magnetic field component is empty !");
        return;
      }
      std::string l_fieldmapname = FileSystem::findFile("/data/" + fieldmapname);
 
      if (interpolname.empty()) {
        B2ERROR("The filename for the triangulation of the beamline magnetic field component is empty !");
        return;
      }
      std::string l_interpolname = FileSystem::findFile("/data/" + interpolname);
 
      m_rmax = validRadius;
 
      B2DEBUG(50, "Delaunay triangulation of the beamline field map: " << l_interpolname);
      B2DEBUG(50, "Beamline field map: " << l_fieldmapname);
 
      // load interpolation triangles
      ifstream INd(l_interpolname);
      int nts; INd >> nts;
      B2DEBUG(50, "Total number of triangles: " << nts);
      vector<triangle_t> ts;
      ts.reserve(nts);
      // load triangle definitions from file
      {
        triangle_t p;
        while (INd >> nts >> p.j0 >> p.j1 >> p.j2 >> p.n0 >> p.n1 >> p.n2) ts.push_back(p);
      }
 
      //Load the map file
      io::filtering_istream IN;
      IN.push(io::gzip_decompressor());
      IN.push(io::file_source(l_fieldmapname));
 
      int nrphi;
      IN >> nrphi >> m_rj >> m_nr >> m_nphi;
      IN >> m_nz >> m_zj >> m_dz0 >> m_dz1;
      m_idphi = m_nphi / M_PI;
      m_rj2 = m_rj * m_rj;
      m_idz0 = 1 / m_dz0;
      m_idz1 = 1 / m_dz1;
      int nz0 = 2 * int(m_zj / m_dz0);
      m_zmax = (m_nz - (nz0 + 1)) / 2 * m_dz1 + m_zj;
      m_nz1 = (m_nz - (nz0 + 1)) / 2;
      m_nz2 = m_nz1 + nz0;
      //    cout<<_zmax<<" "<<nz0<<" "<<m_nz1<<" "<<m_nz2<<endl;
 
      struct cs_t {double c, s;};
      vector<cs_t> cs(nrphi);
      vector<xy_t> pc;
      vector<ROOT::Math::XYZVector> tbc;
      pc.reserve(nrphi);
      char cbuf[256]; IN.getline(cbuf, 256);
      double rmax = 0;
      for (int j = 0; j < nrphi; j++) {
        IN.getline(cbuf, 256);
        char* next = cbuf;
        double r    = strtod(next, &next);
        double phi  = strtod(next, &next);
        strtod(next, &next);
        double Br   = strtod(next, &next);
        double Bphi = strtod(next, &next);
        double Bz   = strtod(next, nullptr);
        r *= 100;
        rmax = std::max(r, rmax);
        if (phi == 0) {
          cs[j] = { 1, 0};
        } else if (phi == 180) {
          cs[j] = { -1, 0};
        } else {
          phi *= M_PI / 180;
          cs[j] = {cos(phi), sin(phi)};
        }
        double x = r * cs[j].c, y = r * cs[j].s;
        pc.push_back({x, y});
        if (cs[j].s == 0) Bphi = 0;
        double Bx = Br * cs[j].c - Bphi * cs[j].s;
        double By = Br * cs[j].s + Bphi * cs[j].c;
        tbc.emplace_back(Bx, By, Bz);
      }
      //    cout<<"Field map has data points up to R="<<rmax<<" cm."<<endl;
 
      m_idr = m_nr / (rmax - m_rj);
      m_rmax = std::min(m_rmax, rmax);
      bool reduce = m_rj > m_rmax;
 
      // if valid radius is within the triangular mesh try to reduce
      // field map keeping only points which participate to the
      // interpolation
      vector<bool> ip;
      if (reduce) {
        ip.resize(nrphi, false);
        vector<bool> it(ts.size(), false);
        auto inside = [this](const xy_t & p)->bool{
          return p.x * p.x + p.y * p.y < m_rmax * m_rmax;
        };
 
        for (int i = 0, imax = ts.size(); i < imax; i++) {
          const triangle_t& p = ts[i];
          const xy_t& p0 = pc[p.j0], &p1 = pc[p.j1], &p2 = pc[p.j2];
          if (inside(p0) || inside(p1) || inside(p2)) {
            it[i] = true;
            ip[p.j0] = true;
            ip[p.j1] = true;
            ip[p.j2] = true;
          }
        }
 
        vector<short int> pindx(nrphi, -1);
        int rnp = 0;
        for (int i = 0, imax = ip.size(); i < imax; i++) {
          if (ip[i]) pindx[i] = rnp++;
        }
        vector<xy_t> rpc;
        rpc.reserve(rnp);
        for (int i = 0, imax = pc.size(); i < imax; i++) {
          if (ip[i]) rpc.push_back(pc[i]);
        }
 
        vector<short int> tindx(ts.size(), -1);
        short int rnt = 0;
        for (int i = 0, imax = it.size(); i < imax; i++) {
          if (it[i]) tindx[i] = rnt++;
        }
        vector<triangle_t> rts;
        rts.reserve(rnt);
        short int nt = ts.size();
        auto newind = [&nt, &tindx, &rnt](short int n) -> short int {return (n < nt) ? tindx[n] : rnt;};
        for (int i = 0, imax = ts.size(); i < imax; i++) {
          if (it[i]) {
            const triangle_t& t = ts[i];
            rts.push_back({pindx[t.j0], pindx[t.j1], pindx[t.j2], newind(t.n0), newind(t.n1), newind(t.n2)});
          }
        }
 
        B2DEBUG(50, "Reduce map size to cover only region R<" << m_rmax << " cm: Ntriangles=" << rnt << " Nxypoints = " << rnp <<
                " Nzslices=" << m_nz << " NBpoints = " << rnp * m_nz);
        std::swap(rpc, pc);
        std::swap(rts, ts);
      } else {
        ip.resize(nrphi, true);
      }
      m_nxy = pc.size();
 
      m_triInterpol.init(pc, ts, 0.1);
 
      vector<ROOT::Math::XYZVector> bc(m_nxy * m_nz);
      unsigned int count = 0;
      for (int i = 0; i < nrphi; i++) {
        if (ip[i]) bc[count++] = ROOT::Math::XYZVector(tbc[i]);
      }
 
      for (int i = 1; i < m_nz; ++i) {
        for (int j = 0; j < nrphi; j++) {
          IN.getline(cbuf, 256);
          if (!ip[j]) continue;
          char* next = cbuf;
          next = strchr(next, ' ');
          next = strchr(next + 1, ' ');
          next = strchr(next + 1, ' ');
          double Br   = strtod(next, &next);
          double Bphi = strtod(next, &next);
          double Bz   = strtod(next, nullptr);
          if (cs[j].s == 0) Bphi = 0;
          double Bx = Br * cs[j].c - Bphi * cs[j].s;
          double By = Br * cs[j].s + Bphi * cs[j].c;
          bc[count++].SetXYZ(Bx, By, Bz);
        }
      }
      assert(count == bc.size());
      swap(bc, m_B);
    }

◆ inRange()

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

inline

Check the space point if the interpolation exists.

Parameters

v	The 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.

    {
      if (std::abs(v.Z()) > m_zmax) return false;
      double R2 = v.X() * v.X() + v.Y() * v.Y();
      if (R2 > m_rmax * m_rmax) return false;
      return true;
    }

◆ interpolateField()

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

inline

Interpolate the magnetic field vector at the specified space point.

Parameters

v	The 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.

    {
      ROOT::Math::XYZVector res = {0, 0, 0};
      double R2 = v.X() * v.X() + v.Y() * v.Y();
      if (R2 > m_rmax * m_rmax) return res;
      double wz1;
      int iz = zIndexAndWeight(v.Z(), wz1);
      if (iz < 0) return res;
      double wz0 = 1 - wz1;
 
      if (R2 < m_rj2) { // triangular interpolation
        xy_t xy = {v.X(), std::abs(v.Y())};
        double w0, w1, w2;
        short int it = m_triInterpol.findTriangle(xy);
        m_triInterpol.weights(it, xy, w0, w1, w2);
        auto t = m_triInterpol.getTriangles().begin() + it;
        int j0 = t->j0, j1 = t->j1, j2 = t->j2;
        const ROOT::Math::XYZVector* B = m_B.data() + m_nxy * iz;
        ROOT::Math::XYZVector b = (B[j0] * w0 + B[j1] * w1 + B[j2] * w2) * wz0;
        B += m_nxy; // next z-slice
        b += (B[j0] * w0 + B[j1] * w1 + B[j2] * w2) * wz1;
        res = b;
      } else {// r-phi grid
        double r = sqrt(R2), phi = atan2(std::abs(v.Y()), v.X());
        double fr = (r - m_rj) * m_idr;
        double fphi = phi * m_idphi;
 
        int ir = static_cast<int>(fr);
        int iphi = static_cast<int>(fphi);
 
        ir -= (ir == m_nr);
        iphi -= (iphi == m_nphi);
 
        double wr1 = fr - ir, wr0 = 1 - wr1;
        double wphi1 = fphi - iphi, wphi0 = 1 - wphi1;
        const int nr1 = m_nr + 1, nphi1 = m_nphi + 1;
        int j00 = m_nxy - nr1 * (nphi1 - iphi) + ir;
        int j01 = j00 + 1;
        int j10 = j00 + nr1;
        int j11 = j01 + nr1;
 
        double w00 = wr0 * wphi0, w01 = wphi0 * wr1, w10 = wphi1 * wr0, w11 = wphi1 * wr1;
        const ROOT::Math::XYZVector* B = m_B.data() + m_nxy * iz;
        ROOT::Math::XYZVector b = (B[j00] * w00 + B[j01] * w01 + B[j10] * w10 + B[j11] * w11) * wz0;
        B += m_nxy; // next z-slice
        b += (B[j00] * w00 + B[j01] * w01 + B[j10] * w10 + B[j11] * w11) * wz1;
        res = b;
      }
      if (v.Y() < 0) res.SetY(-res.Y());
      return res;
    }

◆ zIndexAndWeight()

int zIndexAndWeight	(	double	z,
		double &	w1
	)		const

inline

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

Parameters

z	Z coordinate
w1	weight 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.

    {
      if (std::abs(z) > m_zmax) return -1;
      double fz;
      int iz = 0;
      if (z < - m_zj) {
        fz = (z + m_zmax) * m_idz1;
      } else if (z < m_zj) {
        fz = (z + m_zj) * m_idz0;
        iz = m_nz1;
      } else {
        fz = (z - m_zj) * m_idz1;
        iz = m_nz2;
      }
      int jz = static_cast<int>(fz);
      w1 = fz - jz;
      iz += jz;
      if (iz == m_nz) {
        --iz;
        w1 = 1;
      }
      return iz;
    }