ATLAS Offline Software
CurvilinearParametersT.icc
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 ///////////////////////////////////////////////////////////////////
6 // CurvilinearParametersT.icc, (c) ATLAS Detector software
7 ///////////////////////////////////////////////////////////////////
8 
9 // Gaudi
10 #include "GaudiKernel/MsgStream.h"
11 // Trk
12 #include "TrkEventPrimitives/ParamDefs.h"
13 #include <cmath>
14 
15 namespace Trk {
16 // Constructor with TP arguments
17 template<int DIM, class T, class S>
18 CurvilinearParametersT<DIM, T, S>::CurvilinearParametersT(
19  const AmgVector(DIM + 2) & parameters,
20  std::optional<AmgSymMatrix(DIM)> covariance,
21  unsigned int cIdentifier)
22  : ParametersBase<DIM, T>(std::move(covariance))
23  , m_cIdentifier(cIdentifier)
24 {
25  this->m_position = Amg::Vector3D(parameters[x], parameters[y], parameters[z]);
26  this->m_momentum = Amg::Vector3D(parameters[3], parameters[4], parameters[5]);
27 
28  // flip the charge according to qOverP
29  if (parameters[6] < 0.) {
30  this->m_chargeDef.setCharge(-1);
31  }
32  // assign the parameters
33  this->m_parameters[locX] = 0.;
34  this->m_parameters[locY] = 0.;
35  // get phi & theta from the momentum vector
36  this->m_parameters[phi] = this->momentum().phi();
37  this->m_parameters[theta] = this->momentum().theta();
38  this->m_parameters[qOverP] = parameters[6] / this->momentum().mag();
39 
40  /* we need all the above to be there for the surfac*/
41  this->m_surface = S(this->m_position, curvilinearFrame());
42 }
43 
44 // Constructor with TP arguments
45 template<int DIM, class T, class S>
46 CurvilinearParametersT<DIM, T, S>::CurvilinearParametersT(
47  const Amg::Vector3D& pos,
48  double tphi,
49  double ttheta,
50  double tqOverP,
51  std::optional<AmgSymMatrix(DIM)> cov,
52  unsigned int cIdentifier)
53  : ParametersBase<DIM, T>(std::move(cov))
54  , m_cIdentifier(cIdentifier)
55 {
56  this->m_position = pos;
57  // flip the charge according to qOverP
58  if (tqOverP < 0.) {
59  this->m_chargeDef.setCharge(-1.);
60  } else {
61  this->m_chargeDef.setCharge(1.);
62  }
63 
64  // assign the parameters
65  this->m_parameters[Trk::locX] = 0.;
66  this->m_parameters[Trk::locY] = 0.;
67  this->m_parameters[Trk::phi] = tphi;
68  this->m_parameters[Trk::theta] = ttheta;
69  this->m_parameters[Trk::qOverP] = tqOverP;
70 
71  // make sure that the position & momentum are calculated
72  double p = std::abs(1. / tqOverP);
73  this->m_momentum = Amg::Vector3D(p * std::cos(tphi) * std::sin(ttheta),
74  p * std::sin(tphi) * std::sin(ttheta),
75  p * std::cos(ttheta));
76 
77  /* we need all the above for the surface*/
78  this->m_surface = S(this->m_position, curvilinearFrame());
79 }
80 
81 // full global constructor
82 template<int DIM, class T, class S>
83 CurvilinearParametersT<DIM, T, S>::CurvilinearParametersT(
84  const Amg::Vector3D& pos,
85  const Amg::Vector3D& mom,
86  double charge,
87  std::optional<AmgSymMatrix(DIM)> cov,
88  unsigned int cIdentifier)
89  : ParametersBase<DIM, T>(std::move(cov))
90  , m_cIdentifier(cIdentifier)
91 {
92  this->m_chargeDef.setCharge(charge);
93  // assign the parameters
94  this->m_parameters[Trk::locX] = 0.;
95  this->m_parameters[Trk::locY] = 0.;
96  this->m_parameters[Trk::phi] = mom.phi();
97  this->m_parameters[Trk::theta] = mom.theta();
98 
99  if (charge == 0.) {
100  charge = 1.; // such that below is 1./mom.mag()
101  }
102 
103  this->m_parameters[Trk::qOverP] = charge / mom.mag();
104  this->m_position = pos;
105  this->m_momentum = mom;
106 
107  // we need all the above to create the surface
108  m_surface = S(this->m_position, curvilinearFrame());
109 }
110 
111 /** the curvilinear parameters identifier */
112 template<int DIM, class T, class S>
113 inline unsigned int
114 CurvilinearParametersT<DIM, T, S>::cIdentifier() const
115 {
116  return m_cIdentifier;
117 }
118 
119 template<int DIM, class T, class S>
120 inline void
121 CurvilinearParametersT<DIM, T, S>::setcIdentifier(unsigned int cIdentifier)
122 {
123  m_cIdentifier = cIdentifier;
124 }
125 
126 /** Test to see if there's a surface there. */
127 template<int DIM, class T, class S>
128 inline bool
129 CurvilinearParametersT<DIM, T, S>::hasSurface() const
130 {
131  return true;
132 }
133 
134 /** Access to the Surface method */
135 template<int DIM, class T, class S>
136 inline const S&
137 CurvilinearParametersT<DIM, T, S>::associatedSurface() const
138 {
139  return m_surface;
140 }
141 
142 // equality operator
143 template<int DIM, class T, class S>
144 bool
145 CurvilinearParametersT<DIM, T, S>::operator==(
146  const ParametersBase<DIM, T>& rhs) const
147 {
148  // tolerance for comparing matrices
149  constexpr double tolerance = 1e-8;
150 
151  // make sure we compare objects of same type
152  if (this->type() != rhs.type()){
153  return false;
154  }
155  const auto pCasted = static_cast<decltype(this)>(&rhs);
156  // comparison to myself?
157  if (pCasted == this) {
158  return true;
159  }
160 
161  // compare identifier
162  if (cIdentifier() != pCasted->cIdentifier()) {
163  return false;
164  }
165 
166  // compare UVT frame
167  CurvilinearUVT local_curvilinearFrame = curvilinearFrame();
168  CurvilinearUVT casted_curvilinearFrame = pCasted->curvilinearFrame();
169  if (!local_curvilinearFrame.curvU().isApprox(casted_curvilinearFrame.curvU(),
170  tolerance)) {
171  return false;
172  }
173  if (!local_curvilinearFrame.curvV().isApprox(casted_curvilinearFrame.curvV(),
174  tolerance)) {
175  return false;
176  }
177  if (!local_curvilinearFrame.curvT().isApprox(casted_curvilinearFrame.curvT(),
178  tolerance)) {
179  return false;
180  }
181  // compare surfaces
182  if (associatedSurface() != pCasted->associatedSurface()) {
183  return false;
184  }
185 
186  // return compatibility of base class parts
187  return ParametersBase<DIM, T>::operator==(rhs);
188 }
189 
190 /** clone */
191 template<int DIM, class T, class S>
192 CurvilinearParametersT<DIM, T, S>*
193 CurvilinearParametersT<DIM, T, S>::clone() const
194 {
195  return new CurvilinearParametersT<DIM, T, S>(*this);
196 }
197 
198 /** Return the ParametersType enum */
199 template<int DIM, class T, class S>
200 inline constexpr ParametersType
201 CurvilinearParametersT<DIM, T, S>::type() const
202 {
203  return Trk::Curvilinear;
204 }
205 
206 /** Return the Surface Type (check SurfaceType enums)*/
207 template<int DIM, class T, class S>
208 inline constexpr SurfaceType
209 CurvilinearParametersT<DIM, T, S>::surfaceType() const
210 {
211  return S::staticType;
212 }
213 
214 // Surface return (with on demand construction)
215 template<int DIM, class T, class S>
216 Amg::RotationMatrix3D
217 CurvilinearParametersT<DIM, T, S>::measurementFrame() const
218 {
219  Amg::RotationMatrix3D mFrame;
220  // the columnes
221  CurvilinearUVT local_curvilinearFrame = curvilinearFrame();
222  mFrame.col(0) = local_curvilinearFrame.curvU();
223  mFrame.col(1) = local_curvilinearFrame.curvV();
224  mFrame.col(2) = local_curvilinearFrame.curvT();
225 
226  // return the rotation matrix that defines the curvilinear parameters
227  return mFrame;
228 }
229 
230 template<int DIM, class T, class S>
231 CurvilinearUVT
232 CurvilinearParametersT<DIM, T, S>::curvilinearFrame() const
233 {
234  CurvilinearUVT curvilinFrame(this->momentum().unit());
235  return curvilinFrame;
236 }
237 
238 // Screen output dump
239 template<int DIM, class T, class S>
240 MsgStream&
241 CurvilinearParametersT<DIM, T, S>::dump(MsgStream& out) const
242 {
243  out << "CurvilinearParametersT parameters:" << std::endl;
244  ParametersBase<DIM, T>::dump(out);
245 
246  return out;
247 }
248 
249 // Screen output dump
250 template<int DIM, class T, class S>
251 std::ostream&
252 CurvilinearParametersT<DIM, T, S>::dump(std::ostream& out) const
253 {
254  out << "CurvilinearParametersT parameters:" << std::endl;
255  ParametersBase<DIM, T>::dump(out);
256 
257  return out;
258 }
259 
260 // private updateParametersHelper
261 template<int DIM, class T, class S>
262 void
263 CurvilinearParametersT<DIM, T, S>::updateParametersHelper(
264  const AmgVector(DIM) & updatedParameters)
265 {
266  // valid to use != here, because value is either copied or modified,
267  bool updateMomentum =
268  (updatedParameters[Trk::phi] != this->m_parameters[Trk::phi]) ||
269  (updatedParameters[Trk::theta] != this->m_parameters[Trk::theta]) ||
270  (updatedParameters[Trk::qOverP] != this->m_parameters[Trk::qOverP]);
271 
272  // momentum update is needed
273  if (updateMomentum) {
274  double phi = updatedParameters[Trk::phi];
275  double theta = updatedParameters[Trk::theta];
276  double p = std::abs(1. / updatedParameters[Trk::qOverP]);
277  this->m_chargeDef.setCharge(sgn(updatedParameters[Trk::qOverP]));
278  // assign them and update the momentum 3 vector
279  this->m_parameters[Trk::phi] = phi;
280  this->m_parameters[Trk::theta] = theta;
281  this->m_parameters[Trk::qOverP] = updatedParameters[Trk::qOverP];
282  this->m_momentum = Amg::Vector3D(p * std::cos(phi) * std::sin(theta),
283  p * std::sin(phi) * std::sin(theta),
284  p * std::cos(theta));
285  }
286 
287  // position update if needed - loc1
288  if (updatedParameters[Trk::loc1] != 0.) {
289  this->m_position +=
290  updatedParameters[Trk::loc1] * curvilinearFrame().curvU();
291  }
292  // position update if needed - loc2
293  if (updatedParameters[Trk::loc2] != 0.) {
294  this->m_position +=
295  updatedParameters[Trk::loc2] * curvilinearFrame().curvV();
296  }
297  // Reset also the surface
298  m_surface = S(this->m_position, curvilinearFrame());
299 }
300 
301 } // end of namespace Trk