Belle II Software  release-05-01-25
CDCSZFitter.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2012 - 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/CDCSZFitter.h>
11 
12 #include <tracking/trackFindingCDC/fitting/EigenObservationMatrix.h>
13 #include <tracking/trackFindingCDC/fitting/CDCSZObservations.h>
14 #include <tracking/trackFindingCDC/fitting/CDCObservations2D.h>
15 
16 #include <tracking/trackFindingCDC/eventdata/tracks/CDCTrack.h>
17 #include <tracking/trackFindingCDC/eventdata/tracks/CDCSegmentPair.h>
18 
19 #include <tracking/trackFindingCDC/eventdata/segments/CDCSegment3D.h>
20 #include <tracking/trackFindingCDC/eventdata/segments/CDCSegment2D.h>
21 
22 #include <tracking/trackFindingCDC/eventdata/trajectories/CDCTrajectorySZ.h>
23 #include <tracking/trackFindingCDC/eventdata/trajectories/CDCTrajectory2D.h>
24 
25 #include <tracking/trackFindingCDC/geometry/SZParameters.h>
26 
27 #include <tracking/trackFindingCDC/numerics/EigenView.h>
28 
29 #include <framework/logging/Logger.h>
30 
31 #include <Eigen/Eigen>
32 #include <Eigen/QR>
33 #include <Eigen/Core>
34 
35 using namespace Belle2;
36 using namespace TrackFindingCDC;
37 
39 {
40  static CDCSZFitter szFitter;
41  return szFitter;
42 }
43 
44 namespace {
45  UncertainSZLine fitZ(const Eigen::Matrix<double, 3, 3>& sumMatrixWSZ)
46  {
47  // Solve the normal equation X * n = y
48  Eigen::Matrix<double, 2, 2> sumMatrixWS = sumMatrixWSZ.block<2, 2>(0, 0);
49  Eigen::Matrix<double, 2, 1> sumVectorZOverWS = sumMatrixWSZ.block<2, 1>(0, 2);
50  Eigen::Matrix<double, 2, 1> nZOverWS = sumMatrixWS.llt().solve(sumVectorZOverWS);
51  double chi2 = sumMatrixWSZ(2, 2) - nZOverWS.transpose() * sumVectorZOverWS;
52 
53  using namespace NSZParameterIndices;
54  SZParameters szParameters;
55  szParameters(c_TanL) = nZOverWS(1);
56  szParameters(c_Z0) = nZOverWS(0);
57 
58  SZPrecision szPrecision;
59  szPrecision(c_TanL, c_TanL) = sumMatrixWS(1, 1);
60  szPrecision(c_Z0, c_TanL) = sumMatrixWS(0, 1);
61  szPrecision(c_TanL, c_Z0) = sumMatrixWS(1, 0);
62  szPrecision(c_Z0, c_Z0) = sumMatrixWS(0, 0);
63 
64  SZCovariance szCovariance = SZUtil::covarianceFromFullPrecision(szPrecision);
65  return UncertainSZLine(szParameters, szCovariance, chi2);
66  }
67 
68  // Declare function as currently unused to avoid compiler warning
69  UncertainSZLine fitSZ(const Eigen::Matrix<double, 3, 3>& sumMatrixWSZ) __attribute__((__unused__));
70 
72  UncertainSZLine fitSZ(const Eigen::Matrix<double, 3, 3>& sumMatrixWSZ)
73  {
74 
75  // Matrix of averages
76  Eigen::Matrix<double, 3, 3> averageMatrixWSZ = sumMatrixWSZ / sumMatrixWSZ(0);
77 
78  // Measurement means
79  Eigen::Matrix<double, 3, 1> meansWSZ = averageMatrixWSZ.row(0);
80 
81  // Covariance matrix
82  Eigen::Matrix<double, 3, 3> covMatrixWSZ = averageMatrixWSZ - meansWSZ * meansWSZ.transpose();
83 
84  Eigen::Matrix<double, 2, 2> covMatrixSZ = covMatrixWSZ.block<2, 2>(1, 1);
85  Eigen::Matrix<double, 2, 1> meansSZ = meansWSZ.segment<2>(1);
86 
87  Eigen::SelfAdjointEigenSolver<Eigen::Matrix<double, 2, 2>> eigensolver(covMatrixSZ);
88  if (eigensolver.info() != Eigen::Success) {
89  B2WARNING(
90  "SelfAdjointEigenSolver could not compute the eigen values of the observation matrix");
91  }
92 
93  // the eigenvalues are generated in increasing order
94  // we are interested in the lowest one since we want to compute the normal vector of the line
95  // here
96  Eigen::Matrix<double, 2, 1> nSZ = eigensolver.eigenvectors().col(0);
97  B2INFO("nSZ " << nSZ);
98  nSZ /= std::copysign(nSZ.norm(), -nSZ(1)); // Making n2 negative to normalize to forward along s
99  Eigen::Matrix<double, 3, 1> nWSZ;
100  // cppcheck-suppress constStatement
101  nWSZ << -meansSZ.transpose() * nSZ, nSZ;
102  B2INFO("nWSZ " << nWSZ);
103 
104  double chi2 = nWSZ.transpose() * sumMatrixWSZ * nWSZ;
105 
106  Eigen::Matrix<double, 3, 3> precN = sumMatrixWSZ;
107 
108  using namespace NSZParameterIndices;
109  SZParameters szParameters;
110  szParameters(c_TanL) = -nWSZ(1) / nWSZ(2);
111  szParameters(c_Z0) = -nWSZ(0) / nWSZ(2);
112 
113  Eigen::Matrix<double, 3, 2> ambiguitySZToN = Eigen::Matrix<double, 3, 2>::Zero();
114  ambiguitySZToN(0, c_Z0) = -nWSZ(2);
115  ambiguitySZToN(0, c_TanL) = -nWSZ(1) * nWSZ(0);
116  ambiguitySZToN(1, c_TanL) = -nWSZ(2) * nWSZ(2) * nWSZ(2);
117  ambiguitySZToN(2, c_TanL) = nWSZ(1) * nWSZ(2) * nWSZ(2);
118 
119  SZPrecision szPrecision;
120  mapToEigen(szPrecision) = ambiguitySZToN.transpose() * precN * ambiguitySZToN;
121 
122  SZCovariance szCovariance = SZUtil::covarianceFromFullPrecision(szPrecision);
123  return UncertainSZLine(szParameters, szCovariance, chi2);
124  }
125 }
126 
128 {
129  const bool onlyStereo = true;
130  CDCSZObservations observationsSZ(EFitVariance::c_Proper, onlyStereo);
131  observationsSZ.appendRange(track);
132  if (observationsSZ.size() > 3) {
133  CDCTrajectorySZ szTrajectory;
134  update(szTrajectory, observationsSZ);
135  return szTrajectory;
136  } else {
138  }
139 }
140 
142  const CDCTrajectory2D& axialTrajectory2D) const
143 {
144  B2ASSERT("Expected stereo segment", not stereoSegment.isAxial());
145 
146  CDCTrajectorySZ trajectorySZ;
147  update(trajectorySZ, stereoSegment, axialTrajectory2D);
148  return trajectorySZ;
149 }
150 
152 {
153  CDCSZObservations observationsSZ;
154  observationsSZ.appendRange(segment3D);
155  return fit(std::move(observationsSZ));
156 }
157 
159 {
160  CDCTrajectorySZ trajectorySZ;
161  update(trajectorySZ, observationsSZ);
162  return trajectorySZ;
163 }
164 
165 void CDCSZFitter::update(const CDCSegmentPair& segmentPair) const
166 {
167  const CDCSegment2D* ptrStereoSegment = segmentPair.getStereoSegment();
168  const CDCSegment2D* ptrAxialSegment = segmentPair.getAxialSegment();
169 
170  assert(ptrStereoSegment);
171  assert(ptrAxialSegment);
172 
173  const CDCSegment2D& stereoSegment = *ptrStereoSegment;
174  const CDCSegment2D& axialSegment = *ptrAxialSegment;
175  const CDCTrajectory2D& axialTrajectory2D = axialSegment.getTrajectory2D();
176 
177  CDCTrajectorySZ trajectorySZ;
178  update(trajectorySZ, stereoSegment, axialTrajectory2D);
179 
180  CDCTrajectory3D trajectory3D(axialTrajectory2D, trajectorySZ);
181  segmentPair.setTrajectory3D(trajectory3D);
182 }
183 
185 {
186  CDCSZObservations szObservations;
187  for (size_t i = 0; i < observations2D.size(); ++i) {
188  const double s = observations2D.getX(i);
189  const double z = observations2D.getY(i);
190  szObservations.fill(s, z);
191  }
192  return fit(std::move(szObservations));
193 }
194 
196  const CDCSegment2D& stereoSegment,
197  const CDCTrajectory2D& axialTrajectory2D) const
198 {
199  B2ASSERT("Expected stereo segment", not stereoSegment.isAxial());
200 
201  // recostruct the stereo segment
202  CDCSZObservations observationsSZ;
203  for (const CDCRecoHit2D& recoHit2D : stereoSegment) {
204  CDCRecoHit3D recoHit3D = CDCRecoHit3D::reconstruct(recoHit2D, axialTrajectory2D);
205  observationsSZ.append(recoHit3D);
206  }
207 
208  update(trajectorySZ, observationsSZ);
209 }
210 
211 void CDCSZFitter::update(CDCTrajectorySZ& trajectorySZ, CDCSZObservations& observationsSZ) const
212 {
213  trajectorySZ.clear();
214  if (observationsSZ.size() < 3) {
215  return;
216  }
217 
218  // Determine NDF : Line fit eats up 2 degrees of freedom.
219  size_t ndf = observationsSZ.size() - 2;
220 
221  // Matrix of weighted sums
222  Eigen::Matrix<double, 3, 3> sumMatrixWSZ = getWSZSumMatrix(observationsSZ);
223  UncertainSZLine uncertainSZLine = fitZ(sumMatrixWSZ);
224 
225  uncertainSZLine.setNDF(ndf);
226  trajectorySZ.setSZLine(uncertainSZLine);
227 }
Belle2::TrackFindingCDC::CDCSegmentPair::getStereoSegment
const CDCSegment2D * getStereoSegment() const
Getter for the stereo segment.
Definition: CDCSegmentPair.h:144
Belle2::TrackFindingCDC::CDCObservations2D::size
std::size_t size() const
Returns the number of observations stored.
Definition: CDCObservations2D.h:88
Belle2::TrackFindingCDC::CDCRecoHit3D
Class representing a three dimensional reconstructed hit.
Definition: CDCRecoHit3D.h:62
Belle2::TrackFindingCDC::CDCTrack
Class representing a sequence of three dimensional reconstructed hits.
Definition: CDCTrack.h:51
Belle2::TrackFindingCDC::CDCObservations2D::getX
double getX(int iObservation) const
Getter for the x value of the observation at the given index.
Definition: CDCObservations2D.h:118
Belle2::TrackFindingCDC::CDCSegmentPair
Class representing a pair of one reconstructed axial segement and one stereo segment in adjacent supe...
Definition: CDCSegmentPair.h:44
Belle2::TrackFindingCDC::CDCSZFitter::getFitter
static const CDCSZFitter & getFitter()
Getter for a standard sz line fitter instance.
Definition: CDCSZFitter.cc:38
Belle2::TrackFindingCDC::UncertainSZLine
A line in sz where s is the transverse travel distance as seen in the xy projection with uncertaintie...
Definition: UncertainSZLine.h:38
Belle2::TrackFindingCDC::CDCTrajectorySZ::setSZLine
void setSZLine(const UncertainSZLine &szLine)
Setter for the line in sz space.
Definition: CDCTrajectorySZ.h:146
Belle2::TrackFindingCDC::UncertainSZLine::setNDF
void setNDF(std::size_t ndf)
Setter for the number of degrees of freediom used in the line fit.
Definition: UncertainSZLine.h:162
Belle2::TrackFindingCDC::CDCSZObservations::append
std::size_t append(const CDCRecoHit3D &recoHit3D)
Appends the observed position.
Definition: CDCSZObservations.cc:40
Belle2::TrackFindingCDC::CDCTrajectorySZ::basicAssumption
static CDCTrajectorySZ basicAssumption()
Constucts a basic assumption, what the z0 start position and the sz slope are, including some broad v...
Definition: CDCTrajectorySZ.cc:29
Belle2::TrackFindingCDC::CDCSZObservations::size
std::size_t size() const
Returns the number of observations stored.
Definition: CDCSZObservations.h:55
Belle2::TrackFindingCDC::NSZParameterIndices::c_TanL
@ c_TanL
Constant to address the tan lambda dip out of the xy plane.
Definition: SZParameters.h:36
Belle2::TrackFindingCDC::CDCTrajectorySZ
Linear trajectory in sz space.
Definition: CDCTrajectorySZ.h:41
Belle2::TrackFindingCDC::CDCSegment3D
A segment consisting of three dimensional reconstructed hits.
Definition: CDCSegment3D.h:36
Belle2::TrackFindingCDC::CDCTrajectory2D
Particle trajectory as it is seen in xy projection represented as a circle.
Definition: CDCTrajectory2D.h:46
Belle2::TrackFindingCDC::NSZParameterIndices::c_Z0
@ c_Z0
Constant to address the z start position.
Definition: SZParameters.h:39
Belle2::TrackFindingCDC::CDCTrajectorySZ::clear
void clear()
Clears all information from this trajectory line.
Definition: CDCTrajectorySZ.h:110
Belle2::TrackFindingCDC::CDCObservations2D::getY
double getY(int iObservation) const
Getter for the y value of the observation at the given index.
Definition: CDCObservations2D.h:124
Belle2::TrackFindingCDC::CDCSZFitter::fitWithStereoHits
CDCTrajectorySZ fitWithStereoHits(const CDCTrack &track) const
Returns the fitted sz trajectory of the track with the z-information of all stereo hits of the number...
Definition: CDCSZFitter.cc:127
Belle2::TrackFindingCDC::CDCSegment::isAxial
bool isAxial() const
Indicator if the underlying wires are axial.
Definition: CDCSegment.h:55
Belle2::TrackFindingCDC::CDCObservations2D
Class serving as a storage of observed drift circles to present to the Riemann fitter.
Definition: CDCObservations2D.h:53
Belle2::TrackFindingCDC::CDCRecoHit2D
Class representing a two dimensional reconstructed hit in the central drift chamber.
Definition: CDCRecoHit2D.h:57
Belle2::TrackFindingCDC::CDCSegment::getTrajectory2D
CDCTrajectory2D & getTrajectory2D() const
Getter for the two dimensional trajectory fitted to the segment.
Definition: CDCSegment.h:79
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::TrackFindingCDC::CDCSZFitter::update
void update(const CDCSegmentPair &segmentPair) const
Updates the trajectory of the axial stereo segment pair inplace.
Definition: CDCSZFitter.cc:165
Belle2::TrackFindingCDC::CDCSZFitter
Class implementing the z coordinate over travel distance line fit.
Definition: CDCSZFitter.h:37
Belle2::TrackFindingCDC::CDCSZObservations
Class serving as a storage of observed sz positions to present to the sz line fitters.
Definition: CDCSZObservations.h:36
Belle2::TrackFindingCDC::CDCSZObservations::appendRange
std::size_t appendRange(const std::vector< CDCRecoHit3D > &recoHit3Ds)
Appends all reconstructed hits from the three dimensional track.
Definition: CDCSZObservations.cc:75
Belle2::TrackFindingCDC::CDCSegment2D
A reconstructed sequence of two dimensional hits in one super layer.
Definition: CDCSegment2D.h:40
Belle2::TrackFindingCDC::UncertainParametersUtil< SZUtil, ESZParameter >::covarianceFromFullPrecision
static CovarianceMatrix covarianceFromFullPrecision(const PrecisionMatrix &prec)
Convert the precision matrix to the corresponding covariance matrix.
Definition: UncertainParameters.icc.h:101
Belle2::TrackFindingCDC::PlainMatrix
A matrix implementation to be used as an interface typ through out the track finder.
Definition: PlainMatrix.h:50
Belle2::TrackFindingCDC::CDCRecoHit3D::reconstruct
static CDCRecoHit3D reconstruct(const CDCRecoHit2D &recoHit2D, const CDCTrajectory2D &trajectory2D)
Reconstructs the three dimensional hit from the two dimensional and the two dimensional trajectory.
Definition: CDCRecoHit3D.cc:58
Belle2::TrackFindingCDC::CDCSegmentPair::setTrajectory3D
void setTrajectory3D(const CDCTrajectory3D &trajectory3D) const
Setter for the three dimensional trajectory.
Definition: CDCSegmentPair.h:207
Belle2::TrackFindingCDC::CDCTrajectory3D
Particle full three dimensional trajectory.
Definition: CDCTrajectory3D.h:47
Belle2::TrackFindingCDC::CDCSZFitter::fit
CDCTrajectorySZ fit(const CDCSegment2D &stereoSegment, const CDCTrajectory2D &axialTrajectory2D) const
Returns a fitted trajectory.
Definition: CDCSZFitter.cc:141
Belle2::TrackFindingCDC::CDCSegmentPair::getAxialSegment
const CDCSegment2D * getAxialSegment() const
Getter for the axial segment.
Definition: CDCSegmentPair.h:150
Belle2::TrackFindingCDC::CDCSZObservations::fill
std::size_t fill(double s, double z, double weight=1.0)
Appends the observed position.
Definition: CDCSZObservations.cc:25