Belle II Software  release-08-01-10
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
QualityEstimatorRiemannHelixFit Class Reference
Packages » tracking

Based on R. More...

#include <QualityEstimatorRiemannHelixFit.h>

Inheritance diagram for QualityEstimatorRiemannHelixFit:
Inheritance graph
[legend]
Collaboration diagram for QualityEstimatorRiemannHelixFit:
Collaboration graph
[legend]

Public Member Functions

virtual double estimateQuality (std::vector< SpacePoint const * > const &measurements) final
 Minimal implementation of the quality estimation Calculates quality indicator in range [0,1]. More...
 
void setMagneticFieldStrength (double magneticFieldZ=1.5)
 Setter for z component of magnetic field. More...
 
virtual QualityEstimationResults estimateQualityAndProperties (std::vector< SpacePoint const * > const &measurements)
 Quality estimation providing additional quantities Calculates quality indicator in range [0,1] Optionally returns chi2 and additional informations. More...
 

Protected Member Functions

double calcPt (double const radius) const
 Returns a value for the transverse momentum in GeV calculated from a provided radius. More...
 

Protected Attributes

double m_magneticFieldZ = 1.5
 Member storing the z component of the magnetic field.
 
QualityEstimationResults m_results
 Result of the quality estimation This is stored as a member variable, because some values may be calculated by 'estimateQuality' anyways. More...
 

Detailed Description

Based on R.

Fruehwirth, A. Strandlie, W. Waltenberger, Nuclear instruments and Methods in Physics Research A 490 (2002) 366-378

Definition at line 22 of file QualityEstimatorRiemannHelixFit.h.

Member Function Documentation

◆ calcPt()

double calcPt ( double const  radius) const
inlineprotectedinherited

Returns a value for the transverse momentum in GeV calculated from a provided radius.

Utilizing m_magneticFieldZ and hard coded speed of light

Definition at line 80 of file QualityEstimatorBase.h.

◆ estimateQuality()

double estimateQuality ( std::vector< SpacePoint const * > const &  measurements)
finalvirtual

Minimal implementation of the quality estimation Calculates quality indicator in range [0,1].

measurements - std::vector<SpacePoint const*> ordered from innermost to outermost measurement

Implements QualityEstimatorBase.

Definition at line 20 of file QualityEstimatorRiemannHelixFit.cc.

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

◆ estimateQualityAndProperties()

virtual QualityEstimationResults estimateQualityAndProperties ( std::vector< SpacePoint const * > const &  measurements)
inlinevirtualinherited

Quality estimation providing additional quantities Calculates quality indicator in range [0,1] Optionally returns chi2 and additional informations.

Eg. momentum estimation.

measurements - std::vector<SpacePoint const*> ordered from innermost to outermost measurement

Reimplemented in QualityEstimatorMC, and QualityEstimatorTripletFit.

Definition at line 68 of file QualityEstimatorBase.h.

◆ setMagneticFieldStrength()

void setMagneticFieldStrength ( double  magneticFieldZ = 1.5)
inlineinherited

Setter for z component of magnetic field.

Parameters
magneticFieldZ: value to set it to

Definition at line 53 of file QualityEstimatorBase.h.

Member Data Documentation

◆ m_results

QualityEstimationResults m_results
protectedinherited

Result of the quality estimation This is stored as a member variable, because some values may be calculated by 'estimateQuality' anyways.

Therefore they don't need to be calculated explicitly in 'estimateQualityAndProperties'.

Definition at line 90 of file QualityEstimatorBase.h.

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