14 #include <framework/gearbox/Const.h>
16 #include <svd/geometry/SensorInfo.h>
18 #include "tracking/trackFindingVXD/trackQualityEstimators/QualityEstimatorTripletFit.h"
19 #include <tracking/trackFindingVXD/utilities/CalcCurvatureSignum.h>
26 const int nTriplets = measurements.size() - 2;
28 if (nTriplets < 1)
return 0;
30 double combinedChi2 = 0.;
43 for (
int i = 0; i < nTriplets; i++) {
46 const TVector3 hit0 = measurements.at(i)->getPosition();
47 const TVector3 hit1 = measurements.at(i + 1)->getPosition();
48 const TVector3 hit2 = measurements.at(i + 2)->getPosition();
50 const double d01sq = pow(hit1.X() - hit0.X(), 2) + pow(hit1.Y() - hit0.Y(), 2);
51 const double d12sq = pow(hit2.X() - hit1.X(), 2) + pow(hit2.Y() - hit1.Y(), 2);
52 const double d02sq = pow(hit2.X() - hit0.X(), 2) + pow(hit2.Y() - hit0.Y(), 2);
54 const double d01 = sqrt(d01sq);
55 const double d12 = sqrt(d12sq);
56 const double d02 = sqrt(d02sq);
58 const double z01 = hit1.Z() - hit0.Z();
59 const double z12 = hit2.Z() - hit1.Z();
61 const double R_C = (d01 * d12 * d02) / sqrt(-d01sq * d01sq - d12sq * d12sq - d02sq * d02sq + 2 * d01sq * d12sq + 2 * d12sq * d02sq +
65 const double Phi1C = 2. * asin(d01 / (2. * R_C));
66 const double Phi2C = 2. * asin(d12 / (2. * R_C));
69 const double R3D1C = sqrt(R_C * R_C + (z01 * z01) / (Phi1C * Phi1C));
70 const double R3D2C = sqrt(R_C * R_C + (z12 * z12) / (Phi2C * Phi2C));
72 const double theta1C = acos(z01 / (Phi1C * R3D1C));
73 const double theta2C = acos(z12 / (Phi2C * R3D2C));
74 const double theta = (theta1C + theta2C) / 2.;
76 double alpha1 = R_C * R_C * Phi1C * Phi1C + z01 * z01;
77 alpha1 *= 1. / (0.5 * R_C * R_C * Phi1C * Phi1C * Phi1C / tan(Phi1C / 2.) + z01 * z01);
78 double alpha2 = R_C * R_C * Phi2C * Phi2C + z12 * z12;
79 alpha2 *= 1. / (0.5 * R_C * R_C * Phi2C * Phi2C * Phi2C / tan(Phi2C / 2.) + z12 * z12);
81 const double PhiTilde = - 0.5 * (Phi1C * alpha1 + Phi2C * alpha2);
82 const double eta = 0.5 * Phi1C * alpha1 / R3D1C + 0.5 * Phi2C * alpha2 / R3D2C;
83 const double ThetaTilde = theta2C - theta1C - (1 - alpha2) / tan(theta2C) + (1 - alpha1) / tan(theta1C);
84 const double beta = (1 - alpha2) / (R3D2C * tan(theta2C)) - (1 - alpha1) / (R3D1C * tan(theta1C));
88 double entranceAngle = theta;
90 int detID = measurements.at(i + 1)->getType();
93 TVector3 sensorOrigin = sensor.pointToGlobal(TVector3(0, 0, 0),
true);
94 TVector3 sensoru = sensor.pointToGlobal(TVector3(1, 0, 0),
true);
95 TVector3 sensorv = sensor.pointToGlobal(TVector3(0, 1, 0),
true);
97 TVector3 globalu = sensoru - sensorOrigin;
98 TVector3 globalv = sensorv - sensorOrigin;
99 TVector3 normal = globalu.Cross(globalv);
102 if (sensorOrigin.Angle(normal) > M_PI * 0.5) { normal *= -1.; }
103 entranceAngle = (hit1 - hit0).Angle(normal);
112 double R3D = - (eta * PhiTilde * sin(theta) * sin(theta) + beta * ThetaTilde);
113 R3D *= 1. / (eta * eta * sin(theta) * sin(theta) + beta * beta);
121 double R3D_truncated = R3D > R3DmaxCut ? R3DmaxCut : R3D;
122 const double sigmaMS = b / R3D_truncated;
124 double sigmaR3DSquared = pow(sigmaMS, 2) / (pow(eta * sin(theta), 2) + pow(beta, 2));
126 double Chi2min = pow(beta * PhiTilde - eta * ThetaTilde, 2);
127 Chi2min *= 1. / (sigmaMS * sigmaMS * (eta * eta + beta * beta / pow(sin(theta), 2)));
130 double delta = (beta * PhiTilde - eta * ThetaTilde) / (eta * PhiTilde * sin(theta) * sin(theta) + beta * ThetaTilde);
131 if (8 * delta * delta * sin(theta) * sin(theta) <= 1) {
132 R3D *= (0.75 + sqrt(1 - 8 * delta * delta * sin(theta) * sin(theta)) / 4.);
137 m_alphas.push_back((alpha1 + alpha2) / 2.);
143 combinedChi2 += Chi2min;
147 double numerator = 0;
148 double denominator = 0;
149 for (
short i = 0; i < nTriplets; ++i) {
157 double globalCompatibilityOfR3Ds = 0;
158 for (
short i = 0; i < nTriplets; ++i) {
162 double const finalChi2 = combinedChi2 + globalCompatibilityOfR3Ds;
166 return TMath::Prob(finalChi2, 2 * measurements.size() - 5);
175 if (measurements.size() < 3)
return m_results;
180 double averageThetaPrime = 0;
181 for (
unsigned short i = 0; i <
m_thetas.size(); ++i) {
184 averageThetaPrime /=
m_thetas.size();