Belle II Software light-2406-ragdoll
BFieldComponentQuad.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/bfieldmap/BFieldComponentQuad.h>
10
11#include <framework/utilities/FileSystem.h>
12#include <framework/logging/Logger.h>
13#include <framework/gearbox/Unit.h>
14#include <framework/gearbox/Const.h>
15
16#include <boost/iostreams/filtering_stream.hpp>
17#include <boost/iostreams/device/file.hpp>
18
19#include <cmath>
20
21using namespace std;
22using namespace Belle2;
23namespace io = boost::iostreams;
24
33template <class ForwardIterator, class T, class Compare>
34ForwardIterator linear_sentinel(ForwardIterator it, const T& val, Compare comp)
35{
36 do {
37 if (comp(val, *it)) break;
38 ++it;
39 } while (true);
40 return it;
41}
42
44{
45 // check if input name is not empty
46 if (m_mapFilenameHER.empty()) {
47 B2FATAL("The filename for the HER quadrupole magnetic field component is empty !");
48 return;
49 }
50 if (m_mapFilenameLER.empty()) {
51 B2FATAL("The filename for the LER quadrupole magnetic field component is empty !");
52 return;
53 }
54 if (m_apertFilenameHER.empty()) {
55 B2FATAL("The filename for the HER beam pipe aperture is empty !");
56 return;
57 }
58 if (m_apertFilenameLER.empty()) {
59 B2FATAL("The filename for the LER beam pipe aperture is empty !");
60 return;
61 }
62
63 // check if input files exsits
64 string fullPathMapHER = FileSystem::findFile("/data/" + m_mapFilenameHER);
65 if (!FileSystem::fileExists(fullPathMapHER)) {
66 B2FATAL("The HER quadrupole magnetic field map file '" << m_mapFilenameHER << "' could not be found !");
67 return;
68 }
69 string fullPathMapLER = FileSystem::findFile("/data/" + m_mapFilenameLER);
70 if (!FileSystem::fileExists(fullPathMapLER)) {
71 B2FATAL("The LER quadrupole magnetic field map file '" << m_mapFilenameLER << "' could not be found !");
72 return;
73 }
74 string fullPathMapHERleak = FileSystem::findFile("/data/" + m_mapFilenameHERleak);
75 if (!FileSystem::fileExists(fullPathMapHERleak)) {
76 B2FATAL("The HERleak quadrupole magnetic field map file '" << m_mapFilenameHERleak << "' could not be found !");
77 return;
78 }
79 string fullPathApertHER = FileSystem::findFile("/data/" + m_apertFilenameHER);
80 if (!FileSystem::fileExists(fullPathApertHER)) {
81 B2FATAL("The HER aperture file '" << m_apertFilenameHER << "' could not be found !");
82 return;
83 }
84 string fullPathApertLER = FileSystem::findFile("/data/" + m_apertFilenameLER);
85 if (!FileSystem::fileExists(fullPathApertLER)) {
86 B2FATAL("The LER aperture file '" << m_apertFilenameLER << "' could not be found !");
87 return;
88 }
89
90 //--------------------------------------------------------------
91 // Load the field map files
92 //--------------------------------------------------------------
93
95 struct ParamPoint {
96 double s{0};
97 double L{0};
98 double K0{0};
99 double K1{0};
100 double SK0{0};
101 double SK1{0};
102 double ROTATE{0};
103 double DX{0};
104 double DY{0};
105 /* Note that K parameters used in SAD is multiplied by the element length.
106 * Popular definitions are: K0,SK0[1/m] and K1,SK1[1/m^2]
107 */
108 };
109
110 io::filtering_istream fieldMapFileHER, fieldMapFileLER, fieldMapFileHERleak;
111 fieldMapFileHER.push(io::file_source(fullPathMapHER));
112 fieldMapFileLER.push(io::file_source(fullPathMapLER));
113 fieldMapFileHERleak.push(io::file_source(fullPathMapHERleak));
114
115 //Create the parameter map and read the data from the file
116 auto readLenseParameters = [](istream & IN) -> vector<ParamPoint> {
117 vector<ParamPoint> r;
118 string name;
119 ParamPoint t;
120 while (IN >> name >> t.s >> t.L >> t.K0 >> t.K1 >> t.SK0 >> t.SK1 >> t.ROTATE >> t.DX >> t.DY)
121 {
122 /* Convert parameters in basf2 default unit [cm].*/
123 t.s *= 100;
124 t.L *= 100;
125 t.DX *= 100;
126 t.DY *= 100;
127 t.K1 /= 100;
128 t.SK1 /= 100;
129 r.push_back(t);
130 }
131 return r;
132 };
133
134 std::vector<ParamPoint> params_her, params_ler, params_herl;
135 params_her = readLenseParameters(fieldMapFileHER);
136 params_ler = readLenseParameters(fieldMapFileLER);
137 params_herl = readLenseParameters(fieldMapFileHERleak);
138
139 //--------------------------------------------------------------
140 // Load the aperture map files
141 //--------------------------------------------------------------
142 auto readAperturePoints = [this](istream & IN) -> vector<ApertPoint> {
143 vector<ApertPoint> r;
144 ApertPoint t;
145 while (IN >> t.s >> t.r)
146 {
147 /* s and r are in [mm] in Apert?ER.dat. */
148 /* Transform parameters in basf2 default unit [cm].*/
149 t.s /= 10; t.r /= 10;
150 r.push_back(t);
151 m_maxr2 = max(m_maxr2, t.r);
152 }
153 r.push_back({numeric_limits<double>::infinity(), m_maxr2});
154 return r;
155 };
156
157 io::filtering_istream apertFileHER, apertFileLER;
158 apertFileHER.push(io::file_source(fullPathApertHER));
159 apertFileLER.push(io::file_source(fullPathApertLER));
160
161 m_ah = readAperturePoints(apertFileHER);
162 m_al = readAperturePoints(apertFileLER);
163
164 m_maxr2 *= m_maxr2;
165
172 auto proc3 = [](const vector<ParamPoint>& in, double p0) -> vector<ParamPoint3> {
173 vector<ParamPoint3> out;
174 out.resize(in.size());
175 auto it = out.begin();
176 for (auto b = in.begin(); b != in.end(); ++b, ++it)
177 {
178 ParamPoint3& t = *it;
179 t.s = b->s;
180 t.L = b->L;
181
182 double K0 = b->K0;
183 double K1 = b->K1;
184 double SK0 = b->SK0;
185 double SK1 = b->SK1;
186 double DX = b->DX;
187 double DY = b->DY;
188 double k = p0 / t.L;
189 K0 *= k;
190 K1 *= k;
191 SK0 *= k;
192 SK1 *= k;
193
194#define ROTATE_DIRECTION 1 // 1: default, 0: rotate-off, -1: inversely rotate
195#if ROTATE_DIRECTION==0
196 double sphi = 0, cphi = 1;
197#else
198 double sphi, cphi;
199 sincos(-b->ROTATE * (M_PI / 180), &sphi, &cphi);
200#if ROTATE_DIRECTION==-1
201 sphi = -sphi;
202#endif
203#endif
204 double s2phi = 2 * cphi * sphi;
205 double c2phi = cphi * cphi - sphi * sphi;
206 double tp = DX * SK1 + DY * K1;
207 double tm = DY * SK1 - DX * K1;
208
209 t.mxx = K1 * s2phi + SK1 * c2phi;
210 t.mxy = -SK1 * s2phi + K1 * c2phi;
211 t.mx0 = tm * s2phi - tp * c2phi + SK0 * cphi + K0 * sphi;
212
213 t.myx = -SK1 * s2phi + K1 * c2phi;
214 t.myy = -K1 * s2phi - SK1 * c2phi;
215 t.my0 = tp * s2phi + tm * c2phi + K0 * cphi - SK0 * sphi;
216 }
217 return out;
218 };
219
220 const double c = Const::speedOfLight / (Unit::m / Unit::s) / 100;
221 const double p0_HER = 7.0e+9 / c, p0_LER = 4.0e+9 / c;
222
223 m_h3 = proc3(params_her, p0_HER);
224 vector<ParamPoint3> hleak3 = proc3(params_herl, p0_HER);
225 m_l3 = proc3(params_ler, p0_LER);
226
236 auto merge = [](vector<ParamPoint3>& v, const vector<ParamPoint3>& a) {
237 for (const auto& t : a) {
238 auto i = upper_bound(v.begin(), v.end(), t.s, [](double s, const ParamPoint3 & b) {return s < b.s;});
239 if (i == v.end()) continue;
240 ParamPoint3& p = *(i - 1);
241 if (p.s == t.s && p.L == t.L)
242 p += t;
243 else
244 v.insert(i, t);
245 }
246 };
247 merge(m_h3, hleak3);
248
255 auto getranges = [](const vector<ParamPoint3>& v) {
256 ranges_t r;
257 double s0 = v[0].s, smax = v.back().s + v.back().L;
258 for (int i = 0, N = v.size(); i < N - 1; i++) {
259 const ParamPoint3& t0 = v[i], t1 = v[i + 1];
260 double sn = t0.s + t0.L;
261 if (abs(sn - t1.s) > 1e-10) {
262 r.push_back({s0, sn});
263 s0 = t1.s;
264 }
265 }
266 r.push_back({s0, smax});
267 return r;
268 };
269
270 m_ranges_her = getranges(m_h3);
271 m_ranges_ler = getranges(m_l3);
272
273 const double inf = numeric_limits<double>::infinity();
274 // add sentinels
275 m_ranges_her.insert(m_ranges_her.begin(), { -inf, -inf});
276 m_ranges_her.insert(m_ranges_her.end(), {inf, inf});
277 m_ranges_ler.insert(m_ranges_ler.begin(), { -inf, -inf});
278 m_ranges_ler.insert(m_ranges_ler.end(), {inf, inf});
279
286 auto associate_aperture = [](const vector<ApertPoint>& ap, const ranges_t& v) {
287 auto less = [](double s, const ApertPoint & b) {return s < b.s;};
288 vector<std::vector<ApertPoint>::const_iterator> res;
289 for (auto r : v) {
290 auto i0 = upper_bound(ap.begin(), ap.end(), r.r0, less);
291 if (i0 == ap.begin() || i0 == ap.end()) continue;
292 res.push_back(i0 - 1);
293 }
294 return res;
295 };
296 m_offset_ap_her = associate_aperture(m_ah, m_ranges_her);
297 m_offset_ap_ler = associate_aperture(m_al, m_ranges_ler);
298
305 auto associate_lenses = [](const vector<ParamPoint3>& ap, const ranges_t& v) {
306 vector<std::vector<ParamPoint3>::const_iterator> res;
307 for (auto r : v) {
308 auto i0 = find_if(ap.begin(), ap.end(), [r](const ParamPoint3 & b) {return abs(b.s - r.r0) < 1e-10;});
309 if (i0 == ap.end()) continue;
310 res.push_back(i0);
311 }
312 return res;
313 };
314 m_offset_pp_her = associate_lenses(m_h3, m_ranges_her);
315 m_offset_pp_ler = associate_lenses(m_l3, m_ranges_ler);
316}
317
319{
320 int i = getRange(s, m_ranges_her);
321 if (i < 0) return 0;
322 return getAperture(s, m_offset_ap_her[i]);
323}
324
326{
327 int i = getRange(s, m_ranges_ler);
328 if (i < 0) return 0;
329 return getAperture(s, m_offset_ap_ler[i]);
330}
331
332int BFieldComponentQuad::getRange(double s, const ranges_t& r) const
333{
334 auto it0 = r.begin() + 1;
335 auto it = linear_sentinel(it0, s, [](double s_, const range_t& r_) {return s_ <= r_.r0;});
336 if (s > (--it)->r1) return -1;
337 return it - it0;
338}
339
340double BFieldComponentQuad::getAperture(double s, std::vector<ApertPoint>::const_iterator jt0) const
341{
342 auto jt = linear_sentinel(jt0, s, [](double s_, const ApertPoint & r) {return s_ <= r.s;});
343 const ApertPoint& p1 = *(jt - 1), &p2 = *jt;
344 return p1.r + (p2.r - p1.r) / (p2.s - p1.s) * (s - p1.s);
345}
346
347ROOT::Math::XYZVector BFieldComponentQuad::calculate(const ROOT::Math::XYZVector& point) const
348{
349 const double sa = sin(0.0415), ca = cos(0.0415);
350 //assume point is given in [cm]
351 ROOT::Math::XYZVector B(0, 0, 0);
352
353 const ROOT::Math::XYZVector& v{point};
354 double xc = v.X() * ca, zs = v.Z() * sa, zc = v.Z() * ca, xs = v.X() * sa;
355 ROOT::Math::XYZVector vh{(xc - zs), -v.Y(), -(zc + xs)}; // to the HER beamline frame
356 ROOT::Math::XYZVector vl{(xc + zs), -v.Y(), -(zc - xs)}; // to the LER beamline frame
357
358 double r2h = vh.Perp2(), r2l = vl.Perp2();
359
360 if (r2h < r2l) { /* the point is closer to HER*/
361 if (r2h < m_maxr2) { /* within max radius */
362 double s = vh.Z();
363 int i = getRange(s, m_ranges_her);
364 if (i < 0) return B;
365 double r = getAperture(s, m_offset_ap_her[i]);
366 if (r2h < r * r) {
367 auto kt = m_offset_pp_her[i] + static_cast<unsigned int>(s - m_ranges_her[i + 1].r0);
368 double Bx = kt->getBx(vh.X(), vh.Y());
369 double By = kt->getBy(vh.X(), vh.Y());
370 B.SetXYZ(Bx * ca, -By, -Bx * sa); // to the detector frame
371 }
372 }
373 } else { /* the point is closer to LER*/
374 if (r2l < m_maxr2) { /* within max radius */
375 double s = vl.Z();
376 int i = getRange(s, m_ranges_ler);
377 if (i < 0) return B;
378 double r = getAperture(s, m_offset_ap_ler[i]);
379 if (r2l < r * r) {
380 auto kt = m_offset_pp_ler[i] + static_cast<unsigned int>(s - m_ranges_ler[i + 1].r0);
381 double Bx = kt->getBx(vl.X(), vl.Y());
382 double By = kt->getBy(vl.X(), vl.Y());
383 B.SetXYZ(Bx * ca, -By, Bx * sa); // to the detector frame
384 }
385 }
386 }
387 return B;
388}
int getRange(double a, const ranges_t &b) const
Search for range occupied by optics since now only for ranges are present use linear search.
double getApertureLER(double s) const
Returns the LER beam pipe aperture at given position.
double m_maxr2
The square of maximal aperture for fast rejection.
virtual void initialize() override
Initializes the magnetic field component.
std::vector< ApertPoint > m_al
The the aperture parameters for LER.
std::vector< ApertPoint > m_ah
The the aperture parameters for HER.
std::vector< std::vector< ApertPoint >::const_iterator > m_offset_ap_ler
The vector of pointer to accelerate search in aperture for ler.
std::string m_mapFilenameLER
Magnetic field map of LER
std::string m_mapFilenameHERleak
The filename of the magnetic field map.
virtual ROOT::Math::XYZVector calculate(const ROOT::Math::XYZVector &point) const override
Calculates the magnetic field vector at the specified space point.
std::vector< std::vector< ParamPoint3 >::const_iterator > m_offset_pp_her
The vector of pointer to accelerate search in maps for her.
ranges_t m_ranges_her
ranges vector for HER
std::vector< ParamPoint3 > m_h3
The map for HER.
double getAperture(double s, std::vector< ApertPoint >::const_iterator hint) const
Returns the beam pipe aperture at given position.
std::vector< std::vector< ApertPoint >::const_iterator > m_offset_ap_her
The vector of pointer to accelerate search in aperture for her.
std::vector< ParamPoint3 > m_l3
The map for LER.
std::string m_mapFilenameHER
Magnetic field map of HER
double getApertureHER(double s) const
Returns the HER beam pipe aperture at given position.
std::vector< std::vector< ParamPoint3 >::const_iterator > m_offset_pp_ler
The vector of pointer to accelerate search in maps for ler.
std::vector< range_t > ranges_t
vector of Range data structure.
ranges_t m_ranges_ler
ranges vector for LER
std::string m_apertFilenameHER
Filename of the aperture for HER.
std::string m_apertFilenameLER
The filename of the aperture for LER.
static const double speedOfLight
[cm/ns]
Definition: Const.h:695
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
static bool fileExists(const std::string &filename)
Check if the file with given filename exists.
Definition: FileSystem.cc:32
static const double m
[meters]
Definition: Unit.h:69
static const double s
[second]
Definition: Unit.h:95
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.
Quadrupole lense data structure.