24 const int nTriplets = measurements.size() - 2;
26 if (nTriplets < 1)
return 0;
28 double combinedChi2 = 0.;
41 for (
int i = 0; i < nTriplets; i++) {
44 const B2Vector3D& hit0 = measurements.at(i)->getPosition();
45 const B2Vector3D& hit1 = measurements.at(i + 1)->getPosition();
46 const B2Vector3D& hit2 = measurements.at(i + 2)->getPosition();
48 const double d01sq = pow(hit1.
X() - hit0.
X(), 2) + pow(hit1.
Y() - hit0.
Y(), 2);
49 const double d12sq = pow(hit2.
X() - hit1.
X(), 2) + pow(hit2.
Y() - hit1.
Y(), 2);
50 const double d02sq = pow(hit2.
X() - hit0.
X(), 2) + pow(hit2.
Y() - hit0.
Y(), 2);
52 const double d01 =
sqrt(d01sq);
53 const double d12 =
sqrt(d12sq);
54 const double d02 =
sqrt(d02sq);
56 const double z01 = hit1.
Z() - hit0.
Z();
57 const double z12 = hit2.
Z() - hit1.
Z();
60 const double R_C = (d01 * d12 * d02) /
sqrt(-d01sq * d01sq - d12sq * d12sq - d02sq * d02sq +
61 2 * d01sq * d12sq + 2 * d12sq * d02sq + 2 * d02sq * d01sq);
64 const double Phi1C = 2. * asin(d01 / (2. * R_C));
65 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));
73 const double theta1C = acos(z01 / (Phi1C * R3D1C));
74 const double theta2C = acos(z12 / (Phi2C * R3D2C));
75 const double theta = (theta1C + theta2C) / 2.;
78 double alpha1 = R_C * R_C * Phi1C * Phi1C + z01 * z01;
79 alpha1 *= 1. / (0.5 * R_C * R_C * Phi1C * Phi1C * Phi1C /
tan(Phi1C / 2.) + z01 * z01);
80 double alpha2 = R_C * R_C * Phi2C * Phi2C + z12 * z12;
81 alpha2 *= 1. / (0.5 * R_C * R_C * Phi2C * Phi2C * Phi2C /
tan(Phi2C / 2.) + z12 * z12);
84 const double PhiTilde = - 0.5 * (Phi1C * alpha1 + Phi2C * alpha2);
86 const double eta = 0.5 * (Phi1C * alpha1 / R3D1C + Phi2C * alpha2 / R3D2C);
88 const double tanTheta1C =
tan(theta1C);
89 const double tanTheta2C =
tan(theta2C);
91 const double ThetaTilde = theta2C - theta1C - (1 - alpha2) / tanTheta2C + (1 - alpha1) / tanTheta1C;
93 const double beta = (1 - alpha2) / (R3D2C * tanTheta2C) - (1 - alpha1) / (R3D1C * tanTheta1C);
97 double entranceAngle = theta;
99 int detID = measurements.at(i + 1)->getType();
111 if (sensorOrigin.
Angle(normal) > M_PI * 0.5) { normal *= -1.; }
112 entranceAngle = (hit1 - hit0).Angle(normal);
121 const double sinTheta = sin(theta);
123 double R3D = - (eta * PhiTilde * sinTheta * sinTheta + beta * ThetaTilde);
124 R3D *= 1. / (eta * eta * sinTheta * sinTheta + beta * beta);
132 double R3D_truncated = R3D > R3DmaxCut ? R3DmaxCut : R3D;
133 const double sigmaMS = b / R3D_truncated;
136 double Chi2min = pow(beta * PhiTilde - eta * ThetaTilde, 2);
137 Chi2min *= 1. / (sigmaMS * sigmaMS * (eta * eta + beta * beta / (sinTheta * sinTheta)));
140 double sigmaR3DSquared = sigmaMS * sigmaMS / (eta * eta * sinTheta * sinTheta + beta * beta);
144 double delta = (beta * PhiTilde - eta * ThetaTilde) / (eta * PhiTilde * sinTheta * sinTheta + beta * ThetaTilde);
145 if (8 * delta * delta * sinTheta * sinTheta <= 1) {
147 R3D *= (0.75 +
sqrt(1 - 8 * delta * delta * sinTheta * sinTheta) / 4.);
152 m_alphas.push_back((alpha1 + alpha2) / 2.);
159 combinedChi2 += Chi2min;
163 double numerator = 0;
164 double denominator = 0;
166 for (
short i = 0; i < nTriplets; ++i) {
172 numerator += addValue *
m_R3Ds.at(i);
173 denominator += addValue;
179 double globalCompatibilityOfR3Ds = 0;
180 for (
short i = 0; i < nTriplets; ++i) {
184 double const finalChi2 = combinedChi2 + globalCompatibilityOfR3Ds;
188 return TMath::Prob(finalChi2, 2 * measurements.size() - 5);