Belle II Software  release-05-01-25
RiemannsMethod.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2014 - 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 <tracking/trackFindingCDC/fitting/RiemannsMethod.h>
11 
12 #include <tracking/trackFindingCDC/fitting/EigenObservationMatrix.h>
13 #include <tracking/trackFindingCDC/fitting/CDCObservations2D.h>
14 
15 #include <tracking/trackFindingCDC/eventdata/trajectories/CDCTrajectory2D.h>
16 
17 #include <framework/logging/Logger.h>
18 
19 #include <Eigen/Eigen>
20 #include <Eigen/Core>
21 
22 #include <cassert>
23 
24 using namespace Belle2;
25 using namespace TrackFindingCDC;
26 
28  : m_lineConstrained(false)
29  , m_originConstrained(false)
30 {
31 }
32 
34  CDCObservations2D& observations2D) const
35 {
36 
37  if (observations2D.getNObservationsWithDriftRadius() > 0) {
38  updateWithDriftLength(trajectory2D, observations2D);
39  } else {
40  updateWithoutDriftLength(trajectory2D, observations2D);
41  }
42 
43  EForwardBackward isCoaligned = observations2D.isCoaligned(trajectory2D);
44  if (isCoaligned == EForwardBackward::c_Backward) trajectory2D.reverse();
45 }
46 
47 
48 
50  CDCObservations2D& observations2D) const
51 {
52  EigenObservationMatrix eigenObservation = getEigenObservationMatrix(&observations2D);
53  assert(eigenObservation.cols() == 4);
54  size_t nObservations = observations2D.size();
55 
56  if (isLineConstrained()) {
57 
58  //do a normal line fit
59  Eigen::Matrix<double, Eigen::Dynamic, 2> points = eigenObservation.leftCols<2>();
60 
61  Eigen::Matrix<double, 1, 2> pointMean;
62  //RowVector2d pointMean;
63  // cppcheck-suppress constStatement
64  pointMean << 0.0, 0.0;
65  if (!(isOriginConstrained())) {
66  // subtract the offset from the origin
67  pointMean = points.colwise().mean();
68 
69  points = points.rowwise() - pointMean;
70  }
71  Eigen::Matrix<double, 2, 2> covarianceMatrix = points.transpose() * points;
72 
73  Eigen::SelfAdjointEigenSolver< Eigen::Matrix<double, 2, 2> > eigensolver(covarianceMatrix);
74 
75  if (eigensolver.info() != Eigen::Success)
76  B2WARNING("SelfAdjointEigenSolver could not compute the eigen values of the observation matrix");
77 
78  //the eigenvalues are generated in increasing order
79  //we are interested in the lowest one since we want to compute the normal vector of the plane
80 
81  Eigen::Matrix<double, 2, 1> normalToLine = eigensolver.eigenvectors().col(0);
82 
83  double offset = -pointMean * normalToLine;
84 
85  // set the generalized circle parameters
86  // last set to zero constrains to a line
87  trajectory2D.setGlobalCircle(UncertainPerigeeCircle(PerigeeCircle::fromN(offset, normalToLine(0), normalToLine(1), 0)));
88 
89  } else {
90 
91  //lift the points to the projection space
92  Eigen::Matrix<double, Eigen::Dynamic, 3> projectedPoints(nObservations, 3);
93 
94  projectedPoints.col(0) = eigenObservation.col(0);
95  projectedPoints.col(1) = eigenObservation.col(1);
96  projectedPoints.col(2) = eigenObservation.leftCols<2>().rowwise().squaredNorm();
97 
98 
99  Eigen::Matrix<double, 1, 3> pointMean;
100  // cppcheck-suppress constStatement
101  pointMean << 0.0, 0.0, 0.0;
102  if (!(isOriginConstrained())) {
103  // subtract the offset from the origin
104  pointMean = projectedPoints.colwise().mean();
105 
106  projectedPoints = projectedPoints.rowwise() - pointMean;
107  }
108 
109  Eigen::Matrix<double, 3, 3> covarianceMatrix = projectedPoints.transpose() * projectedPoints;
110 
111  Eigen::SelfAdjointEigenSolver< Eigen::Matrix<double, 3, 3> > eigensolver(covarianceMatrix);
112 
113  if (eigensolver.info() != Eigen::Success)
114  B2WARNING("Eigen::SelfAdjointEigenSolver could not compute the eigen values of the observation matrix");
115 
116  //the eigenvalues are generated in increasing order
117  //we are interested in the lowest one since we want to compute the normal vector of the plane
118 
119  Eigen::Matrix<double, 3, 1> normalToPlane = eigensolver.eigenvectors().col(0);
120 
121  double offset = -pointMean * normalToPlane;
122 
123  trajectory2D.setGlobalCircle(UncertainPerigeeCircle(PerigeeCircle::fromN(offset, normalToPlane(0), normalToPlane(1),
124  normalToPlane(2))));
125  //fit.setParameters();
126 
127  }
128 
129  //check if the orientation is alright
130  Vector2D directionAtCenter = trajectory2D.getFlightDirection2D(Vector2D(0.0, 0.0));
131 
132 
133  size_t voteForChangeSign = 0;
134  for (size_t iPoint = 0; iPoint < nObservations; ++iPoint) {
135  double pointInSameDirection = eigenObservation(iPoint, 0) * directionAtCenter.x() +
136  eigenObservation(iPoint, 1) * directionAtCenter.y();
137  if (pointInSameDirection < 0) ++voteForChangeSign;
138  }
139 
140  if (voteForChangeSign > nObservations / 2.0) trajectory2D.reverse();
141 
142 }
143 
144 
145 
147 {
148  EigenObservationMatrix eigenObservation = getEigenObservationMatrix(&observations2D);
149  assert(eigenObservation.cols() == 4);
150  size_t nObservations = observations2D.size();
151 
152  //cout << "updateWithRightLeft : " << endl;
153  //observations always have the structure
154  /*
155  observables[0][0] == x of first point
156  observables[0][1] == y of first point
157  observables[0][2] == desired distance of first point
158  */
159 
160  Eigen::Matrix< double, Eigen::Dynamic, 1 > distances = eigenObservation.col(2);
161 
162  //cout << "distances : " << endl << distances << endl;
163 
164  if ((isLineConstrained()) && (isOriginConstrained())) {
165 
166 
167  Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > projectedPoints(nObservations, 2);
168 
169  //all coordiates
170  //projectedPoints.col(0) = Eigen::Matrix<double,Eigen::Dynamic,1>::Constant(nObservations,1.0);
171  projectedPoints.col(0) = eigenObservation.col(0);
172  projectedPoints.col(1) = eigenObservation.col(1);
173  //projectedPoints.col(2) = eigenObservation.leftCols<2>().rowwise().squaredNorm() - distances.rowwise().squaredNorm();
174 
175  Eigen::Matrix<double, 2, 1> parameters =
176  projectedPoints.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(distances);
177 
178  trajectory2D.setGlobalCircle(UncertainPerigeeCircle(PerigeeCircle::fromN(0.0, parameters(0), parameters(1), 0.0)));
179  }
180 
181  else if ((! isLineConstrained()) && (isOriginConstrained())) {
182 
183  Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > projectedPoints(nObservations, 3);
184 
185  //all coordiates
186  //projectedPoints.col(0) = Eigen::Matrix<double,Eigen::Dynamic,1>::Constant(1.0);
187  projectedPoints.col(0) = eigenObservation.col(0);
188  projectedPoints.col(1) = eigenObservation.col(1);
189  projectedPoints.col(2) = eigenObservation.leftCols<2>().rowwise().squaredNorm() - distances.rowwise().squaredNorm();
190 
191  Eigen::Matrix<double, 3, 1> parameters =
192  projectedPoints.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(distances);
193 
194  trajectory2D.setGlobalCircle(UncertainPerigeeCircle(PerigeeCircle::fromN(0.0, parameters(0), parameters(1), parameters(2))));
195  }
196 
197  else if ((isLineConstrained()) && (! isOriginConstrained())) {
198 
199  Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > projectedPoints(nObservations, 3);
200 
201  //all coordiates
202  projectedPoints.col(0) = Eigen::Matrix<double, Eigen::Dynamic, 1>::Constant(nObservations, 1.0);
203  projectedPoints.col(1) = eigenObservation.col(0);
204  projectedPoints.col(2) = eigenObservation.col(1);
205  //projectedPoints.col(3) = eigenObservation.leftCols<2>().rowwise().squaredNorm()- distances.rowwise().squaredNorm();
206 
207  Eigen::Matrix<double, 3, 1> parameters =
208  projectedPoints.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(distances);
209 
210  trajectory2D.setGlobalCircle(UncertainPerigeeCircle(PerigeeCircle::fromN(parameters(0), parameters(1), parameters(2), 0.0)));
211  //fit.setParameters(parameters(0),parameters(1),parameters(2),0.0);
212 
213  }
214 
215  else if ((! isLineConstrained()) && (! isOriginConstrained())) {
216 
217  Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > projectedPoints(nObservations, 4);
218 
219  //all coordiates
220  projectedPoints.col(0) = Eigen::Matrix<double, Eigen::Dynamic, 1>::Constant(nObservations, 1.0);
221  projectedPoints.col(1) = eigenObservation.col(0);
222  projectedPoints.col(2) = eigenObservation.col(1);
223  projectedPoints.col(3) = eigenObservation.leftCols<2>().rowwise().squaredNorm() - distances.rowwise().squaredNorm();
224 
225  //cout << "points : " << endl << projectedPoints << endl;
226 
227  Eigen::Matrix<double, 4, 1> parameters =
228  projectedPoints.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(distances);
229 
230  trajectory2D.setGlobalCircle(UncertainPerigeeCircle(PerigeeCircle::fromN(parameters(0), parameters(1), parameters(2),
231  parameters(3))));
232 
233  }
234 }
Belle2::TrackFindingCDC::CDCObservations2D::size
std::size_t size() const
Returns the number of observations stored.
Definition: CDCObservations2D.h:88
Belle2::TrackFindingCDC::CDCObservations2D::isCoaligned
EForwardBackward isCoaligned(const CDCTrajectory2D &trajectory2D) const
Checks if the last observation in the vector lies at greater or lower travel distance than the last o...
Definition: CDCObservations2D.h:305
Belle2::TrackFindingCDC::RiemannsMethod::updateWithDriftLength
void updateWithDriftLength(CDCTrajectory2D &trajectory2D, CDCObservations2D &observations2D) const
Executes the fit using the drift length information.
Definition: RiemannsMethod.cc:146
Belle2::TrackFindingCDC::Vector2D
A two dimensional vector which is equipped with functions for correct handeling of orientation relat...
Definition: Vector2D.h:37
Belle2::TrackFindingCDC::PerigeeCircle::fromN
static PerigeeCircle fromN(double n0, double n1, double n2, double n3=0)
Constructor with the four parameters of the generalized circle.
Belle2::TrackFindingCDC::Vector2D::y
double y() const
Getter for the y coordinate.
Definition: Vector2D.h:619
Belle2::TrackFindingCDC::CDCTrajectory2D::reverse
void reverse()
Reverses the trajectory in place.
Definition: CDCTrajectory2D.cc:97
Belle2::TrackFindingCDC::RiemannsMethod::isOriginConstrained
bool isOriginConstrained() const
Getter for the indictor that curves through the origin should be fitted by this fitter.
Definition: RiemannsMethod.h:64
Belle2::TrackFindingCDC::NForwardBackward::EForwardBackward
EForwardBackward
Enumeration to represent the distinct possibilities of the right left passage information.
Definition: EForwardBackward.h:35
Belle2::TrackFindingCDC::CDCTrajectory2D::setGlobalCircle
void setGlobalCircle(const UncertainPerigeeCircle &perigeeCircle)
Setter for the generalized circle that describes the trajectory.
Definition: CDCTrajectory2D.h:459
Belle2::TrackFindingCDC::CDCTrajectory2D
Particle trajectory as it is seen in xy projection represented as a circle.
Definition: CDCTrajectory2D.h:46
Belle2::TrackFindingCDC::RiemannsMethod::RiemannsMethod
RiemannsMethod()
Default constructor.
Definition: RiemannsMethod.cc:27
Belle2::TrackFindingCDC::UncertainPerigeeCircle
Adds an uncertainty matrix to the circle in perigee parameterisation.
Definition: UncertainPerigeeCircle.h:39
Belle2::TrackFindingCDC::CDCObservations2D
Class serving as a storage of observed drift circles to present to the Riemann fitter.
Definition: CDCObservations2D.h:53
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TrackFindingCDC::CDCTrajectory2D::getFlightDirection2D
Vector2D getFlightDirection2D(const Vector2D &point) const
Get the unit direction of flight at the given point, where arcLength2D = 0.
Definition: CDCTrajectory2D.h:320
Belle2::TrackFindingCDC::RiemannsMethod::update
void update(CDCTrajectory2D &trajectory2D, CDCObservations2D &observations2D) const
Executes the fit and updates the trajectory parameters. This may render the information in the observ...
Definition: RiemannsMethod.cc:33
Belle2::TrackFindingCDC::Vector2D::x
double x() const
Getter for the x coordinate.
Definition: Vector2D.h:609
Belle2::TrackFindingCDC::CDCObservations2D::getNObservationsWithDriftRadius
std::size_t getNObservationsWithDriftRadius() const
Returns the number of observations having a drift radius radius.
Definition: CDCObservations2D.cc:316
Belle2::TrackFindingCDC::RiemannsMethod::isLineConstrained
bool isLineConstrained() const
Getter for the indictor that lines should be fitted by this fitter.
Definition: RiemannsMethod.h:58
Belle2::TrackFindingCDC::RiemannsMethod::updateWithoutDriftLength
void updateWithoutDriftLength(CDCTrajectory2D &trajectory2D, CDCObservations2D &observations2D) const
Executes the fit without using the drift length information.
Definition: RiemannsMethod.cc:49