15 #include "GaudiKernel/SystemOfUnits.h"
17 #include "GaudiKernel/MsgStream.h"
25 explicit RestoreIOSFlags (std::ostream &
os)
27 m_precision(m_os->precision())
30 m_os->precision(m_precision);
34 std::streamsize m_precision;
114 double smallOffset = 0.0000000000001;
119 constexpr
int n = 200;
124 for (
int binTheta = 0; binTheta <
s_maxBinTheta - 1; ++binTheta) {
134 for (
int k = 0;
k <
n; ++
k){
137 double w = (
n -
k)*(
n -
k);
138 double zLocal =
z - zAtAxis;
211 msg << __func__<<
"\n"
212 << std::setiosflags(std::ios::fixed)
213 <<
" eta rEnd mean(Bz) max(dBz/dR) mean(Bt) max(dBt/dR) "
214 <<
"min(Bt) max(Bt) reverse-bend(z) integrals: Bt.dR Bl.dR"
231 for (
int i = 0;
i != 31; ++
i)
237 if (std::abs(
cotTheta) > maxZ/maxR) rEnd = maxZ/direction.z();
238 double step = rEnd/
static_cast<double>(numSteps);
244 double minBT = 9999.;
245 double maxGradBT = 0.;
246 double maxGradBZ = 0.;
249 double reverseZ = 0.;
252 for (
int j = 0; j != numSteps; ++j)
254 position += 0.5*
step*direction;
258 position += 0.5*
step*direction;
259 double BZ =
field.z();
260 double BT = vCrossB.x()*direction.y() - vCrossB.y()*direction.x();
261 double BL = std::sqrt(vCrossB.mag2() - BT*BT);
265 if (BT > maxBT) maxBT = BT;
266 if (BT < minBT) minBT = BT;
269 if (BT*prevBT < 0.) reverseZ = position.z() -
step*direction.z();
270 double grad = std::abs(BT - prevBT);
271 if (grad > maxGradBT) maxGradBT = grad;
272 grad = std::abs(BZ - prevBZ);
273 if (grad > maxGradBZ) maxGradBZ = grad;
280 maxGradBT *=
static_cast<double>(numSteps)/
step;
281 maxGradBZ *=
static_cast<double>(numSteps)/
step;
282 meanBL /=
static_cast<double>(numSteps);
283 meanBT /=
static_cast<double>(numSteps);
284 meanBZ /=
static_cast<double>(numSteps);
285 double integralBL = meanBL*rEnd;
286 double integralBT = meanBT*rEnd;
288 msg << std::setw(6) << std::setprecision(2) << eta
311 for (
int k = 0;
k != 3; ++
k)
319 for (
int j = 0; j != numSteps; ++j)
321 position += 0.5*
step*direction;
325 position += 0.5*
step*direction;
326 double BT = vCrossB.x()*direction.y() - vCrossB.y()*direction.x();
329 asymm = asymm/
static_cast<double>(numSteps) - meanBT;
353 << std::setiosflags(std::ios::fixed)
354 <<
"SolenoidParametrization: line with eta " << std::setw(6) << std::setprecision(2) << eta
355 <<
" from (r,z) 0.0," << std::setw(6) << std::setprecision(1) << z_origin
356 <<
" inner terms: z0 "<< std::setw(6) << std::setprecision(2)
358 <<
" z^2 "<< std::setw(6) << std::setprecision(3)
360 <<
" z^3 " << std::setw(6) << std::setprecision(3)
362 <<
" outer terms: z0 "<< std::setw(6) << std::setprecision(3)
364 <<
" z^2 "<< std::setw(6) << std::setprecision(3)
366 <<
" z^3 " << std::setw(6) << std::setprecision(3)
368 << std::resetiosflags(std::ios::fixed) <<
"\n";
387 double chiSquareIn = 0.;
388 double chiSquareOut = 0.;
391 double worstBCalc = 0.;
392 double worstBTrue = 0.;
393 double worstDiff = -1.;
398 for (
int k = 0;
k <
n; ++
k)
415 if (std::abs(
diff) > worstDiff)
417 worstDiff = std::abs(
diff);
428 << std::setiosflags(std::ios::fixed)
429 <<
"SolenoidParametrization: line with eta " << std::setw(6) << std::setprecision(2) << eta
430 <<
" from (r,z) 0.0, " << std::setw(6) << std::setprecision(1) << zOrigin
431 <<
" rms diff inner/outer " << std::setw(6) << std::setprecision(3)
432 << std::sqrt(chiSquareIn/nIn) /
Gaudi::Units::tesla <<
" " << std::setw(6) << std::setprecision(3)
434 <<
" worst residual at: (r,z) "
435 << std::setw(6) << std::setprecision(1) << worstR
436 <<
", " << std::setw(6) << std::setprecision(1) << worstZ
437 <<
" with B true/calc " << std::setw(6) << std::setprecision(3)
440 << std::resetiosflags(std::ios::fixed) <<
"\n";