Belle II Software light-2406-ragdoll
UncertainHelix.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8#include <framework/dataobjects/UncertainHelix.h>
9
10#include <TMatrixD.h>
11
12#include <cassert>
13
14using namespace Belle2;
15using namespace HelixParameterIndex;
16
18 Helix(),
19 m_covariance(5),
20 m_pValue(0)
21{
22}
23
24UncertainHelix::UncertainHelix(const ROOT::Math::XYZVector& position,
25 const ROOT::Math::XYZVector& momentum,
26 const short int charge,
27 const double bZ,
28 const TMatrixDSym& cartesianCovariance,
29 const double pValue) :
30 Helix(ROOT::Math::XYZVector(0.0, 0.0, position.Z()), momentum, charge, bZ),
31 m_covariance(cartesianCovariance), // Initialize the covariance matrix to the 6x6 covariance and reduce it afterwards
32 m_pValue(pValue)
33{
34 // Maybe push these out of this function:
35 // Indices of the cartesian coordinates
36 const int iX = 0;
37 const int iY = 1;
38 const int iZ = 2;
39 const int iPx = 3;
40 const int iPy = 4;
41 const int iPz = 5;
42
43 // We initialised the m_covariance to the cartesian covariance and
44 // reduce it now to the real 5x5 matrix that should be there.
45
46 // 1. Rotate to a system where phi0 = 0
47 TMatrixD jacobianRot(6, 6);
48 jacobianRot.Zero();
49
50 const double px = momentum.X();
51 const double py = momentum.Y();
52 const double pt = hypot(px, py);
53 const double cosPhi0 = px / pt;
54 const double sinPhi0 = py / pt;
55
56 // Passive rotation matrix by phi0:
57 jacobianRot(iX, iX) = cosPhi0;
58 jacobianRot(iX, iY) = sinPhi0;
59 jacobianRot(iY, iX) = -sinPhi0;
60 jacobianRot(iY, iY) = cosPhi0;
61 jacobianRot(iZ, iZ) = 1.0;
62
63 jacobianRot(iPx, iPx) = cosPhi0;
64 jacobianRot(iPx, iPy) = sinPhi0;
65 jacobianRot(iPy, iPx) = -sinPhi0;
66 jacobianRot(iPy, iPy) = cosPhi0;
67 jacobianRot(iPz, iPz) = 1.0;
68
69 m_covariance.Similarity(jacobianRot);
70
71 // 2. Translate to perigee parameters on the position
72 const double pz = momentum.Z();
73 const double invPt = 1 / pt;
74 const double invPtSquared = invPt * invPt;
75 const double alpha = getAlpha(bZ);
76
77 TMatrixD jacobianToHelixParameters(5, 6);
78 jacobianToHelixParameters.Zero();
79
80 jacobianToHelixParameters(iD0, iY) = -1;
81 jacobianToHelixParameters(iPhi0, iX) = charge * invPt / alpha;
82 jacobianToHelixParameters(iPhi0, iPy) = invPt;
83 jacobianToHelixParameters(iOmega, iPx) = -charge * invPtSquared / alpha;
84 jacobianToHelixParameters(iTanLambda, iPx) = - pz * invPtSquared;
85 jacobianToHelixParameters(iTanLambda, iPz) = invPt;
86 jacobianToHelixParameters(iZ0, iX) = - pz * invPt;
87 jacobianToHelixParameters(iZ0, iZ) = 1;
88 m_covariance.Similarity(jacobianToHelixParameters);
89
90 // The covariance m_covariance is now the correct 5x5 covariance matrix.
91 assert(m_covariance.GetNrows() == 5);
92
93 // 3. Extrapolate to the origin.
94 /* const double arcLength2D = */ passiveMoveBy(-position.X(), -position.Y(), 0.0);
95}
96
97
99 const double& phi0,
100 const double& omega,
101 const double& z0,
102 const double& tanLambda,
103 const TMatrixDSym& covariance,
104 const double pValue) :
105 Helix(d0, phi0, omega, z0, tanLambda),
106 m_covariance(covariance),
107 m_pValue(pValue)
108{
109}
110
111
112
113TMatrixDSym UncertainHelix::getCartesianCovariance(const double bZ_tesla) const
114{
115 // 0. Define indices
116 // Maybe push these out of this function:
117 // Indices of the cartesian coordinates
118 const int iX = 0;
119 const int iY = 1;
120 const int iZ = 2;
121 const int iPx = 3;
122 const int iPy = 4;
123 const int iPz = 5;
124
125 // Transform covariance matrix
126 TMatrixD jacobianInflate(6, 5);
127 jacobianInflate.Zero();
128
129 // 1. Inflate the perigee covariance to a cartesian covariance where phi0 == 0 is assumed
130 // The real phi0 is a simple rotation which can be handled in the next step.
131 // Jacobian matrix for the translation
132
133 const double& d0 = getD0();
134 const double& omega = getOmega();
135 const double& tanLambda = getTanLambda();
136
137 const double alpha = getAlpha(bZ_tesla);
138 const double absAlphaOmega = alpha * std::fabs(omega);
139 const double signedAlphaOmega2 = absAlphaOmega * omega;
140
141 const double invAbsAlphaOmega = 1.0 / absAlphaOmega;
142 const double invSignedAlphaOmega2 = 1.0 / signedAlphaOmega2;
143
144 // Position after the move.
145 jacobianInflate(iX, iPhi0) = d0;
146 jacobianInflate(iY, iD0) = -1.0;
147 jacobianInflate(iZ, iZ0) = 1.0;
148
149 // Momentum
150 jacobianInflate(iPx, iOmega) = -invSignedAlphaOmega2;
151 jacobianInflate(iPy, iPhi0) = invAbsAlphaOmega;
152 jacobianInflate(iPz, iOmega) = -tanLambda * invSignedAlphaOmega2;
153 jacobianInflate(iPz, iTanLambda) = invAbsAlphaOmega;
154
155 TMatrixDSym cov6 = m_covariance; //copy
156 cov6.Similarity(jacobianInflate);
157
159 const double& cosPhi0 = getCosPhi0();
160 const double& sinPhi0 = getSinPhi0();
161
162 TMatrixD jacobianRot(6, 6);
163 jacobianRot.Zero();
164
165 // Active rotation matrix by phi0:
166 jacobianRot(iX, iX) = cosPhi0;
167 jacobianRot(iX, iY) = -sinPhi0;
168 jacobianRot(iY, iX) = sinPhi0;
169 jacobianRot(iY, iY) = cosPhi0;
170 jacobianRot(iZ, iZ) = 1.0;
171
172 jacobianRot(iPx, iPx) = cosPhi0;
173 jacobianRot(iPx, iPy) = -sinPhi0;
174 jacobianRot(iPy, iPx) = sinPhi0;
175 jacobianRot(iPy, iPy) = cosPhi0;
176 jacobianRot(iPz, iPz) = 1.0;
177
178 cov6.Similarity(jacobianRot);
179 return cov6;
180}
181
183{
185
186 // D0, omega and tan lambda have to be taken to their opposites
187 // Phi0 is augmented by pi which does not change its covariances
188 // Z0 stays the same
189 TMatrixD jacobianReverse(5, 5);
190 jacobianReverse.UnitMatrix();
191 jacobianReverse(iD0, iD0) = -1;
192 jacobianReverse(iOmega, iOmega) = -1;
193 jacobianReverse(iTanLambda, iTanLambda) = -1;
194
195 m_covariance.Similarity(jacobianReverse);
196}
197
198double UncertainHelix::passiveMoveBy(const double& byX,
199 const double& byY,
200 const double& byZ)
201{
202 // Move the covariance matrix first to have access to the original parameters
203 TMatrixD jacobianPassiveMove(5, 5);
204 calcPassiveMoveByJacobian(byX, byY, jacobianPassiveMove);
205 m_covariance.Similarity(jacobianPassiveMove);
206 return Helix::passiveMoveBy(byX, byY, byZ);
207}
This class represents an ideal helix in perigee parameterization.
Definition: Helix.h:59
double getSinPhi0() const
Getter for the cosine of the azimuth angle of travel direction at the perigee.
Definition: Helix.h:384
void reverse()
Reverses the direction of travel of the helix in place.
Definition: Helix.cc:119
double getCosPhi0() const
Getter for the cosine of the azimuth angle of travel direction at the perigee.
Definition: Helix.h:381
static double getAlpha(const double bZ)
Calculates the alpha value for a given magnetic field in Tesla.
Definition: Helix.cc:109
double getOmega() const
Getter for omega, which is a signed curvature measure of the track.
Definition: Helix.h:387
double getD0() const
Getter for d0, which is the signed distance to the perigee in the r-phi plane.
Definition: Helix.h:372
double passiveMoveBy(const ROOT::Math::XYZVector &by)
Moves origin of the coordinate system (passive transformation) by the given vector.
Definition: Helix.h:245
double getTanLambda() const
Getter for tan lambda, which is the z over two dimensional arc length slope of the track.
Definition: Helix.h:393
TMatrixD calcPassiveMoveByJacobian(const ROOT::Math::XYZVector &by, const double expandBelowChi=M_PI/8) const
Calculate the 5x5 jacobian matrix for the transport of the helix parameters, when moving the origin o...
Definition: Helix.cc:301
UncertainHelix()
Default constuctor initialising all members to zero.
void reverse()
Reverses the direction of travel of the helix in place.
TMatrixDSym getCartesianCovariance(const double bZ_tesla=1.5) const
Getter for the position and momentum covariance matrix.
double passiveMoveBy(const ROOT::Math::XYZVector &by)
Moves origin of the coordinate system (passive transformation) by the given vector.
TMatrixDSym m_covariance
5x5 covariance of the perigee parameters.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24