Belle II Software  release-05-01-25
QualityEstimatorRiemannHelixFit.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2017 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Felix Metzner, Jonas Wagner *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include "tracking/trackFindingVXD/trackQualityEstimators/QualityEstimatorRiemannHelixFit.h"
12 #include <tracking/trackFindingVXD/utilities/CalcCurvatureSignum.h>
13 #include <math.h>
14 #include <Eigen/Dense>
15 #include <framework/logging/Logger.h>
16 #include <TMath.h>
17 
18 using namespace Belle2;
19 // Set precision to be used for Eigen Library
20 typedef double Precision;
21 
22 double QualityEstimatorRiemannHelixFit::estimateQuality(std::vector<SpacePoint const*> const& measurements)
23 {
24  const int nHits = measurements.size();
25  if (nHits < 3) return 0;
26 
27  // Circle Fit
28 
29  Eigen::Matrix<Precision, Eigen::Dynamic, Eigen::Dynamic> W = Eigen::Matrix<Precision, Eigen::Dynamic, Eigen::Dynamic>::Zero(nHits,
30  nHits);
31  Eigen::Matrix<Precision, Eigen::Dynamic, Eigen::Dynamic> Wz = Eigen::Matrix<Precision, Eigen::Dynamic, Eigen::Dynamic>::Zero(nHits,
32  nHits);
33  Eigen::Matrix<Precision, Eigen::Dynamic, 3> X = Eigen::Matrix<Precision, Eigen::Dynamic, 3>::Zero(nHits, 3);
34  Eigen::Matrix<Precision, Eigen::Dynamic, 1> Z = Eigen::Matrix<Precision, Eigen::Dynamic, 1>::Zero(nHits, 1);
35  Precision traceOfW = 0.;
36 
37  short index = 0;
38  for (SpacePoint const* hit : measurements) {
39  auto position = hit->getPosition();
40  auto positionError = hit->getPositionError();
41  double x = position.X();
42  double y = position.Y();
43  double z = position.Z();
44  double sigmaX = positionError.X();
45  double sigmaY = positionError.Y();
46  double sigmaZ = positionError.Z();
47 
48  double r2 = x * x + y * y;
49  double inverseVarianceXY = 1. / sqrt(sigmaX * sigmaX + sigmaY * sigmaY);
50 
51  // The following weight matrix W can be improved for cases with multiple scattering
52  // by adding a correction term which will make the matrix non-diagonal.
53  // However, this requires prior knowledge of the curvature of the track and thus a
54  // second iteration (see page 368 of above mentioned source).
55  W(index, index) = inverseVarianceXY;
56  traceOfW += inverseVarianceXY;
57 
58  X(index, 0) = x;
59  X(index, 1) = y;
60  X(index, 2) = r2;
61 
62  // Values for z line fit for extrended Riemann fit
63  Wz(index, index) = 1. / sigmaZ;
64  Z(index, 0) = z;
65 
66  index++;
67  }
68 
69  Eigen::Matrix<Precision, 1, 3> xBar = Eigen::Matrix<Precision, Eigen::Dynamic, 1>::Ones(nHits, 1).transpose() * W * X / traceOfW;
70  Eigen::Matrix<Precision, 3, 3> Vx = X.transpose() * W * X - xBar.transpose() * xBar * traceOfW;
71 
72  // Find eigenvector to smallest eigenvalue
73  Eigen::EigenSolver<Eigen::Matrix<Precision, 3, 3>> eigencollection(Vx);
74  Eigen::Matrix<Precision, 3, 1> eigenvalues = eigencollection.eigenvalues().real();
75  Eigen::Matrix<std::complex<Precision>, 3, 3> eigenvectors = eigencollection.eigenvectors();
76  Eigen::Matrix<Precision, 3, 1>::Index minRow, minCol;
77  eigenvalues.minCoeff(&minRow, &minCol);
78 
79  Eigen::Matrix<Precision, 3, 1> n = eigenvectors.col(minRow).real();
80 
81  // Calculate results with this eigenvector
82  Precision c = - xBar * n;
83  Precision x0 = - 0.5 * n(0) / n(2);
84  Precision y0 = - 0.5 * n(1) / n(2);
85  Precision rho2 = (1 - n(2) * (n(2) + 4 * c)) / (4 * n(2) * n(2));
86  Precision rho = sqrt(rho2);
87 
88  // calculation of chi2 for circle fit using Karimaeki circle fit
89  Precision divisor = 1. / traceOfW;
90  Eigen::Matrix<Precision, Eigen::Dynamic, 1> unitvec = Eigen::Matrix<Precision, Eigen::Dynamic, 1>::Ones(nHits, 1);
91  Precision meanX = unitvec.transpose() * W * X.col(0);
92  meanX *= divisor;
93  Precision meanY = unitvec.transpose() * W * X.col(1);
94  meanY *= divisor;
95  Precision meanXY = unitvec.transpose() * W * (X.col(0).cwiseProduct(X.col(1)));
96  meanXY *= divisor;
97  Precision meanX2 = unitvec.transpose() * W * (X.col(0).cwiseProduct(X.col(0)));
98  meanX2 *= divisor;
99  Precision meanY2 = unitvec.transpose() * W * (X.col(1).cwiseProduct(X.col(1)));
100  meanY2 *= divisor;
101  Precision meanXR2 = unitvec.transpose() * W * (X.col(0).cwiseProduct(X.col(2)));
102  meanXR2 *= divisor;
103  Precision meanYR2 = unitvec.transpose() * W * (X.col(1).cwiseProduct(X.col(2)));
104  meanYR2 *= divisor;
105  Precision meanR2 = unitvec.transpose() * W * X.col(2);
106  meanR2 *= divisor;
107  Precision meanR4 = unitvec.transpose() * W * (X.col(2).cwiseProduct(X.col(2)));
108  meanR4 *= divisor;
109 
110  // covariances:
111  Precision covXX = meanX2 - meanX * meanX;
112  Precision covXY = meanXY - meanX * meanY;
113  Precision covYY = meanY2 - meanY * meanY;
114  Precision covXR2 = meanXR2 - meanX * meanR2;
115  Precision covYR2 = meanYR2 - meanY * meanR2;
116  Precision covR2R2 = meanR4 - meanR2 * meanR2;
117 
118  // q1, q2: helping variables, to make the code more readable
119  Precision q1 = covR2R2 * covXY - covXR2 * covYR2;
120  Precision q2 = covR2R2 * (covXX - covYY) - covXR2 * covXR2 + covYR2 * covYR2;
121 
122  Precision pocaPhi = 0.5 * atan2(2. * q1, q2);
123 
124  Precision sinPhi = sin(pocaPhi);
125  Precision cosPhi = cos(pocaPhi);
126  Precision kappa = (sinPhi * covXR2 - cosPhi * covYR2) / covR2R2;
127  Precision delta = -kappa * meanR2 + sinPhi * meanX - cosPhi * meanY;
128  Precision rootTerm = sqrt(1. - 4.*delta * kappa);
129  Precision curvature = 2.*kappa / (rootTerm);
130  Precision pocaD = 2.*delta / (1. + rootTerm);
131  short curvatureSign = calcCurvatureSignum(measurements);
132 
133  if ((curvature < 0 && curvatureSign >= 0) || (curvature > 0 && curvatureSign < 0)) {
134  curvature = -curvature;
135  // pocaPhi = pocaPhi + M_PI; //
136  pocaD = -pocaD;
137  }
138 
139  Precision chi2 = traceOfW * (1. + pocaD / rho) * (1. + curvature * pocaD) *
140  (sinPhi * sinPhi * covXX - 2.*sinPhi * cosPhi * covXY + cosPhi * cosPhi * covYY - kappa * kappa * covR2R2);
141 
142  // Arc length calculation
143  Precision x_first = X.col(0)(0) - x0;
144  Precision y_first = X.col(1)(0) - y0;
145  Precision r_mag_first = sqrt(x_first * x_first + y_first * y_first);
146  Eigen::Matrix<Precision, Eigen::Dynamic, 1> x0s = Eigen::Matrix<Precision, Eigen::Dynamic, 1>::Ones(nHits, 1) * x0;
147  Eigen::Matrix<Precision, Eigen::Dynamic, 1> y0s = Eigen::Matrix<Precision, Eigen::Dynamic, 1>::Ones(nHits, 1) * y0;
148  Eigen::Matrix<Precision, Eigen::Dynamic, 1> r_mags = ((X.col(0) - x0s).cwiseProduct(X.col(0) - x0s) + (X.col(1) - y0s).cwiseProduct(
149  X.col(1) - y0s)).cwiseSqrt();
150 
151  Eigen::Matrix<Precision, Eigen::Dynamic, 1> arc_angles = (x_first * (X.col(0) - x0s) + y_first * (X.col(1) - y0s)).cwiseQuotient(
152  r_mags) / r_mag_first;
153  Eigen::Matrix<Precision, Eigen::Dynamic, 1> arc_lengths = rho * arc_angles.array().acos().matrix();
154 
155  Eigen::Matrix<Precision, Eigen::Dynamic, 2> A = Eigen::Matrix<Precision, Eigen::Dynamic, 2>::Ones(nHits, 2);
156  arc_lengths(0) = 0.;
157  A.col(1) = arc_lengths;
158 
159  // Linear regression of z on arg_lengths with weight matrix Wz
160  Eigen::Matrix<Precision, 2, 1> p = (A.transpose() * Wz * A).inverse() * A.transpose() * Wz * Z;
161 
162  // Construction of momentum vector with innermost hit and reconstructed circle center
163  Eigen::Matrix<Precision, 3, 1> momVec = Eigen::Matrix<Precision, 3, 1>::Zero();
164  momVec(0) = y0 - X.col(1)(0);
165  momVec(1) = - (x0 - X.col(0)(0));
166 
167  Precision pT = Precision(calcPt(rho));
168  momVec = pT * momVec.normalized();
169 
170  Eigen::Matrix<Precision, 3, 1> vec01 = X.row(1) - X.row(0);
171  vec01(2) = Z(1) - Z(0);
172  Precision angle01 = std::acos(vec01.dot(momVec) / momVec.norm() / vec01.norm());
173  if (angle01 > 0.5 * M_PI) { momVec *= -1.; }
174 
175  // Calculation of a chi2 distributed quantity for the quality of fit of the z component fit.
176  Eigen::Matrix<Precision, Eigen::Dynamic, 1> ones = Eigen::Matrix<Precision, Eigen::Dynamic, 1>::Ones(nHits, 1);
177  Precision chi2_z = ((Z - p(0) * ones - p(1) * arc_lengths).cwiseQuotient(Wz * ones)).transpose() * ((Z - p(0) * ones - p(
178  1) * arc_lengths).cwiseQuotient(Wz * ones));
179 
180  // Adding chi2 and chi2_z, thus creating a chi2 distribution with (n-3) + (n-2) = 2n-5 degrees of freedom
181  m_results.chiSquared = chi2 + chi2_z;
182  B2DEBUG(75, "Chi Squared of extended Riemann = " << * (m_results.chiSquared) << std::endl);
183 
184  Precision pZ = pT * p(1);
185  momVec(2) = pZ;
186  m_results.pt = pT;
187  m_results.p = TVector3(momVec(0), momVec(1), momVec(2));
188  m_results.curvatureSign = curvatureSign;
189 
190  return TMath::Prob(*(m_results.chiSquared), 2 * measurements.size() - 5);
191 }
192 
Belle2::QualityEstimationResults::p
boost::optional< B2Vector3< double > > p
momentum vector estimate from the QE
Definition: QualityEstimatorBase.h:42
Belle2::QualityEstimationResults::chiSquared
boost::optional< double > chiSquared
chi squared value obtained by the fit of the QE
Definition: QualityEstimatorBase.h:38
Belle2::calcCurvatureSignum
short calcCurvatureSignum(std::vector< SpacePoint const * > const &measurements)
Calculate curvature based on triplets of measurements.
Definition: CalcCurvatureSignum.h:32
Belle2::QualityEstimationResults::pt
boost::optional< double > pt
transverse momentum estimate from the QE
Definition: QualityEstimatorBase.h:40
Belle2::SpacePoint
SpacePoint typically is build from 1 PXDCluster or 1-2 SVDClusters.
Definition: SpacePoint.h:52
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::QualityEstimatorRiemannHelixFit::estimateQuality
virtual double estimateQuality(std::vector< SpacePoint const * > const &measurements) final
Minimal implementation of the quality estimation Calculates quality indicator in range [0,...
Definition: QualityEstimatorRiemannHelixFit.cc:22
Belle2::QualityEstimatorBase::calcPt
double calcPt(double const radius) const
Returns a value for the transverse momentum in GeV calculated from a provided radius.
Definition: QualityEstimatorBase.h:91
Belle2::QualityEstimationResults::curvatureSign
boost::optional< short > curvatureSign
direction of curvature as obtained by the QE
Definition: QualityEstimatorBase.h:39
Belle2::QualityEstimatorBase::m_results
QualityEstimationResults m_results
Result of the quality estimation This is stored as a member variable, because some values may be calc...
Definition: QualityEstimatorBase.h:101