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