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