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