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 26 of file RiemannsMethod.cc.

27 : m_lineConstrained(false)
28 , m_originConstrained(false)
29{
30}
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 32 of file RiemannsMethod.cc.

34{
35
36 if (observations2D.getNObservationsWithDriftRadius() > 0) {
37 updateWithDriftLength(trajectory2D, observations2D);
38 } else {
39 updateWithoutDriftLength(trajectory2D, observations2D);
40 }
41
42 EForwardBackward isCoaligned = observations2D.isCoaligned(trajectory2D);
43 if (isCoaligned == EForwardBackward::c_Backward) trajectory2D.reverse();
44}
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 145 of file RiemannsMethod.cc.

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

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

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: