Initializes the magnetic field component.
This method opens the magnetic field map file.
Magnetic field data structure.
fold rotation to/from lense frame to a single matrix multiplicaton and vector addition
In case several maps in the same position we can simply sum up all matrices since magnetic field has superposition properties as well as keep only one vector of parameters for each beamline
44{
45
47 B2FATAL("The filename for the HER quadrupole magnetic field component is empty !");
48 return;
49 }
51 B2FATAL("The filename for the LER quadrupole magnetic field component is empty !");
52 return;
53 }
55 B2FATAL("The filename for the HER beam pipe aperture is empty !");
56 return;
57 }
59 B2FATAL("The filename for the LER beam pipe aperture is empty !");
60 return;
61 }
62
63
66 B2FATAL(
"The HER quadrupole magnetic field map file '" <<
m_mapFilenameHER <<
"' could not be found !");
67 return;
68 }
71 B2FATAL(
"The LER quadrupole magnetic field map file '" <<
m_mapFilenameLER <<
"' could not be found !");
72 return;
73 }
76 B2FATAL(
"The HERleak quadrupole magnetic field map file '" <<
m_mapFilenameHERleak <<
"' could not be found !");
77 return;
78 }
82 return;
83 }
87 return;
88 }
89
90
91
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
106
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
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
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
141
142 auto readAperturePoints = [this](istream & IN) -> vector<ApertPoint> {
143 vector<ApertPoint> r;
145 while (IN >> t.s >> t.r)
146 {
147
148
149 t.s /= 10; t.r /= 10;
150 r.push_back(t);
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
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 {
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
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
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;
241 if (p.s == t.s && p.L == t.L)
242 p += t;
243 else
244 v.insert(i, t);
245 }
246 };
248
255 auto getranges = [](const vector<ParamPoint3>& v) {
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++) {
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
272
273 const double inf = numeric_limits<double>::infinity();
274
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 };
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 };
316}
std::vector< ApertPoint > m_al
The the aperture parameters for LER.
std::vector< ApertPoint > m_ah
The the aperture parameters for HER.
std::string m_mapFilenameLER
Magnetic field map of LER.
std::string m_mapFilenameHERleak
The filename of the magnetic field map.
std::vector< ParamPoint3 > m_h3
The map for HER.
std::vector< ParamPoint3 > m_l3
The map for LER.
std::string m_mapFilenameHER
Magnetic field map of HER.
std::vector< range_t > ranges_t
vector of Range data structure.
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]
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...
static bool fileExists(const std::string &filename)
Check if the file with given filename exists.
static const double m
[meters]
static const double s
[second]
std::vector< std::vector< double > > merge(std::vector< std::vector< std::vector< double > > > toMerge)
merge { vector<double> a, vector<double> b} into {a, b}
Quadrupole lense data structure.