ATLAS Offline Software
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 
18 namespace Trk {
19 
20 // Constructor with local arguments - uses global <-> local for parameters
21 template<int DIM, class T, class S>
22 ParametersT<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
64 template<int DIM, class T, class S>
65 ParametersT<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 */
101 template<int DIM, class T, class S>
102 ParametersT<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
135 template<int DIM, class T, class S>
136 Trk::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. */
182 template<int DIM, class T, class S>
183 inline bool
184 ParametersT<DIM, T, S>::hasSurface() const
185 {
186  return (this->m_associatedSurface != nullptr);
187 }
188 
189 /** Access to the Surface method */
190 template<int DIM, class T, class S>
191 inline const S&
192 ParametersT<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 */
202 template<int DIM, class T, class S>
203 bool
204 ParametersT<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 }
225 template<int DIM, class T, class S>
226 bool
227 ParametersT<DIM, T, S>::operator==(const ParametersT& rhs) const
228 {
229  return *this == static_cast<const ParametersBase<DIM, T>&>(rhs);
230 }
231 /** clone */
232 template<int DIM, class T, class S>
233 ParametersT<DIM, T, S>*
234 ParametersT<DIM, T, S>::clone() const
235 {
236  return new ParametersT<DIM, T, S>(*this);
237 }
238 
239 /** Return the ParametersType enum */
240 template<int DIM, class T, class S>
241 inline constexpr ParametersType
242 ParametersT<DIM, T, S>::type() const
243 {
244  return Trk::AtaSurface;
245 }
246 
247 /** Return the Surface Type (check SurfaceType enums)*/
248 template<int DIM, class T, class S>
249 inline constexpr SurfaceType
250 ParametersT<DIM, T, S>::surfaceType() const
251 {
252  return S::staticType;
253 }
254 
255 // return the measurementFrame
256 template<int DIM, class T, class S>
257 Amg::RotationMatrix3D
258 ParametersT<DIM, T, S>::measurementFrame() const
259 {
260  return associatedSurface().measurementFrame(this->position(),
261  this->momentum());
262 }
263 
264 // private updateParametersHelper
265 template<int DIM, class T, class S>
266 ATH_FLATTEN
267 void
268 ParametersT<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
317 template<int DIM, class T, class S>
318 ParametersT<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
352 template<int DIM, class T, class S>
353 MsgStream&
354 ParametersT<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
363 template<int DIM, class T, class S>
364 std::ostream&
365 ParametersT<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 Sun Jun 30 2024 21:23:54 for ATLAS Offline Software by   doxygen 1.8.18