Belle II Software development
RiemannsMethod Class Reference

Class implementing the Riemann fit for two dimensional trajectory circle. More...

#include <RiemannsMethod.h>

Public Member Functions

 RiemannsMethod ()
 Default constructor.
 
void update (TrackingUtilities::CDCTrajectory2D &trajectory2D, CDCObservations2D &observations2D) const
 Executes the fit and updates the trajectory parameters. This may render the information in the observation object.
 
bool isLineConstrained () const
 Getter for the indicator that lines should be fitted by this fitter.
 
bool isOriginConstrained () const
 Getter for the indicator that curves through the origin should be fitted by this fitter.
 
void setLineConstrained (bool constrained=true)
 Indicator if this fitter is setup to fit lines.
 
void setOriginConstrained (bool constrained=true)
 Indicator if this fitter is setup to fit curves through the origin.
 

Private Member Functions

void updateWithoutDriftLength (TrackingUtilities::CDCTrajectory2D &trajectory2D, CDCObservations2D &observations2D) const
 Executes the fit without using the drift length information.
 
void updateWithDriftLength (TrackingUtilities::CDCTrajectory2D &trajectory2D, CDCObservations2D &observations2D) const
 Executes the fit using the drift length information.
 

Private Attributes

bool m_lineConstrained
 Memory for the flag indicating that lines should be fitted.
 
bool m_originConstrained
 Memory for the flag indicating that curves through the origin shall be fitter.
 

Detailed Description

Class implementing the Riemann fit for two dimensional trajectory circle.

Definition at line 22 of file RiemannsMethod.h.

Constructor & Destructor Documentation

◆ RiemannsMethod()

Default constructor.

Definition at line 28 of file RiemannsMethod.cc.

29 : m_lineConstrained(false)
30 , m_originConstrained(false)
31{
32}
bool m_originConstrained
Memory for the flag indicating that curves through the origin shall be fitter.
bool m_lineConstrained
Memory for the flag indicating that lines should be fitted.

Member Function Documentation

◆ isLineConstrained()

bool isLineConstrained ( ) const
inline

Getter for the indicator that lines should be fitted by this fitter.

Definition at line 50 of file RiemannsMethod.h.

51 {
52 return m_lineConstrained;
53 }

◆ isOriginConstrained()

bool isOriginConstrained ( ) const
inline

Getter for the indicator that curves through the origin should be fitted by this fitter.

Definition at line 56 of file RiemannsMethod.h.

57 {
58 return m_originConstrained;
59 }

◆ setLineConstrained()

void setLineConstrained ( bool constrained = true)
inline

Indicator if this fitter is setup to fit lines.

Definition at line 62 of file RiemannsMethod.h.

63 {
64 m_lineConstrained = constrained;
65 }

◆ setOriginConstrained()

void setOriginConstrained ( bool constrained = true)
inline

Indicator if this fitter is setup to fit curves through the origin.

Definition at line 68 of file RiemannsMethod.h.

69 {
70 m_originConstrained = constrained;
71 }

◆ update()

void update ( TrackingUtilities::CDCTrajectory2D & trajectory2D,
CDCObservations2D & observations2D ) const

Executes the fit and updates the trajectory parameters. This may render the information in the observation object.

Definition at line 34 of file RiemannsMethod.cc.

36{
37
38 if (observations2D.getNObservationsWithDriftRadius() > 0) {
39 updateWithDriftLength(trajectory2D, observations2D);
40 } else {
41 updateWithoutDriftLength(trajectory2D, observations2D);
42 }
43
44 EForwardBackward isCoaligned = observations2D.isCoaligned(trajectory2D);
45 if (isCoaligned == EForwardBackward::c_Backward) trajectory2D.reverse();
46}
TrackingUtilities::EForwardBackward isCoaligned(const TrackingUtilities::CDCTrajectory2D &trajectory2D) const
Checks if the last observation in the vector lies at greater or lower travel distance than the last o...
std::size_t getNObservationsWithDriftRadius() const
Returns the number of observations having a drift radius radius.
void updateWithoutDriftLength(TrackingUtilities::CDCTrajectory2D &trajectory2D, CDCObservations2D &observations2D) const
Executes the fit without using the drift length information.
void updateWithDriftLength(TrackingUtilities::CDCTrajectory2D &trajectory2D, CDCObservations2D &observations2D) const
Executes the fit using the drift length information.
EForwardBackward
Enumeration to represent the distinct possibilities of the right left passage information.

◆ updateWithDriftLength()

void updateWithDriftLength ( TrackingUtilities::CDCTrajectory2D & trajectory2D,
CDCObservations2D & observations2D ) const
private

Executes the fit using the drift length information.

This method is used if there is no drift length information is available from the observations.

Definition at line 147 of file RiemannsMethod.cc.

