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