Belle II Software  release-05-01-25
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 32 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 91 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 22 of file QualityEstimatorRiemannHelixFit.cc.

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 }

◆ 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 79 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 64 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 101 of file QualityEstimatorBase.h.

The documentation for this class was generated from the following files:
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::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