2 Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
5 ///////////////////////////////////////////////////////////////////
6 // CurvilinearParametersT.icc, (c) ATLAS Detector software
7 ///////////////////////////////////////////////////////////////////
10 #include "GaudiKernel/MsgStream.h"
12 #include "TrkEventPrimitives/ParamDefs.h"
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)
25 this->m_position = Amg::Vector3D(parameters[x], parameters[y], parameters[z]);
26 this->m_momentum = Amg::Vector3D(parameters[3], parameters[4], parameters[5]);
28 // flip the charge according to qOverP
29 if (parameters[6] < 0.) {
30 this->m_chargeDef.setCharge(-1);
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();
40 /* we need all the above to be there for the surfac*/
41 this->m_surface = S(this->m_position, curvilinearFrame());
44 // Constructor with TP arguments
45 template<int DIM, class T, class S>
46 CurvilinearParametersT<DIM, T, S>::CurvilinearParametersT(
47 const Amg::Vector3D& pos,
51 std::optional<AmgSymMatrix(DIM)> cov,
52 unsigned int cIdentifier)
53 : ParametersBase<DIM, T>(std::move(cov))
54 , m_cIdentifier(cIdentifier)
56 this->m_position = pos;
57 // flip the charge according to qOverP
59 this->m_chargeDef.setCharge(-1.);
61 this->m_chargeDef.setCharge(1.);
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;
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));
77 /* we need all the above for the surface*/
78 this->m_surface = S(this->m_position, curvilinearFrame());
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,
87 std::optional<AmgSymMatrix(DIM)> cov,
88 unsigned int cIdentifier)
89 : ParametersBase<DIM, T>(std::move(cov))
90 , m_cIdentifier(cIdentifier)
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();
100 charge = 1.; // such that below is 1./mom.mag()
103 this->m_parameters[Trk::qOverP] = charge / mom.mag();
104 this->m_position = pos;
105 this->m_momentum = mom;
107 // we need all the above to create the surface
108 m_surface = S(this->m_position, curvilinearFrame());
111 /** the curvilinear parameters identifier */
112 template<int DIM, class T, class S>
114 CurvilinearParametersT<DIM, T, S>::cIdentifier() const
116 return m_cIdentifier;
119 template<int DIM, class T, class S>
121 CurvilinearParametersT<DIM, T, S>::setcIdentifier(unsigned int cIdentifier)
123 m_cIdentifier = cIdentifier;
126 /** Test to see if there's a surface there. */
127 template<int DIM, class T, class S>
129 CurvilinearParametersT<DIM, T, S>::hasSurface() const
134 /** Access to the Surface method */
135 template<int DIM, class T, class S>
137 CurvilinearParametersT<DIM, T, S>::associatedSurface() const
143 template<int DIM, class T, class S>
145 CurvilinearParametersT<DIM, T, S>::operator==(
146 const ParametersBase<DIM, T>& rhs) const
148 // tolerance for comparing matrices
149 constexpr double tolerance = 1e-8;
151 // make sure we compare objects of same type
152 if (this->type() != rhs.type()){
155 const auto pCasted = static_cast<decltype(this)>(&rhs);
156 // comparison to myself?
157 if (pCasted == this) {
161 // compare identifier
162 if (cIdentifier() != pCasted->cIdentifier()) {
167 CurvilinearUVT local_curvilinearFrame = curvilinearFrame();
168 CurvilinearUVT casted_curvilinearFrame = pCasted->curvilinearFrame();
169 if (!local_curvilinearFrame.curvU().isApprox(casted_curvilinearFrame.curvU(),
173 if (!local_curvilinearFrame.curvV().isApprox(casted_curvilinearFrame.curvV(),
177 if (!local_curvilinearFrame.curvT().isApprox(casted_curvilinearFrame.curvT(),
182 if (associatedSurface() != pCasted->associatedSurface()) {
186 // return compatibility of base class parts
187 return ParametersBase<DIM, T>::operator==(rhs);
191 template<int DIM, class T, class S>
192 CurvilinearParametersT<DIM, T, S>*
193 CurvilinearParametersT<DIM, T, S>::clone() const
195 return new CurvilinearParametersT<DIM, T, S>(*this);
198 /** Return the ParametersType enum */
199 template<int DIM, class T, class S>
200 inline constexpr ParametersType
201 CurvilinearParametersT<DIM, T, S>::type() const
203 return Trk::Curvilinear;
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
211 return S::staticType;
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
219 Amg::RotationMatrix3D mFrame;
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();
226 // return the rotation matrix that defines the curvilinear parameters
230 template<int DIM, class T, class S>
232 CurvilinearParametersT<DIM, T, S>::curvilinearFrame() const
234 CurvilinearUVT curvilinFrame(this->momentum().unit());
235 return curvilinFrame;
238 // Screen output dump
239 template<int DIM, class T, class S>
241 CurvilinearParametersT<DIM, T, S>::dump(MsgStream& out) const
243 out << "CurvilinearParametersT parameters:" << std::endl;
244 ParametersBase<DIM, T>::dump(out);
249 // Screen output dump
250 template<int DIM, class T, class S>
252 CurvilinearParametersT<DIM, T, S>::dump(std::ostream& out) const
254 out << "CurvilinearParametersT parameters:" << std::endl;
255 ParametersBase<DIM, T>::dump(out);
260 // private updateParametersHelper
261 template<int DIM, class T, class S>
263 CurvilinearParametersT<DIM, T, S>::updateParametersHelper(
264 const AmgVector(DIM) & updatedParameters)
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]);
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));
287 // position update if needed - loc1
288 if (updatedParameters[Trk::loc1] != 0.) {
290 updatedParameters[Trk::loc1] * curvilinearFrame().curvU();
292 // position update if needed - loc2
293 if (updatedParameters[Trk::loc2] != 0.) {
295 updatedParameters[Trk::loc2] * curvilinearFrame().curvV();
297 // Reset also the surface
298 m_surface = S(this->m_position, curvilinearFrame());
301 } // end of namespace Trk