Minimal implementation of the quality estimation Calculates quality indicator in range [0,1].
measurements - std::vector<SpacePoint const*> ordered from innermost to outermost measurement
21{
22 const int nHits = measurements.size();
23 if (nHits < 3) return 0;
24
25
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;
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
50
51
52
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
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
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
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
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
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
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);
130
131 if ((curvature < 0 && curvatureSign >= 0) || (curvature > 0 && curvatureSign < 0)) {
132 curvature = -curvature;
133
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
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
158 Eigen::Matrix<Precision, 2, 1> p = (A.transpose() * Wz * A).inverse() * A.transpose() * Wz * Z;
159
160
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
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
180 B2DEBUG(25,
"Chi Squared of extended Riemann = " << * (
m_results.
chiSquared) << std::endl);
181
182 Precision pZ = pT * p(1);
183 momVec(2) = pZ;
188
190}
QualityEstimationResults m_results
Result of the quality estimation This is stored as a member variable, because some values may be calc...
double calcPt(double const radius) const
Returns a value for the transverse momentum in GeV calculated from a provided radius.
SpacePoint typically is build from 1 PXDCluster or 1-2 SVDClusters.
B2Vector3< double > B2Vector3D
typedef for common usage with double
double sqrt(double a)
sqrt for double
short calcCurvatureSignum(std::vector< SpacePoint const * > const &measurements)
Calculate curvature based on triplets of measurements.
std::optional< double > pt
transverse momentum estimate from the QE
std::optional< short > curvatureSign
direction of curvature as obtained by the QE
std::optional< double > pocaD
distance to the z-axis of the POCA
std::optional< double > chiSquared
chi squared value obtained by the fit of the QE
std::optional< B2Vector3D > p
momentum vector estimate from the QE
1.9.6