47 B2FATAL(
"The filename for the HER quadrupole magnetic field component is empty !");
51 B2FATAL(
"The filename for the LER quadrupole magnetic field component is empty !");
55 B2FATAL(
"The filename for the HER beam pipe aperture is empty !");
59 B2FATAL(
"The filename for the LER beam pipe aperture is empty !");
66 B2FATAL(
"The HER quadrupole magnetic field map file '" <<
m_mapFilenameHER <<
"' could not be found !");
71 B2FATAL(
"The LER quadrupole magnetic field map file '" <<
m_mapFilenameLER <<
"' could not be found !");
76 B2FATAL(
"The HERleak quadrupole magnetic field map file '" <<
m_mapFilenameHERleak <<
"' could not be found !");
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));
116 auto readLenseParameters = [](istream & IN) -> vector<ParamPoint> {
117 vector<ParamPoint> r;
120 while (IN >> name >> t.s >> t.L >> t.K0 >> t.K1 >> t.SK0 >> t.SK1 >> t.ROTATE >> t.DX >> t.DY)
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);
142 auto readAperturePoints = [
this](istream & IN) -> vector<ApertPoint> {
143 vector<ApertPoint> r;
145 while (IN >> t.s >> t.r)
149 t.s /= 10; t.r /= 10;
153 r.push_back({numeric_limits<double>::infinity(),
m_maxr2});
157 io::filtering_istream apertFileHER, apertFileLER;
158 apertFileHER.push(io::file_source(fullPathApertHER));
159 apertFileLER.push(io::file_source(fullPathApertLER));
161 m_ah = readAperturePoints(apertFileHER);
162 m_al = readAperturePoints(apertFileLER);
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)
194#define ROTATE_DIRECTION 1
195#if ROTATE_DIRECTION==0
196 double sphi = 0, cphi = 1;
199 sincos(-b->ROTATE * (M_PI / 180), &sphi, &cphi);
200#if ROTATE_DIRECTION==-1
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;
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;
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;
221 const double p0_HER = 7.0e+9 / c, p0_LER = 4.0e+9 / c;
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);
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)
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});
266 r.push_back({s0, smax});
273 const double inf = numeric_limits<double>::infinity();
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;
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);
305 auto associate_lenses = [](
const vector<ParamPoint3>& ap,
const ranges_t& v) {
306 vector<std::vector<ParamPoint3>::const_iterator> res;
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;
349 const double sa = sin(0.0415), ca = cos(0.0415);
351 ROOT::Math::XYZVector B(0, 0, 0);
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)};
356 ROOT::Math::XYZVector vl{(xc + zs), -v.Y(), -(zc - xs)};
358 double r2h = vh.Perp2(), r2l = vl.Perp2();
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);
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);