148{
149 EigenObservationMatrix eigenObservation = getEigenObservationMatrix(&observations2D);
150 assert(eigenObservation.cols() == 4);
151 size_t nObservations = observations2D.size();
152
153 //cout << "updateWithRightLeft : " << endl;
154 //observations always have the structure
155 /*
156 observables[0][0] == x of first point
157 observables[0][1] == y of first point
158 observables[0][2] == desired distance of first point
159 */
160
161 Eigen::Matrix< double, Eigen::Dynamic, 1 > distances = eigenObservation.col(2);
162
163 //cout << "distances : " << endl << distances << endl;
164
165 if ((isLineConstrained()) && (isOriginConstrained())) {
166
167
168 Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > projectedPoints(nObservations, 2);
169
170 //all coordinates
171 //projectedPoints.col(0) = Eigen::Matrix<double,Eigen::Dynamic,1>::Constant(nObservations,1.0);
172 projectedPoints.col(0) = eigenObservation.col(0);
173 projectedPoints.col(1) = eigenObservation.col(1);
174 //projectedPoints.col(2) = eigenObservation.leftCols<2>().rowwise().squaredNorm() - distances.rowwise().squaredNorm();
175
176 Eigen::Matrix<double, 2, 1> parameters =
177 projectedPoints.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(distances);
178
179 trajectory2D.setGlobalCircle(UncertainPerigeeCircle(PerigeeCircle::fromN(0.0, parameters(0), parameters(1), 0.0)));
180 }
181
182 else if ((! isLineConstrained()) && (isOriginConstrained())) {
183
184 Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > projectedPoints(nObservations, 3);
185
186 //all coordinates
187 //projectedPoints.col(0) = Eigen::Matrix<double,Eigen::Dynamic,1>::Constant(1.0);
188 projectedPoints.col(0) = eigenObservation.col(0);
189 projectedPoints.col(1) = eigenObservation.col(1);
190 projectedPoints.col(2) = eigenObservation.leftCols<2>().rowwise().squaredNorm() - distances.rowwise().squaredNorm();
191
192 Eigen::Matrix<double, 3, 1> parameters =
193 projectedPoints.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(distances);
194
195 trajectory2D.setGlobalCircle(UncertainPerigeeCircle(PerigeeCircle::fromN(0.0, parameters(0), parameters(1), parameters(2))));
196 }
197
198 else if ((isLineConstrained()) && (! isOriginConstrained())) {
199
200 Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > projectedPoints(nObservations, 3);
201
202 //all coordinates
203 projectedPoints.col(0) = Eigen::Matrix<double, Eigen::Dynamic, 1>::Constant(nObservations, 1.0);
204 projectedPoints.col(1) = eigenObservation.col(0);
205 projectedPoints.col(2) = eigenObservation.col(1);
206 //projectedPoints.col(3) = eigenObservation.leftCols<2>().rowwise().squaredNorm()- distances.rowwise().squaredNorm();
207
208 Eigen::Matrix<double, 3, 1> parameters =
209 projectedPoints.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(distances);
210
211 trajectory2D.setGlobalCircle(UncertainPerigeeCircle(PerigeeCircle::fromN(parameters(0), parameters(1), parameters(2), 0.0)));
212 //fit.setParameters(parameters(0),parameters(1),parameters(2),0.0);
213
214 }
215
216 else if ((! isLineConstrained()) && (! isOriginConstrained())) {
217
218 Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > projectedPoints(nObservations, 4);
219
220 //all coordinates
221 projectedPoints.col(0) = Eigen::Matrix<double, Eigen::Dynamic, 1>::Constant(nObservations, 1.0);
222 projectedPoints.col(1) = eigenObservation.col(0);
223 projectedPoints.col(2) = eigenObservation.col(1);
224 projectedPoints.col(3) = eigenObservation.leftCols<2>().rowwise().squaredNorm() - distances.rowwise().squaredNorm();
225
226 //cout << "points : " << endl << projectedPoints << endl;
227
228 Eigen::Matrix<double, 4, 1> parameters =
229 projectedPoints.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(distances);
230
231 trajectory2D.setGlobalCircle(UncertainPerigeeCircle(PerigeeCircle::fromN(parameters(0), parameters(1), parameters(2),
232 parameters(3))));
233
234 }
235}
std::size_t size() const
Returns the number of observations stored.
bool isLineConstrained() const
Getter for the indicator that lines should be fitted by this fitter.
bool isOriginConstrained() const
Getter for the indicator that curves through the origin should be fitted by this fitter.
static PerigeeCircle fromN(double n0, double n1, double n2, double n3=0)
Constructor with the four parameters of the generalized circle.

◆ updateWithoutDriftLength()

void updateWithoutDriftLength ( TrackingUtilities::CDCTrajectory2D & trajectory2D,
CDCObservations2D & observations2D ) const
private

