ATLAS Offline Software
Loading...
Searching...
No Matches
SolenoidParametrization.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2017, 2019 CERN for the benefit of the ATLAS collaboration
3*/
4
5/***************************************************************************
6 Fast (approximate) methods for solenoidal field properties
7 ----------------------------------------------------------
8 ***************************************************************************/
9
10//<<<<<< INCLUDES >>>>>>
11
12#include <algorithm>
13#include <iomanip>
14#include "EventPrimitives/EventPrimitives.h"
15#include "GaudiKernel/SystemOfUnits.h"
16#include "TrkExSolenoidalIntersector/SolenoidParametrization.h"
17#include "GaudiKernel/MsgStream.h"
18
19#include "MagFieldElements/AtlasFieldCache.h"
20
21namespace {
22 class RestoreIOSFlags
23 {
24 public:
25 explicit RestoreIOSFlags (std::ostream &os)
26 : m_os(&os),
27 m_precision(m_os->precision())
28 {}
29 ~RestoreIOSFlags() {
30 m_os->precision(m_precision);
31 }
32 private:
33 std::ostream *m_os;
34 std::streamsize m_precision;
35 };
36}
37
38
39namespace Trk
40{
41
42//<<<<<< PRIVATE VARIABLE DEFINITIONS >>>>>>
43
44const double SolenoidParametrization::s_binInvSizeTheta = 1./0.1;
45const double SolenoidParametrization::s_binInvSizeZ = 1./20.*Gaudi::Units::mm;
46const double SolenoidParametrization::s_binZeroTheta = 0.;
47const double SolenoidParametrization::s_binZeroZ = -160.*Gaudi::Units::mm;
48const double SolenoidParametrization::s_lightSpeed = -1.*299792458*Gaudi::Units::m/Gaudi::Units::s;
49const int SolenoidParametrization::s_maxBinTheta = 72;
50const int SolenoidParametrization::s_maxBinZ = 17;
51const double SolenoidParametrization::s_maximumImpactAtOrigin= 30.*Gaudi::Units::mm;
52const double SolenoidParametrization::s_maximumZatOrigin = 250.*Gaudi::Units::mm;
53const int SolenoidParametrization::s_numberParameters = 6;
54const double SolenoidParametrization::s_rInner = 570.*Gaudi::Units::mm;
55const double SolenoidParametrization::s_rOuter = 1050.*Gaudi::Units::mm;
56const double SolenoidParametrization::s_zInner = 2150.0*Gaudi::Units::mm; // just after wheel #7
57const double SolenoidParametrization::s_zOuter = 2800.0*Gaudi::Units::mm; // just after wheel #9
58
59
60SolenoidParametrization::BinParameters::BinParameters (const double r,
61 const double z,
62 const double cotTheta)
63{
64 if (cotTheta > 0) {
65 m_signTheta = 1;
66 m_cotTheta = cotTheta;
67 m_zAtAxis = z - r*cotTheta;
68 }
69 else {
70 m_signTheta = -1;
71 m_cotTheta = -cotTheta;
72 m_zAtAxis = r*cotTheta - z;
73 }
74}
75
76
77SolenoidParametrization::Parameters::Parameters (const SolenoidParametrization& spar,
78 const double r,
79 const double z,
80 const double cotTheta)
81 : BinParameters (r, z, cotTheta)
82{
83 int key = fieldKey(*this);
84 if (r > s_rInner || m_signTheta*z > s_zInner)
85 {
86 key += s_numberParameters/2;
87 }
88 spar.setTerms (key, *this);
89}
90
91
92
93 SolenoidParametrization::SolenoidParametrization(const AtlasFieldCacheCondObj &field_cond_obj)
94 : m_fieldCondObj (&field_cond_obj),
95 m_parameters ()
96{
97
98 MagField::AtlasFieldCache fieldCache;
99 m_fieldCondObj->getInitializedCache (fieldCache);
100 m_centralField = fieldComponent(0.,0.,0., fieldCache);
101 // now parametrise field - if requested
102 {
103 parametrizeSolenoid();
104 }
105}
106
107//<<<<<< PRIVATE MEMBER FUNCTION DEFINITIONS >>>>>>
108
109void
110SolenoidParametrization::parametrizeSolenoid(void){
111 // set parametrisation granularity (up to cotTheta = 7.)
112 // get value of cubic term for approx: Bz = Bcentral*(1 - term * z^3)
113 // 'fit' to average over cotTheta lines
114 double smallOffset = 0.0000000000001; // avoid FPE
115 double zAtAxis = s_binZeroZ; // + smallOffset ?
116 MagField::AtlasFieldCache fieldCache;
117 m_fieldCondObj->getInitializedCache (fieldCache);
118
119 constexpr int n = 200;
120 Amg::VectorX difference(n); // it is filled below in the loop over n
121 Amg::MatrixX derivative(n,s_numberParameters); //set zero below
122 for (int binZ = 0; binZ < s_maxBinZ; ++binZ){
123 double cotTheta = smallOffset;
124 for (int binTheta = 0; binTheta < s_maxBinTheta - 1; ++binTheta) {
125 double r = 0.;
126 double z = zAtAxis;
127 double dr;
128 if (cotTheta < s_zOuter/s_rOuter){
129 dr = s_rOuter/double(n);
130 } else {
131 dr = s_zOuter/(cotTheta*double(n));
132 }
133 derivative.setZero();
134 for (int k = 0; k < n; ++k){
135 r += dr;
136 z += dr*cotTheta;
137 double w = (n - k)*(n - k);
138 double zLocal = z - zAtAxis;
139 if (r > s_rInner || z > s_zInner){
140 derivative(k,3) = w;
141 derivative(k,4) = w*zLocal*zLocal;
142 derivative(k,5) = w*zLocal*zLocal*zLocal;
143 } else {
144 derivative(k,0) = w;
145 derivative(k,1) = w*zLocal*zLocal;
146 derivative(k,2) = w*zLocal*zLocal*zLocal;
147 }
148 difference(k) = w*(fieldComponent(r,z,cotTheta, fieldCache) - m_centralField);
149 }
150 // solve for parametrization coefficients
151 Amg::VectorX solution = derivative.colPivHouseholderQr().solve(difference);
152 BinParameters parms (zAtAxis, cotTheta);
153 int key = fieldKey(parms);
154 assert (m_parameters[key] == 0.);
155 m_parameters[key++] = m_centralField + solution(0);
156 m_parameters[key++] = solution(1);
157 m_parameters[key++] = solution(2);
158 m_parameters[key++] = m_centralField + solution(3);
159 m_parameters[key++] = solution(4);
160 m_parameters[key++] = solution(5);
161 // duplicate last z-bin for contiguous neighbour lookup
162 if (binZ == s_maxBinZ - 1){
163 assert (m_parameters[key] == 0.);
164 m_parameters[key++] = m_centralField + solution(0);
165 m_parameters[key++] = solution(1);
166 m_parameters[key++] = solution(2);
167 m_parameters[key++] = m_centralField + solution(3);
168 m_parameters[key++] = solution(4);
169 m_parameters[key++] = solution(5);
170 key -= s_numberParameters;
171 }
172
173 // duplicate next to previous z-bin for contiguous neighbour lookup
174 if (binZ > 0){
175 key -= 2*s_numberParameters*s_maxBinTheta;
176 assert (m_parameters[key] == 0.);
177 m_parameters[key++] = m_centralField + solution(0);
178 m_parameters[key++] = solution(1);
179 m_parameters[key++] = solution(2);
180 m_parameters[key++] = m_centralField + solution(3);
181 m_parameters[key++] = solution(4);
182 m_parameters[key++] = solution(5);
183 }
184 cotTheta += 1./s_binInvSizeTheta;
185 } // Loop over binTheta
186 zAtAxis += 1./s_binInvSizeZ;
187 }
188 //
189 // duplicate end theta-bins for contiguous neighbour lookup
190 zAtAxis = s_binZeroZ; // + smallOffset ??
191 for (int binZ = 0; binZ < s_maxBinZ; ++binZ){
192 double cotTheta = double(s_maxBinTheta)/s_binInvSizeTheta;
193 BinParameters parms (zAtAxis, cotTheta);
194 int key = fieldKey(parms);
195 for (int k = 0; k < 2*s_numberParameters; ++k){
196 assert (m_parameters[key+2*s_numberParameters] == 0.);
197 m_parameters[key+2*s_numberParameters] = m_parameters[key];
198 ++key;
199 }
200 zAtAxis += 1./s_binInvSizeZ;
201 }
202}
203
204//<<<<<< PUBLIC MEMBER FUNCTION DEFINITIONS >>>>>>
205
206void
207SolenoidParametrization::printFieldIntegrals (MsgStream& msg) const
208{
209 // integrate along lines of const eta from origin to r = 1m or |z| = 2.65m
210 // direction normalised st transverse component = 1 (equiv to a fixed pt)
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"
215 << " asymm: x y z"
216 << " " << std::endl
217 << " m T T/m T T/m "
218 << " T T m T.m T.m"
219 << " T T T"
220 << "/n";
221
222 double maxR = 1000.*Gaudi::Units::mm;
223 double maxZ = 2650.*Gaudi::Units::mm;
224 int numSteps = 1000;
225
226 MagField::AtlasFieldCache fieldCache;
227 m_fieldCondObj->getInitializedCache (fieldCache);
228
229 // step through eta-range
230 double eta = 0.;
231 for (int i = 0; i != 31; ++i)
232 {
233 double phi = 0.;
234 double cotTheta = std::sinh(eta);
235 Amg::Vector3D direction(std::cos(phi),std::sin(phi),cotTheta);
236 double rEnd = maxR;
237 if (std::abs(cotTheta) > maxZ/maxR) rEnd = maxZ/direction.z();
238 double step = rEnd/static_cast<double>(numSteps); // radial step in mm
239 Amg::Vector3D position(0.,0.,0.);
240 double meanBL = 0.;
241 double meanBT = 0.;
242 double meanBZ = 0.;
243 double maxBT = 0.;
244 double minBT = 9999.;
245 double maxGradBT = 0.;
246 double maxGradBZ = 0.;
247 double prevBT = 0.;
248 double prevBZ = 0.;
249 double reverseZ = 0.;
250
251 // look up field along eta-line
252 for (int j = 0; j != numSteps; ++j)
253 {
254 position += 0.5*step*direction;
255 Amg::Vector3D field;
256 fieldCache.getField(position.data(),field.data());
257 Amg::Vector3D vCrossB = direction.cross(field);
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);
262 meanBL += BL;
263 meanBT += BT;
264 meanBZ += BZ;
265 if (BT > maxBT) maxBT = BT;
266 if (BT < minBT) minBT = BT;
267 if (j > 0)
268 {
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;
274 }
275 prevBT = BT;
276 prevBZ = BZ;
277 }
278
279 // normalize
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;
287
288 msg << std::setw(6) << std::setprecision(2) << eta
289 << std::setw(8) << std::setprecision(3) << rEnd /Gaudi::Units::meter
290 << std::setw(11) << std::setprecision(4) << meanBZ /Gaudi::Units::tesla
291 << std::setw(12) << std::setprecision(3) << maxGradBZ /Gaudi::Units::tesla
292 << std::setw(13) << std::setprecision(4) << meanBT /Gaudi::Units::tesla
293 << std::setw(12) << std::setprecision(3) << maxGradBT /Gaudi::Units::tesla
294 << std::setw(13) << std::setprecision(4) << minBT /Gaudi::Units::tesla
295 << std::setw(8) << std::setprecision(4) << maxBT /Gaudi::Units::tesla;
296 if (reverseZ > 0.)
297 {
298 msg << std::setw(17) << std::setprecision(2) << reverseZ /Gaudi::Units::meter
299 << std::setw(18) << std::setprecision(4) << integralBT /(Gaudi::Units::tesla*Gaudi::Units::meter)
300 << std::setw(8) << std::setprecision(4) << integralBL /(Gaudi::Units::tesla*Gaudi::Units::meter)
301 << " ";
302 }
303 else
304 {
305 msg << std::setw(35) << std::setprecision(4) << integralBT /(Gaudi::Units::tesla*Gaudi::Units::meter)
306 << std::setw(8) << std::setprecision(4) << integralBL /(Gaudi::Units::tesla*Gaudi::Units::meter)
307 << " ";
308 }
309
310 // check symmetry (reflect in each axis)
311 for (int k = 0; k != 3; ++k)
312 {
313 if (k == 0) direction = Amg::Vector3D(-std::cos(phi),std::sin(phi),cotTheta);
314 if (k == 1) direction = Amg::Vector3D(std::cos(phi),-std::sin(phi),cotTheta);
315 if (k == 2) direction = Amg::Vector3D(std::cos(phi),std::sin(phi),-cotTheta);
316 position = Amg::Vector3D(0.,0.,0.);
317 double asymm = 0.;
318 // look up field along eta-line
319 for (int j = 0; j != numSteps; ++j)
320 {
321 position += 0.5*step*direction;
322 Amg::Vector3D field;
323 fieldCache.getField(position.data(),field.data());
324 Amg::Vector3D vCrossB = direction.cross(field);
325 position += 0.5*step*direction;
326 double BT = vCrossB.x()*direction.y() - vCrossB.y()*direction.x();
327 asymm += BT;
328 }
329 asymm = asymm/static_cast<double>(numSteps) - meanBT;
330 msg << std::setw(9) << std::setprecision(4) << asymm /Gaudi::Units::tesla;
331 }
332 msg<<"/n";
333 eta += 0.1;
334 }
335}
336
337void
338SolenoidParametrization::printParametersForEtaLine (double eta, double z_origin, MsgStream & msg) const
339{
340 double cotTheta = 1./std::tan(2.*std::atan(1./std::exp(eta)));
341 BinParameters parms (z_origin, cotTheta);
342 int key = fieldKey(parms);
343 double z_max;
344 if (cotTheta < s_zInner/s_rInner)
345 {
346 z_max = s_rInner*cotTheta;
347 }
348 else
349 {
350 z_max = s_zInner;
351 }
352 msg <<__func__<<"\n"
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)
357 << m_parameters[key]/m_centralField
358 << " z^2 "<< std::setw(6) << std::setprecision(3)
359 << m_parameters[key+1]*z_max*z_max/m_centralField
360 << " z^3 " << std::setw(6) << std::setprecision(3)
361 << m_parameters[key+2]*z_max*z_max*z_max/m_centralField
362 << " outer terms: z0 "<< std::setw(6) << std::setprecision(3)
363 << m_parameters[key+3]/m_centralField
364 << " z^2 "<< std::setw(6) << std::setprecision(3)
365 << m_parameters[key+4]*z_max*z_max/m_centralField
366 << " z^3 " << std::setw(6) << std::setprecision(3)
367 << m_parameters[key+5]*z_max*z_max*z_max/m_centralField
368 << std::resetiosflags(std::ios::fixed) << "\n";
369}
370
371void
372SolenoidParametrization::printResidualForEtaLine (double eta, double zOrigin, MsgStream & msg) const
373{
374 double cotTheta = 1./std::tan(2.*std::atan(1./std::exp(std::abs(eta))));
375 double z = zOrigin;
376 double r = 0.;
377 int n = 200;
378 double dr;
379 if (cotTheta < s_zOuter/s_rOuter)
380 {
381 dr = s_rOuter/double(n);
382 }
383 else
384 {
385 dr = s_zOuter/(cotTheta*double(n));
386 }
387 double chiSquareIn = 0.;
388 double chiSquareOut = 0.;
389 double nIn = 0.;
390 double nOut = 0.;
391 double worstBCalc = 0.;
392 double worstBTrue = 0.;
393 double worstDiff = -1.;
394 double worstR = 0.;
395 double worstZ = 0.;
396 MagField::AtlasFieldCache fieldCache;
397 m_fieldCondObj->getInitializedCache (fieldCache);
398 for (int k = 0; k < n; ++k)
399 {
400 double b = fieldComponent(r,z,cotTheta,fieldCache);
401 Parameters parms (*this, r, z, cotTheta);
402 double diff = (fieldComponent(z, parms) - b)/s_lightSpeed;
403
404 if (r > s_rInner || z > s_zInner)
405 {
406 chiSquareOut+= diff*diff;
407 nOut += 1.;
408 }
409 else
410 {
411 chiSquareIn += diff*diff;
412 nIn += 1.;
413 }
414
415 if (std::abs(diff) > worstDiff)
416 {
417 worstDiff = std::abs(diff);
418 worstBCalc = fieldComponent(z, parms);
419 worstBTrue = b;
420 worstR = r;
421 worstZ = z;
422 }
423 r += dr;
424 z += dr*cotTheta;
425 }
426
427 msg <<__func__<<"\n"
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)
433 << std::sqrt(chiSquareOut/nOut) /Gaudi::Units::tesla << std::endl
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)
438 << worstBTrue/s_lightSpeed /Gaudi::Units::tesla
439 << " " << std::setw(6) << std::setprecision(3) << worstBCalc/s_lightSpeed /Gaudi::Units::tesla
440 << std::resetiosflags(std::ios::fixed) << "\n";
441}
442
443
444
445
446} // end of namespace
447
448
eta
Scalar eta() const
pseudorapidity method
Definition AmgMatrixBasePlugin.h:83
diff
void diff(const Jet &rJet1, const Jet &rJet2, std::map< std::string, double > varDiff)
Difference between jets - Non-Class function required by trigger.
Definition Jet.cxx:631
MagField::AtlasFieldCache
Local cache for magnetic field (based on MagFieldServices/AtlasFieldSvcTLS.h)
Definition AtlasFieldCache.h:43
MagField::AtlasFieldCache::getField
void getField(const double *ATH_RESTRICT xyz, double *ATH_RESTRICT bxyz, double *ATH_RESTRICT deriv=nullptr)
get B field value at given position xyz[3] is in mm, bxyz[3] is in kT if deriv[9] is given,...
Definition AtlasFieldCache.cxx:42
Trk::SolenoidParametrization::BinParameters::BinParameters
BinParameters(const double zAtAxis, const double cotTheta)
Definition SolenoidParametrization.h:37
Trk::SolenoidParametrization::Parameters::Parameters
Parameters(const SolenoidParametrization &spar, const double r, const double z, const double cotTheta)
Definition SolenoidParametrization.cxx:77
Trk::SolenoidParametrization::fieldKey
static int fieldKey(BinParameters &parms)
Definition SolenoidParametrization.h:131
Trk::SolenoidParametrization::printParametersForEtaLine
void printParametersForEtaLine(double eta, double z_origin, MsgStream &msg) const
Definition SolenoidParametrization.cxx:338
Trk::SolenoidParametrization::setTerms
void setTerms(int, Parameters &parms) const
Definition SolenoidParametrization.h:196
Trk::SolenoidParametrization::SolenoidParametrization
SolenoidParametrization(const AtlasFieldCacheCondObj &field_cond_obj)
Definition SolenoidParametrization.cxx:93
Trk::SolenoidParametrization::m_fieldCondObj
const AtlasFieldCacheCondObj * m_fieldCondObj
Definition SolenoidParametrization.h:117
Trk::SolenoidParametrization::printFieldIntegrals
void printFieldIntegrals(MsgStream &m) const
Definition SolenoidParametrization.cxx:207
Trk::SolenoidParametrization::fieldComponent
double fieldComponent(double z, const Parameters &parms) const
Definition SolenoidParametrization.h:225
Trk::SolenoidParametrization::printResidualForEtaLine
void printResidualForEtaLine(double eta, double zOrigin, MsgStream &msg) const
Definition SolenoidParametrization.cxx:372
r
int r
Definition globals.cxx:22
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition EventPrimitives.h:27
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
Amg::VectorX
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.
Definition EventPrimitives.h:30
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition FakeTrackBuilder.h:9
Trk::z
@ z
global position (cartesian)
Definition ParamDefs.h:57
Trk::phi
@ phi
Definition ParamDefs.h:75
Trk::binZ
@ binZ
Definition BinningType.h:49
msg
MsgStream & msg
Definition testRead.cxx:32
Generated on for ATLAS Offline Software by doxygen 1.14.0