ATLAS Offline Software
Loading...
Searching...
No Matches
ParametersT.icc
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3*/
4
5///////////////////////////////////////////////////////////////////
6// ParametersT.icc, (c) ATLAS Detector software
7///////////////////////////////////////////////////////////////////
8
9// Trk
10#include "TrkEventPrimitives/ParamDefs.h"
11#include "CxxUtils/trapping_fp.h"
12#include "CxxUtils/inline_hints.h"
13#include <cmath>
14// STD
15#include <limits>
16#include <utility>
17
18namespace Trk {
19
20// Constructor with local arguments - uses global <-> local for parameters
21template<int DIM, class T, class S>
22ParametersT<DIM, T, S>::ParametersT(double loc1,
23 double loc2,
24 double phi,
25 double theta,
26 double qop,
27 const S& surface,
28 std::optional<AmgSymMatrix(DIM)> cov)
29 : ParametersBase<DIM, T>(AmgVector(DIM)::Zero(), std::move(cov), sgn(qop))
30 , SurfaceUniqHolderImpl<S>(surface)
31{
32 // Tell clang to optimize assuming that FP exceptions can trap.
33 // Otherwise, it can vectorize the division, which can lead to
34 // spurious division-by-zero traps from unused vector lanes.
35 CXXUTILS_TRAPPING_FP;
36 this->m_position = Amg::Vector3D::Zero();
37 this->m_momentum = Amg::Vector3D::Zero();
38
39 // check qoverp is physical
40 double p = 0.;
41 if (qop != 0) {
42 p = std::abs(1. / qop);
43 } else {
44 // qop is unphysical. No momentum measurement.
45 p = InvalidParam::INVALID_P;
46 qop = InvalidParam::INVALID_QOP;
47 }
48
49 // fill the parameters
50 // cppcheck-suppress constStatement
51 this->m_parameters << loc1, loc2, phi, theta, qop;
52
53 // now calculate the momentum
54 double sinPhi, cosPhi, cosTheta, sinTheta;
55 sincos(phi, &sinPhi, &cosPhi);
56 sincos(theta, &sinTheta, &cosTheta);
57 this->m_momentum =
58 Amg::Vector3D(p * cosPhi * sinTheta, p * sinPhi * sinTheta, p * cosTheta);
59 this->m_associatedSurface->localToGlobal(this->localPosition(),
60 this->m_momentum, this->m_position);
61}
62
63// Constructor with local arguments - uses global <-> local for parameters
64template<int DIM, class T, class S>
65ParametersT<DIM, T, S>::ParametersT(const AmgVector(DIM) & parameters,
66 const S& surface,
67 std::optional<AmgSymMatrix(DIM)> cov)
68 : ParametersBase<DIM, T>(parameters,
69 std::move(cov),
70 sgn(parameters[Trk::qOverP]))
71 , SurfaceUniqHolderImpl<S>(surface)
72{
73 this->m_position = Amg::Vector3D::Zero();
74 this->m_momentum = Amg::Vector3D::Zero();
75 // decide the sign of the charge
76 double qop = this->m_parameters[Trk::qOverP];
77
78 // check qoverp is physical
79 double p = 0.;
80 if (qop != 0.) {
81 p = std::abs(1. / qop);
82 } else {
83 // qop is unphysical. No momentum measurement.
84 p = InvalidParam::INVALID_P;
85 qop = InvalidParam::INVALID_QOP;
86 }
87
88 // fill momentum & then position using the surface
89 this->m_momentum =
90 Amg::Vector3D(p * std::cos(this->m_parameters[Trk::phi]) *
91 std::sin(this->m_parameters[Trk::theta]),
92 p * std::sin(this->m_parameters[Trk::phi]) *
93 std::sin(this->m_parameters[Trk::theta]),
94 p * std::cos(this->m_parameters[Trk::theta]));
95
96 this->m_associatedSurface->localToGlobal(this->localPosition(),
97 this->m_momentum, this->m_position);
98}
99
100// Constructor with global arguments - uses global <-> local for parameters */
101template<int DIM, class T, class S>
102ParametersT<DIM, T, S>::ParametersT(const Amg::Vector3D& pos,
103 const Amg::Vector3D& mom,
104 double charge,
105 const S& surface,
106 std::optional<AmgSymMatrix(DIM)> cov)
107 : ParametersBase<DIM, T>(AmgVector(DIM)::Zero(), std::move(cov), charge)
108 , SurfaceUniqHolderImpl<S>(surface)
109{
110
111 this->m_position = pos;
112 this->m_momentum = mom;
113
114 // get the local parameters via the surface
115 Amg::Vector2D lPosition;
116 const bool ok = this->m_associatedSurface->globalToLocal(
117 this->position(), this->momentum(), lPosition);
118 if (not ok) {
119 lPosition = Amg::Vector2D(InvalidParam::INVALID, InvalidParam::INVALID);
120 }
121
122 // For a neutral particle, last parm should be 1/p rather than q/p.
123 double qopnum = this->charge();
124 if (qopnum == 0) {
125 qopnum = 1;
126 }
127
128 // fill the vector now
129 this->m_parameters << lPosition[Trk::loc1], lPosition[Trk::loc2],
130 this->momentum().phi(), this->momentum().theta(),
131 qopnum / this->momentum().norm();
132}
133
134// Constructor with mixed arguments 1 - uses global <-> local for parameters
135template<int DIM, class T, class S>
136Trk::ParametersT<DIM, T, S>::ParametersT(const Amg::Vector3D& pos,
137 double phi,
138 double theta,
139 double qop,
140 const S& surface,
141 std::optional<AmgSymMatrix(DIM)> cov)
142 : ParametersBase<DIM, T>(AmgVector(DIM)::Zero(), std::move(cov), 1.)
143 , SurfaceUniqHolderImpl<S>(surface)
144{
145
146 this->m_position = pos;
147 this->m_momentum = Amg::Vector3D::Zero();
148
149 // decide the sign of the charge
150 if (qop < 0.) {
151 this->m_chargeDef.setCharge(-1);
152 }
153
154 // fill momentum & then position using the surface
155 double p = 0.0;
156 if (qop != 0.) {
157 p = std::abs(1. / qop);
158 } else {
159 // qop is unphysical. No momentum measurement.
160 p = InvalidParam::INVALID_P;
161 qop = InvalidParam::INVALID_QOP;
162 }
163 double sinPhi, cosPhi, cosTheta, sinTheta;
164 sincos(phi, &sinPhi, &cosPhi);
165 sincos(theta, &sinTheta, &cosTheta);
166 this->m_momentum =
167 Amg::Vector3D(p * cosPhi * sinTheta, p * sinPhi * sinTheta, p * cosTheta);
168
169 // get the local parameters via the surface
170 Amg::Vector2D lPosition;
171 const bool ok = this->m_associatedSurface->globalToLocal(
172 this->position(), this->momentum(), lPosition);
173 if (not ok) {
174 lPosition = Amg::Vector2D(InvalidParam::INVALID, InvalidParam::INVALID);
175 }
176 // fill the vector now
177 // cppcheck-suppress constStatement
178 this->m_parameters << lPosition[Trk::loc1], lPosition[Trk::loc2], phi, theta, qop;
179}
180
181/** Test to see if there's a surface there. */
182template<int DIM, class T, class S>
183inline bool
184ParametersT<DIM, T, S>::hasSurface() const
185{
186 return (this->m_associatedSurface != nullptr);
187}
188
189/** Access to the Surface method */
190template<int DIM, class T, class S>
191inline const S&
192ParametersT<DIM, T, S>::associatedSurface() const
193{
194 if (not this->m_associatedSurface){
195 throw std::runtime_error(
196 "No associated surface in ParametersT::associatedSurface");
197 }
198 return *(this->m_associatedSurface);
199}
200
201/** equality operator */
202template<int DIM, class T, class S>
203bool
204ParametersT<DIM, T, S>::operator==(const ParametersBase<DIM, T>& rhs) const
205{
206 // make sure we compare objects of same type
207 decltype(this) pCasted = dynamic_cast<decltype(this)>(&rhs);
208 if (!pCasted) {
209 return false;
210 }
211
212 // comparison to myself?
213 if (pCasted == this) {
214 return true;
215 }
216
217 // compare surfaces
218 if (associatedSurface() != pCasted->associatedSurface()) {
219 return false;
220 }
221
222 // return compatibility of base class parts
223 return ParametersBase<DIM, T>::operator==(rhs);
224}
225template<int DIM, class T, class S>
226bool
227ParametersT<DIM, T, S>::operator==(const ParametersT& rhs) const
228{
229 return *this == static_cast<const ParametersBase<DIM, T>&>(rhs);
230}
231/** clone */
232template<int DIM, class T, class S>
233ParametersT<DIM, T, S>*
234ParametersT<DIM, T, S>::clone() const
235{
236 return new ParametersT<DIM, T, S>(*this);
237}
238
239/** Return the ParametersType enum */
240template<int DIM, class T, class S>
241inline constexpr ParametersType
242ParametersT<DIM, T, S>::type() const
243{
244 return Trk::AtaSurface;
245}
246
247/** Return the Surface Type (check SurfaceType enums)*/
248template<int DIM, class T, class S>
249inline constexpr SurfaceType
250ParametersT<DIM, T, S>::surfaceType() const
251{
252 return S::staticType;
253}
254
255// return the measurementFrame
256template<int DIM, class T, class S>
257Amg::RotationMatrix3D
258ParametersT<DIM, T, S>::measurementFrame() const
259{
260 return associatedSurface().measurementFrame(this->position(),
261 this->momentum());
262}
263
264// private updateParametersHelper
265template<int DIM, class T, class S>
266ATH_FLATTEN
267void
268ParametersT<DIM, T, S>::updateParametersHelper(const AmgVector(DIM) &
269 updatedParameters)
270{
271 // valid to use != here, because value is either copied or modified,
272 bool updatePosition =
273 (updatedParameters[Trk::loc1] != this->m_parameters[Trk::loc1]) ||
274 (updatedParameters[Trk::loc2] != this->m_parameters[Trk::loc2]);
275
276 bool updateMomentum =
277 (updatedParameters[Trk::phi] != this->m_parameters[Trk::phi]) ||
278 (updatedParameters[Trk::theta] != this->m_parameters[Trk::theta]) ||
279 (updatedParameters[qOverP] != this->m_parameters[qOverP]);
280
281 if constexpr (S::staticType == Trk::SurfaceType::Line ||
282 S::staticType == Trk::SurfaceType::Perigee) {
283 updatePosition |= updateMomentum;
284 }
285
286 // update the parameters vector
287 this->m_parameters = updatedParameters;
288
289 // momentum update is needed
290 if (updateMomentum) {
291 double phi = this->m_parameters[Trk::phi];
292 double theta = this->m_parameters[Trk::theta];
293 double p = InvalidParam::INVALID_P;
294 if (this->m_parameters[Trk::qOverP] != 0) {
295 this->m_chargeDef = sgn(this->m_parameters[Trk::qOverP]);
296 p = std::abs(1. / this->m_parameters[Trk::qOverP]);
297 }
298 double sinPhi, cosPhi, cosTheta, sinTheta;
299 sincos(phi, &sinPhi, &cosPhi);
300 sincos(theta, &sinTheta, &cosTheta);
301 this->m_momentum = Amg::Vector3D(p * cosPhi * sinTheta,
302 p * sinPhi * sinTheta, p * cosTheta);
303 }
304
305 // position or momentum update needed
306 if (updatePosition) {
307 if (this->m_associatedSurface) {
308 this->m_associatedSurface->localToGlobal(
309 this->localPosition(), this->m_momentum, this->m_position);
310 } else {
311 this->m_position.setZero();
312 }
313 }
314}
315
316// Protected Constructor with local arguments - persistency only
317template<int DIM, class T, class S>
318ParametersT<DIM, T, S>::ParametersT(const AmgVector(DIM) & pars,
319 const S* surface,
320 std::optional<AmgSymMatrix(DIM)> cov)
321 : ParametersBase<DIM, T>(pars, std::move(cov))
322 , SurfaceUniqHolderImpl<S>(surface)
323{
324 float qop = this->m_parameters[Trk::qOverP];
325 // decide the sign of the charge
326 if (qop < 0.) {
327 this->m_chargeDef.setCharge(-1);
328 }
329 double p = 0.0;
330 if (qop != 0.) {
331 p = std::abs(1. / qop);
332 } else {
333 // qop is unphysical. No momentum measurement.
334 p = InvalidParam::INVALID_P;
335 qop = InvalidParam::INVALID_QOP;
336 }
337 if (this->m_associatedSurface) {
338 // fill momentum & then position using the surface
339 this->m_momentum = Amg::Vector3D(
340 p * std::cos(this->m_parameters[Trk::phi]) * std::sin(this->m_parameters[Trk::theta]),
341 p * std::sin(this->m_parameters[Trk::phi]) * std::sin(this->m_parameters[Trk::theta]),
342 p * std::cos(this->m_parameters[Trk::theta]));
343 this->m_associatedSurface->localToGlobal(
344 this->localPosition(), this->m_momentum, this->m_position);
345 } else {
346 this->m_momentum.setZero();
347 this->m_position.setZero();
348 }
349}
350
351// Screen output dump
352template<int DIM, class T, class S>
353MsgStream&
354ParametersT<DIM, T, S>::dump(MsgStream& out) const
355{
356 out << "ParametersT parameters " << std::endl;
357 ParametersBase<DIM, T>::dump(out);
358
359 return out;
360}
361
362// Screen output dump
363template<int DIM, class T, class S>
364std::ostream&
365ParametersT<DIM, T, S>::dump(std::ostream& out) const
366{
367 out << "ParametersT parameters:" << std::endl;
368 ParametersBase<DIM, T>::dump(out);
369
370 return out;
371}
372
373} // end of namespace Trk
Generated on for ATLAS Offline Software by doxygen 1.14.0