12 using namespace Eigen;
14 static inline float sign(
float x) {
15 return ((
float)(x > 0.)) - ((float)(x < 0.));
23 Matrix<float, 2, 1>&
ha,
24 Matrix<float, 2, 1>& diff) {
29 Matrix<float, 2, 5>
H = Matrix<float, 2, 5>::Zero(2, 5);
30 Matrix<float, 2, 1>
ha = Matrix<float, 2, 1>::Zero(2, 1);
33 Matrix<float, 2, 1>
m = Matrix<float, 2, 1>::Zero(2, 1);
34 Matrix<float, 2, 2>
G = Matrix<float, 2, 2>::Zero(2, 2);
37 Matrix<float, 5, 5>
Q = Matrix<float, 5, 5>::Zero(5, 5);
40 Matrix<float, 5, 2> Ht = H.transpose();
42 Matrix<float, 5, 5> C_proj = state.
C +
Q;
44 Matrix<float, 5, 5> C_proj_inv = C_proj.fullPivLu().inverse();
46 Matrix<float, 5, 5> Cadd = Ht * G * H;
48 Matrix<float, 5, 5> Cinv = (C_proj_inv + Cadd);
50 state.
C = Cinv.fullPivLu().inverse();
52 Matrix<float, 5, 2> K = state.
C * Ht * G;
54 Matrix<float, 2, 1> proj_diff = Matrix<float, 2, 1>::Zero(2, 1);
57 Matrix<float, 5, 1> state_add = K * proj_diff;
59 Matrix<float, 5, 1> old_state = Matrix<float, 5, 1>::Zero(5, 1);
60 old_state(0, 0) = state.
phi;
61 old_state(1, 0) = state.
d;
63 old_state(2, 0) = state.
nu;
64 old_state(3, 0) = state.
z0;
65 old_state(4, 0) = state.
dzdl;
67 Matrix<float, 5, 1> new_state = old_state + state_add;
69 state.
phi = new_state(0, 0);
70 state.
d = new_state(1, 0);
71 state.
nu = new_state(2, 0);
73 state.
z0 = new_state(3, 0);
74 state.
dzdl = new_state(4, 0);
78 Matrix<float, 5, 1> state_diff = (new_state - old_state);
79 Matrix<float, 1, 5> state_diff_T = state_diff.transpose();
81 Matrix<float, 1, 1> chi2_add = (state_diff_T * C_proj_inv * state_diff);
82 state.
chi2 += chi2_add(0, 0);
89 Matrix<float, 5, 3>& dAdAp,
float& phi_p,
90 float& cosphi_p,
float& sinphi_p) {
97 float cosphi = cos(phi);
98 float sinphi = sin(phi);
105 float x0 = state.
x_int;
106 float y0 = state.
y_int;
108 phi_p = atan2((1. + k * d) * sinphi - k * y0, (1. + k * d) * cosphi - k * x0);
109 cosphi_p = cos(phi_p);
110 sinphi_p = sin(phi_p);
130 float tx = cosphi_p + k * x0;
131 float ty = sinphi_p + k * y0;
132 float tx2ty2_inv = 1. / (tx * tx + ty * ty);
133 float dphidtx = -ty * tx2ty2_inv;
134 float dphidty = tx * tx2ty2_inv;
135 float dtxdA_p = -sinphi_p;
136 float dtydA_p = cosphi_p;
137 dAdAp(0, 0) = dphidtx * dtxdA_p + dphidty * dtydA_p;
141 dAdAp(0, 1) = dphidtx * dtxdA_p + dphidty * dtydA_p;
144 float k_inv = 1. /
k;
145 float tx2ty2sqrtinv = sqrt(tx2ty2_inv);
146 float dddtx = tx2ty2sqrtinv * tx * k_inv;
147 float dddty = tx2ty2sqrtinv * ty * k_inv;
151 dAdAp(1, 0) = dddtx * dtxdA_p + dddty * dtydA_p;
155 dAdAp(1, 1) = dphidtx * dtxdA_p + dphidty * dtydA_p -
156 (sqrt(tx * tx + ty * ty) - 1.) * k_inv * k_inv;
159 dAdAp(1, 0) = y0 * cosphi_p - x0 * sinphi_p;
160 dAdAp(1, 1) = (dAdAp(1, 0)) * (dAdAp(1, 0)) * 0.5;
167 float sign_dzdl =
sign(dzdl);
169 float dx = d * cosphi;
170 float dy = d * sinphi;
171 float D = sqrt((x0 - dx) * (x0 - dx) + (y0 - dy) * (y0 - dy));
172 float D_inv = 1. /
D;
173 float v = 0.5 * k *
D;
175 float k_inv = 1. /
k;
176 float s = 2. * asin(v) * k_inv;
177 float s_inv = 1. /
s;
178 float dsdv = 2. / (k * sqrt(1 - v * v));
179 float dsdk = -s * k_inv;
183 float dz = sqrt(tmp1 / (1. - dzdl * dzdl));
184 float dz3 = dz * dz *
dz;
187 float ddxdA = -d * sinphi * dAdAp(0, 0) + cosphi * dAdAp(1, 0);
188 float ddydA = d * cosphi * dAdAp(0, 0) + sinphi * dAdAp(1, 0);
193 0.5 * D_inv * (2. * (dx - x0) * ddxdA + 2. * (dy - y0) * ddydA);
194 float dvdA = 0.5 * (k * dDdA + D * dkdA);
195 float dsdA = dsdv * dvdA + dsdk * dkdA;
196 float ddzdA = (dz * s_inv) * dsdA + (dz3 * tmp3) * ddzdldA;
197 dAdAp(3, 0) = -sign_dzdl * ddzdA;
205 dDdA = 0.5 * D_inv * (2. * (dx - x0) * ddxdA + 2. * (dy - y0) * ddydA);
206 dvdA = 0.5 * (k * dDdA + D * dkdA);
207 dsdA = dsdv * dvdA + dsdk * dkdA;
208 ddzdA = (dz * s_inv) * dsdA + (dz3 * tmp3) * ddzdldA;
209 dAdAp(3, 1) = -sign_dzdl * ddzdA;
217 dDdA = 0.5 * D_inv * (2. * (dx - x0) * ddxdA + 2. * (dy - y0) * ddydA);
218 dvdA = 0.5 * (k * dDdA + D * dkdA);
219 dsdA = dsdv * dvdA + dsdk * dkdA;
220 ddzdA = (dz * s_inv) * dsdA + (dz3 * tmp3) * ddzdldA;
221 dAdAp(3, 2) = -sign_dzdl * ddzdA;
224 float temp1 = k * D * 0.5;
226 float temp2 = D * 0.5;
231 s += (3. / 20.) *
temp2;
233 s += (5. / 56.) *
temp2;
234 float s_inv = 1. /
s;
235 float onedzdl2_inv = 1. / (1. - dzdl *
dzdl);
236 float dz = sqrt(s * s * dzdl * dzdl * onedzdl2_inv);
237 float dz_inv = 1. /
dz;
253 float ddxdA = -d * sinphi * dAdAp(0, 0) + cosphi * dAdAp(1, 0);
254 float ddydA = d * cosphi * dAdAp(0, 0) + sinphi * dAdAp(1, 0);
259 0.5 * D_inv * (2. * (dx - x0) * ddxdA + 2. * (dy - y0) * ddydA);
261 dsdA += (1. / 24.) * (3. * k2 * D2 * dDdA + 2. * k * D3 * dkdA);
262 dsdA += (3. / 640.) * (5. * D4 * k4 * dDdA + 4. * k3 * D5 * dkdA);
263 dsdA += (5. / 7168.) * (7. * D6 * k6 * dDdA + 6. * k5 * D7 * dkdA);
264 float ddzdA = 0.5 * dz_inv *
265 (2. * (dsdA)*dz2 * s_inv +
266 s * s * (2. * dzdl * ddzdldA * onedzdl2_inv * onedzdl2_inv));
267 dAdAp(3, 0) = -sign_dzdl * ddzdA;
275 dDdA = 0.5 * D_inv * (2. * (dx - x0) * ddxdA + 2. * (dy - y0) * ddydA);
277 dsdA += (1. / 24.) * (3. * k2 * D2 * dDdA + 2. * k * D3 * dkdA);
278 dsdA += (3. / 640.) * (5. * D4 * k4 * dDdA + 4. * k3 * D5 * dkdA);
279 dsdA += (5. / 7168.) * (7. * D6 * k6 * dDdA + 6. * k5 * D7 * dkdA);
280 ddzdA = 0.5 * dz_inv *
281 (2. * (dsdA)*dz2 * s_inv +
282 s * s * (2. * dzdl * ddzdldA * onedzdl2_inv * onedzdl2_inv));
283 dAdAp(3, 1) = -sign_dzdl * ddzdA;
291 dDdA = 0.5 * D_inv * (2. * (dx - x0) * ddxdA + 2. * (dy - y0) * ddydA);
293 dsdA += (1. / 24.) * (3. * k2 * D2 * dDdA + 2. * k * D3 * dkdA);
294 dsdA += (3. / 640.) * (5. * D4 * k4 * dDdA + 4. * k3 * D5 * dkdA);
295 dsdA += (5. / 7168.) * (7. * D6 * k6 * dDdA + 6. * k5 * D7 * dkdA);
296 ddzdA = 0.5 * dz_inv *
297 (2. * (dsdA)*dz2 * s_inv +
298 s * s * (2. * dzdl * ddzdldA * onedzdl2_inv * onedzdl2_inv));
299 dAdAp(3, 2) = -sign_dzdl * ddzdA;
307 dAdAp(0, 1) *= 2. * state.
nu;
308 dAdAp(1, 1) *= 2. * state.
nu;
309 dAdAp(3, 1) *= 2. * state.
nu;
310 dAdAp(4, 1) *= 2. * state.
nu;
315 Matrix<float, 3, 3>& dApdp, Vector3f&
p,
316 float ,
float cosphi,
float sinphi) {
325 float pT_inv = 1. /
pT;
326 float pT_inv2 = pT_inv * pT_inv;
328 float px = -sinphi *
pT;
329 float py = cosphi *
pT;
330 float pz = dzdl * pT / sqrt(1. - dzdl * dzdl);
333 dApdp(0, 0) = -py * pT_inv2;
335 dApdp(0, 1) = px * pT_inv2;
338 float pTneg3 = (pT_inv * pT_inv2);
339 dApdp(1, 0) = -0.003 *
Bfield * px * pTneg3;
341 dApdp(1, 1) = -0.003 *
Bfield * py * pTneg3;
343 float pmag = sqrt(pT * pT + pz * pz);
344 float pmag_inv = 1. / pmag;
345 float pmag_inv_3 = pmag_inv * pmag_inv * pmag_inv;
347 dApdp(2, 0) = -pz * pmag_inv_3 * px;
348 dApdp(2, 1) = -pz * pmag_inv_3 * py;
349 dApdp(2, 2) = -pz * pmag_inv_3 * pz + pmag_inv;
357 dApdp(1, 0) *= 0.5 / state.
nu;
358 dApdp(1, 1) *= 0.5 / state.
nu;
365 Vector3f
b = Vector3f::Zero(3);
368 float p_mag = sqrt(p.dot(p));
370 Vector3f p_unit = p / p_mag;
372 Vector3f axis = b.cross(p_unit);
373 float angle = acos(p_unit.dot(b));
375 axis /= sqrt(axis.dot(axis));
377 Matrix<float, 3, 3> rot_p = Matrix<float, 3, 3>::Zero(3, 3);
378 float cos_p = cos(angle);
379 float sin_p = sin(angle);
380 rot_p(0, 0) = cos_p + axis(0) * axis(0) * (1. - cos_p);
381 rot_p(0, 1) = axis(0) * axis(1) * (1. - cos_p) - axis(2) * sin_p;
382 rot_p(0, 2) = axis(0) * axis(2) * (1. - cos_p) + axis(1) * sin_p;
383 rot_p(1, 0) = axis(1) * axis(0) * (1. - cos_p) + axis(2) * sin_p;
384 rot_p(1, 1) = cos_p + axis(1) * axis(1) * (1. - cos_p);
385 rot_p(1, 2) = axis(1) * axis(2) * (1. - cos_p) - axis(0) * sin_p;
386 rot_p(2, 0) = axis(2) * axis(0) * (1. - cos_p) - axis(1) * sin_p;
387 rot_p(2, 1) = axis(2) * axis(1) * (1. - cos_p) + axis(0) * sin_p;
388 rot_p(2, 2) = cos_p + axis(2) * axis(2) * (1. - cos_p);
390 dpdb = rot_p * p_mag;
437 Matrix<float, 5, 5>&
Q) {
438 Matrix<float, 5, 3> dAdAp = Matrix<float, 5, 3>::Zero(5, 3);
450 Matrix<float, 3, 3> dApdp = Matrix<float, 3, 3>::Zero(3, 3);
451 Vector3f
p = Vector3f::Zero(3);
454 Matrix<float, 3, 3> dpdb = Matrix<float, 3, 3>::Zero(3, 3);
457 Matrix<float, 3, 2> dbdt = Matrix<float, 3, 2>::Zero(3, 2);
460 Matrix<float, 5, 2> dAdt = dAdAp * dApdp * dpdb * dbdt;
462 for (
unsigned int j = 0; j < 5; ++j) {
463 for (
unsigned int i = 0; i < 5; ++i) {
464 Q(i, j) = (dAdt(i, 0) * dAdt(j, 0) + dAdt(i, 1) * dAdt(j, 1)) * var;