Executes the fit without using the drift length information.

This method is used if there is drift length information is available from the observations.

Definition at line 50 of file RiemannsMethod.cc.

52{
53 EigenObservationMatrix eigenObservation = getEigenObservationMatrix(&observations2D);
54 assert(eigenObservation.cols() == 4);
55 size_t nObservations = observations2D.size();
56
57 if (isLineConstrained()) {
58
59 //do a normal line fit
60 Eigen::Matrix<double, Eigen::Dynamic, 2> points = eigenObservation.leftCols<2>();
61
62 Eigen::Matrix<double, 1, 2> pointMean;
63 //RowVector2d pointMean;
64 // cppcheck-suppress constStatement
65 pointMean << 0.0, 0.0;
66 if (!(isOriginConstrained())) {
67 // subtract the offset from the origin
68 pointMean = points.colwise().mean();
69
70 points = points.rowwise() - pointMean;
71 }
72 Eigen::Matrix<double, 2, 2> covarianceMatrix = points.transpose() * points;
73
74 Eigen::SelfAdjointEigenSolver< Eigen::Matrix<double, 2, 2> > eigensolver(covarianceMatrix);
75
76 if (eigensolver.info() != Eigen::Success)
77 B2WARNING("SelfAdjointEigenSolver could not compute the eigen values of the observation matrix");
78
79 //the eigenvalues are generated in increasing order
80 //we are interested in the lowest one since we want to compute the normal vector of the plane
81
82 Eigen::Matrix<double, 2, 1> normalToLine = eigensolver.eigenvectors().col(0);
83
84 double offset = -pointMean * normalToLine;
85
86 // set the generalized circle parameters
87 // last set to zero constrains to a line
88 trajectory2D.setGlobalCircle(UncertainPerigeeCircle(PerigeeCircle::fromN(offset, normalToLine(0), normalToLine(1), 0)));
89
90 } else {
91
92 //lift the points to the projection space
93 Eigen::Matrix<double, Eigen::Dynamic, 3> projectedPoints(nObservations, 3);
94
95 projectedPoints.col(0) = eigenObservation.col(0);
96 projectedPoints.col(1) = eigenObservation.col(1);
97 projectedPoints.col(2) = eigenObservation.leftCols<2>().rowwise().squaredNorm();
98
99
100 Eigen::Matrix<double, 1, 3> pointMean;
101 // cppcheck-suppress constStatement
102 pointMean << 0.0, 0.0, 0.0;
103 if (!(isOriginConstrained())) {
104 // subtract the offset from the origin
105 pointMean = projectedPoints.colwise().mean();
106
107 projectedPoints = projectedPoints.rowwise() - pointMean;
108 }
109
110 Eigen::Matrix<double, 3, 3> covarianceMatrix = projectedPoints.transpose() * projectedPoints;
111
112 Eigen::SelfAdjointEigenSolver< Eigen::Matrix<double, 3, 3> > eigensolver(covarianceMatrix);
113
114 if (eigensolver.info() != Eigen::Success)
115 B2WARNING("Eigen::SelfAdjointEigenSolver could not compute the eigen values of the observation matrix");
116
117 //the eigenvalues are generated in increasing order
118 //we are interested in the lowest one since we want to compute the normal vector of the plane
119
120 Eigen::Matrix<double, 3, 1> normalToPlane = eigensolver.eigenvectors().col(0);
121
122 double offset = -pointMean * normalToPlane;
123
124 trajectory2D.setGlobalCircle(UncertainPerigeeCircle(PerigeeCircle::fromN(offset, normalToPlane(0), normalToPlane(1),
125 normalToPlane(2))));
126 //fit.setParameters();
127
128 }
129
130 //check if the orientation is alright
131 const ROOT::Math::XYVector& directionAtCenter = trajectory2D.getFlightDirection2D(ROOT::Math::XYVector(0.0, 0.0));
132
133
134 size_t voteForChangeSign = 0;
135 for (size_t iPoint = 0; iPoint < nObservations; ++iPoint) {
136 double pointInSameDirection = eigenObservation(iPoint, 0) * directionAtCenter.x() +
137 eigenObservation(iPoint, 1) * directionAtCenter.y();
138 if (pointInSameDirection < 0) ++voteForChangeSign;
139 }
140
141 if (voteForChangeSign > nObservations / 2.0) trajectory2D.reverse();
142
143}

Member Data Documentation

◆ m_lineConstrained

bool m_lineConstrained
private

Memory for the flag indicating that lines should be fitted.

Definition at line 75 of file RiemannsMethod.h.

◆ m_originConstrained

bool m_originConstrained
private

Memory for the flag indicating that curves through the origin shall be fitter.

Definition at line 78 of file RiemannsMethod.h.


The documentation for this class was generated from the following files: