Belle II Software  release-05-01-25
MagneticFieldComponent3D.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2017 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Martin Ritter *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <geometry/dbobjects/MagneticFieldComponent3D.h>
12 
13 namespace {
29  template<int ORDER = 4> constexpr double fast_atan2_minimax(double y, double x)
30  {
31  static_assert(ORDER >= 1 && ORDER <= 6, "fast_atan2_minimax: Only orders 1-6 are supported");
32  constexpr double pi4 = M_PI / 4, pi2 = M_PI / 2;
33  const double ax = std::abs(x), ay = std::abs(y);
34  const double z = (ax - ay) / (ay + ax);
35  const double v = std::abs(z), v2 = v * v;
36  double p0{0}, p1{0};
37  if (ORDER == 2) {
38  p0 = 0.273;
39  } else if (ORDER == 3) {
40  p0 = 0.2447;
41  p1 = 0.0663;
42  } else if (ORDER == 4) {
43  // 4th order polynomial minimax approximation
44  // max error <=2.57373e-05 rad (0.00147464 deg)
45  constexpr double c4[] = {
46  +0.2132675884368258,
47  +0.23481662556227,
48  -0.2121597649928347,
49  +0.0485854027042442
50  };
51  p0 = c4[0] + v2 * c4[2];
52  p1 = c4[1] + v2 * c4[3];
53  } else if (ORDER == 5) {
54  // 5th order polynomial minimax approximation
55  // max error <=7.57429e-06 rad (0.000433975 deg)
56  constexpr double c5[] = {
57  +0.2141439315495107,
58  +0.2227491783659812,
59  -0.1628994411740733,
60  -0.02778537455524869,
61  +0.03962954416153075
62  };
63  p0 = c5[0] + v2 * (c5[2] + v2 * c5[4]);
64  p1 = c5[1] + v2 * (c5[3]);
65  } else if (ORDER == 6) {
66  // 6th order polynomial minimax approximation
67  // max error <=4.65134e-07 rad (2.66502e-05 deg)
68  constexpr double c6[] = {
69  +0.2145843590230225,
70  +0.2146820003230985,
71  -0.116250549812964,
72  -0.1428509550362758,
73  +0.1660612278047719,
74  -0.05086851503449636
75  };
76  p0 = c6[0] + v2 * (c6[2] + v2 * (c6[4]));
77  p1 = c6[1] + v2 * (c6[3] + v2 * (c6[5]));
78  }
79  double phi = (z + 1) * pi4 - (z * (v - 1)) * (p0 + v * p1);
80  phi = pi2 - copysign(phi, x);
81  return copysign(phi, y);
82  }
83 }
84 
85 
86 namespace Belle2 {
92  {
93  using std::get;
94 
108  auto getIndexWeight = [this](double value, int coordinate) -> std::tuple<unsigned int, double> {
109  double weight = value * m_invgridPitch[coordinate];
110  auto index = static_cast<unsigned int>(weight);
111  index = std::min(index, static_cast<unsigned int>(m_mapSize[coordinate] - 2));
112  weight -= index;
113  return std::make_tuple(index, weight);
114  };
115 
116  const double z = pos.z();
117  const double r2 = pos.Perp2();
118 
119  // Calculate the lower index of the point in the Z grid
120  // Note that z coordinate is inverted to match ANSYS frame
121  const auto iz = getIndexWeight(m_maxZ - z, 2);
122 
123  // directly on z axis: get B-field values from map in cartesian system assuming phi=0, so Bx = Br and By = Bphi
124  if (r2 == 0) return interpolate(0, 0, get<0>(iz), 0, 0, get<1>(iz));
125 
126  const double r = std::sqrt(r2);
127  const auto ir = getIndexWeight(r - m_minR, 0);
128 
129  // Calculate the lower index of the point in the Phi grid
130  const double ay = std::abs(pos.y());
131  const auto iphi = getIndexWeight(fast_atan2_minimax<4>(ay, pos.x()), 1);
132 
133  // Get B-field values from map in cylindrical coordinates
134  B2Vector3D b = interpolate(get<0>(ir), get<0>(iphi), get<0>(iz), get<1>(ir), get<1>(iphi), get<1>(iz));
135  // and convert it to cartesian
136  const double norm = 1 / r;
137  const double s = ay * norm;
138  const double c = pos.x() * norm;
139  // Flip sign of By if y<0
140  const double sgny = (pos.y() >= 0) - (pos.y() < 0);
141  // in cartesian system
142  return B2Vector3D(-(b.x() * c - b.y() * s), -sgny * (b.x() * s + b.y() * c), b.z());
143  }
144 
145  B2Vector3D MagneticFieldComponent3D::interpolate(unsigned int ir, unsigned int iphi, unsigned int iz,
146  double wr1, double wphi1, double wz1) const
147  {
148  const unsigned int strideZ = m_mapSize[0] * m_mapSize[1];
149  const unsigned int strideR = m_mapSize[1];
150 
151  const double wz0 = 1 - wz1;
152  const double wr0 = 1 - wr1;
153  const double wphi0 = 1 - wphi1;
154  const unsigned int j000 = iz * strideZ + ir * strideR + iphi;
155  const unsigned int j001 = j000 + 1;
156  const unsigned int j010 = j000 + strideR;
157  const unsigned int j011 = j001 + strideR;
158  const unsigned int j100 = j000 + strideZ;
159  const unsigned int j101 = j001 + strideZ;
160  const unsigned int j110 = j010 + strideZ;
161  const unsigned int j111 = j011 + strideZ;
162  const double w00 = wphi0 * wr0;
163  const double w10 = wphi0 * wr1;
164  const double w01 = wphi1 * wr0;
165  const double w11 = wphi1 * wr1;
166  const std::vector<B2Vector3F>& B = m_bmap;
167  return
168  (B[j000] * w00 + B[j001] * w01 + B[j010] * w10 + B[j011] * w11) * wz0 +
169  (B[j100] * w00 + B[j101] * w01 + B[j110] * w10 + B[j111] * w11) * wz1;
170  }
172 }
VXDTFFilterTest::v2
const std::vector< double > v2
MATLAB generated random vector.
Definition: decorrelationMatrix.cc:44
Belle2::MagneticFieldComponent3D::interpolate
B2Vector3D interpolate(unsigned int ir, unsigned int iphi, unsigned int iz, double wr, double wphi, double wz) const
Linear interpolate the magnetic field inside a bin.
Definition: MagneticFieldComponent3D.cc:145
Belle2::MagneticFieldComponent3D::m_maxZ
float m_maxZ
maximal Z for which this field is present
Definition: MagneticFieldComponent3D.h:70
Belle2::MagneticFieldComponent3D::getField
B2Vector3D getField(const B2Vector3D &pos) const override
return the field assuming we are inside the active region as returned by inside()
Definition: MagneticFieldComponent3D.cc:91
Belle2::MagneticFieldComponent3D::m_invgridPitch
float m_invgridPitch[3]
inverted grid pitch in r, phi and z
Definition: MagneticFieldComponent3D.h:76
Belle2::MagneticFieldComponent3D::m_bmap
std::vector< B2Vector3F > m_bmap
magnetic field strength
Definition: MagneticFieldComponent3D.h:78
Belle2::B2Vector3< double >
Belle2::MagneticFieldComponent3D::m_minR
float m_minR
minimal R=sqrt(x^2+y^2) for which this field is present
Definition: MagneticFieldComponent3D.h:64
Belle2::B2Vector3D
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:507
Belle2::MagneticFieldComponent3D::m_mapSize
int m_mapSize[3]
number of bins in r, phi and z
Definition: MagneticFieldComponent3D.h:72
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19