CbmRoot
Loading...
Searching...
No Matches
KfTrackKalmanFilter.cxx
Go to the documentation of this file.
1/* Copyright (C) 2017-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Maksym Zyzak [committer], Valentina Akishina */
4
5#include "KfTrackKalmanFilter.h"
6
7#include "AlgoFairloggerCompat.h"
8#include "KfMeasurementU.h"
9#include "KfMeasurementXy.h"
10
11namespace cbm::algo::kf
12{
13
14 template<typename DataT>
15 void TrackKalmanFilter<DataT>::Filter1d(const kf::MeasurementU<DataT>& m)
16 {
17 DataT zeta, HCH;
18 DataT F0, F1, F2, F3, F4, F5, F6;
19 DataT K1, K2, K3, K4, K5, K6;
20
21 zeta = m.CosPhi() * fTr.X() + m.SinPhi() * fTr.Y() - m.U();
22
23 // F = CH'
24 F0 = m.CosPhi() * fTr.C00() + m.SinPhi() * fTr.C10();
25 F1 = m.CosPhi() * fTr.C10() + m.SinPhi() * fTr.C11();
26
27 HCH = (F0 * m.CosPhi() + F1 * m.SinPhi());
28
29 F2 = m.CosPhi() * fTr.C20() + m.SinPhi() * fTr.C21();
30 F3 = m.CosPhi() * fTr.C30() + m.SinPhi() * fTr.C31();
31 F4 = m.CosPhi() * fTr.C40() + m.SinPhi() * fTr.C41();
32 F5 = m.CosPhi() * fTr.C50() + m.SinPhi() * fTr.C51();
33 F6 = m.CosPhi() * fTr.C60() + m.SinPhi() * fTr.C61();
34
35 constexpr bool doProtect = !std::is_same<DataTscal, double>::value;
36
37 const DataTmask maskDoFilter = doProtect ? (HCH < m.Du2() * 16.f) : DataTmask(true);
38
39 // correction to HCH is needed for the case when sigma2 is so small
40 // with respect to HCH that it disappears due to the roundoff error
41 //
42 DataT w = m.Du2() + (doProtect ? (DataT(1.0000001) * HCH) : HCH);
43 DataT wi = kf::utils::iif(fMask, DataT(1.) / w, DataT(0.));
44
45 DataT zetawi = zeta / (kf::utils::iif(maskDoFilter, m.Du2(), DataT(0.)) + HCH);
46 zetawi = kf::utils::iif(fMask, zetawi, DataT(0.));
47
48 wi = kf::utils::iif(m.Du2() > DataT(0.), wi, DataT(0.));
49
50 fTr.ChiSq() += zeta * zeta * wi;
51 fTr.Ndf() += kf::utils::iif(fMask, m.Ndf(), DataT(0.));
52
53 K1 = F1 * wi;
54 K2 = F2 * wi;
55 K3 = F3 * wi;
56 K4 = F4 * wi;
57 K5 = F5 * wi;
58 K6 = F6 * wi;
59
60 fTr.X() -= F0 * zetawi;
61 fTr.Y() -= F1 * zetawi;
62 fTr.Tx() -= F2 * zetawi;
63 fTr.Ty() -= F3 * zetawi;
64 fTr.Qp() -= F4 * zetawi;
65 fTr.Time() -= F5 * zetawi;
66 fTr.Vi() -= F6 * zetawi;
67
68 fTr.C00() -= F0 * F0 * wi;
69
70 fTr.C10() -= K1 * F0;
71 fTr.C11() -= K1 * F1;
72
73 fTr.C20() -= K2 * F0;
74 fTr.C21() -= K2 * F1;
75 fTr.C22() -= K2 * F2;
76
77 fTr.C30() -= K3 * F0;
78 fTr.C31() -= K3 * F1;
79 fTr.C32() -= K3 * F2;
80 fTr.C33() -= K3 * F3;
81
82 fTr.C40() -= K4 * F0;
83 fTr.C41() -= K4 * F1;
84 fTr.C42() -= K4 * F2;
85 fTr.C43() -= K4 * F3;
86 fTr.C44() -= K4 * F4;
87
88 fTr.C50() -= K5 * F0;
89 fTr.C51() -= K5 * F1;
90 fTr.C52() -= K5 * F2;
91 fTr.C53() -= K5 * F3;
92 fTr.C54() -= K5 * F4;
93 fTr.C55() -= K5 * F5;
94
95 fTr.C60() -= K6 * F0;
96 fTr.C61() -= K6 * F1;
97 fTr.C62() -= K6 * F2;
98 fTr.C63() -= K6 * F3;
99 fTr.C64() -= K6 * F4;
100 fTr.C65() -= K6 * F5;
101 fTr.C66() -= K6 * F6;
102 }
103
104 template<typename DataT>
105 void TrackKalmanFilter<DataT>::FilterTime(DataT t, DataT dt2, const DataTmask& timeInfo)
106 {
107 // filter track with a time measurement
108
109 // F = CH'
110 DataT F0 = fTr.C50();
111 DataT F1 = fTr.C51();
112 DataT F2 = fTr.C52();
113 DataT F3 = fTr.C53();
114 DataT F4 = fTr.C54();
115 DataT F5 = fTr.C55();
116 DataT F6 = fTr.C65();
117
118 DataT HCH = fTr.C55();
119
120 DataTmask mask = fMask && timeInfo;
121
122 // when dt0 is much smaller than current time error,
123 // set track time exactly to the measurement value without filtering
124 // it helps to keep the initial time errors reasonably small
125 // the calculations in the covariance matrix are not affected
126
127 const DataTmask maskDoFilter = mask && (HCH < dt2 * 16.f);
128
129 DataT wi = kf::utils::iif(mask, DataT(1.) / (dt2 + DataT(1.0000001) * HCH), DataT(0.));
130 DataT zeta = kf::utils::iif(mask, fTr.Time() - t, DataT(0.));
131 DataT zetawi = kf::utils::iif(mask, zeta / (kf::utils::iif(maskDoFilter, dt2, DataT(0.)) + HCH), DataT(0.));
132
133 fTr.ChiSqTime() += kf::utils::iif(maskDoFilter, zeta * zeta * wi, DataT(0.));
134 fTr.NdfTime() += kf::utils::iif(mask, DataT(1.), DataT(0.));
135
136 DataT K1 = F1 * wi;
137 DataT K2 = F2 * wi;
138 DataT K3 = F3 * wi;
139 DataT K4 = F4 * wi;
140 DataT K5 = F5 * wi;
141 DataT K6 = F6 * wi;
142
143 fTr.X() -= F0 * zetawi;
144 fTr.Y() -= F1 * zetawi;
145 fTr.Tx() -= F2 * zetawi;
146 fTr.Ty() -= F3 * zetawi;
147 fTr.Qp() -= F4 * zetawi;
148 fTr.Time() -= F5 * zetawi;
149 fTr.Vi() -= F6 * zetawi;
150
151 fTr.C00() -= F0 * F0 * wi;
152
153 fTr.C10() -= K1 * F0;
154 fTr.C11() -= K1 * F1;
155
156 fTr.C20() -= K2 * F0;
157 fTr.C21() -= K2 * F1;
158 fTr.C22() -= K2 * F2;
159
160 fTr.C30() -= K3 * F0;
161 fTr.C31() -= K3 * F1;
162 fTr.C32() -= K3 * F2;
163 fTr.C33() -= K3 * F3;
164
165 fTr.C40() -= K4 * F0;
166 fTr.C41() -= K4 * F1;
167 fTr.C42() -= K4 * F2;
168 fTr.C43() -= K4 * F3;
169 fTr.C44() -= K4 * F4;
170
171 fTr.C50() -= K5 * F0;
172 fTr.C51() -= K5 * F1;
173 fTr.C52() -= K5 * F2;
174 fTr.C53() -= K5 * F3;
175 fTr.C54() -= K5 * F4;
176 fTr.C55() -= K5 * F5;
177
178 fTr.C60() -= K6 * F0;
179 fTr.C61() -= K6 * F1;
180 fTr.C62() -= K6 * F2;
181 fTr.C63() -= K6 * F3;
182 fTr.C64() -= K6 * F4;
183 fTr.C65() -= K6 * F5;
184 fTr.C66() -= K6 * F6;
185 }
186
187
188 template<typename DataT>
189 void TrackKalmanFilter<DataT>::FilterXY(const kf::MeasurementXy<DataT>& mxy, bool skipUnmeasuredCoordinates)
190 {
191 {
192 kf::MeasurementU<DataT> mx;
193 mx.SetCosPhi(DataT(1.));
194 mx.SetSinPhi(DataT(0.));
195 mx.SetU(mxy.X());
196 mx.SetDu2(mxy.Dx2());
197 mx.SetNdf(mxy.NdfX());
198
199 kf::MeasurementU<DataT> mu;
200 mu.SetCosPhi(-mxy.Dxy() / mxy.Dx2());
201 mu.SetSinPhi(DataT(1.));
202 mu.SetU(mu.CosPhi() * mxy.X() + mxy.Y());
203 mu.SetDu2(mxy.Dy2() - mxy.Dxy() * mxy.Dxy() / mxy.Dx2());
204 mu.SetNdf(mxy.NdfY());
205
206 auto maskOld = fMask;
207 if (skipUnmeasuredCoordinates) {
208 fMask = maskOld & (mxy.NdfX() > DataT(0.));
209 }
210 Filter1d(mx);
211 if (skipUnmeasuredCoordinates) {
212 fMask = maskOld & (mxy.NdfY() > DataT(0.));
213 }
214 Filter1d(mu);
215 fMask = maskOld;
216
217 return;
218 }
219
220 //----------------------------------------------------------------------------------------------
221 // the other way: filter 2 dimensions at once
222
223 const DataT TWO(2.);
224
225 DataT zeta0, zeta1, S00, S10, S11, si;
226 DataT F00, F10, F20, F30, F40, F50, F60;
227 DataT F01, F11, F21, F31, F41, F51, F61;
228 DataT K00, K10, K20, K30, K40, K50, K60;
229 DataT K01, K11, K21, K31, K41, K51, K61;
230
231 zeta0 = fTr.X() - mxy.X();
232 zeta1 = fTr.Y() - mxy.Y();
233
234 // F = CH'
235 F00 = fTr.C00();
236 F10 = fTr.C10();
237 F20 = fTr.C20();
238 F30 = fTr.C30();
239 F40 = fTr.C40();
240 F50 = fTr.C50();
241 F60 = fTr.C60();
242
243 F01 = fTr.C10();
244 F11 = fTr.C11();
245 F21 = fTr.C21();
246 F31 = fTr.C31();
247 F41 = fTr.C41();
248 F51 = fTr.C51();
249 F61 = fTr.C61();
250
251 S00 = F00 + mxy.Dx2();
252 S10 = F10 + mxy.Dxy();
253 S11 = F11 + mxy.Dy2();
254
255 si = 1.f / (S00 * S11 - S10 * S10);
256 DataT S00tmp = S00;
257 S00 = si * S11;
258 S10 = -si * S10;
259 S11 = si * S00tmp;
260
261 fTr.ChiSq() += zeta0 * zeta0 * S00 + 2.f * zeta0 * zeta1 * S10 + zeta1 * zeta1 * S11;
262 fTr.Ndf() += TWO;
263
264 K00 = F00 * S00 + F01 * S10;
265 K01 = F00 * S10 + F01 * S11;
266
267 K10 = F10 * S00 + F11 * S10;
268 K11 = F10 * S10 + F11 * S11;
269
270 K20 = F20 * S00 + F21 * S10;
271 K21 = F20 * S10 + F21 * S11;
272
273 K30 = F30 * S00 + F31 * S10;
274 K31 = F30 * S10 + F31 * S11;
275
276 K40 = F40 * S00 + F41 * S10;
277 K41 = F40 * S10 + F41 * S11;
278
279 K50 = F50 * S00 + F51 * S10;
280 K51 = F50 * S10 + F51 * S11;
281
282 K60 = F60 * S00 + F61 * S10;
283 K61 = F60 * S10 + F61 * S11;
284
285 fTr.X() -= K00 * zeta0 + K01 * zeta1;
286 fTr.Y() -= K10 * zeta0 + K11 * zeta1;
287 fTr.Tx() -= K20 * zeta0 + K21 * zeta1;
288 fTr.Ty() -= K30 * zeta0 + K31 * zeta1;
289 fTr.Qp() -= K40 * zeta0 + K41 * zeta1;
290 fTr.Time() -= K50 * zeta0 + K51 * zeta1;
291 fTr.Vi() -= K60 * zeta0 + K61 * zeta1;
292
293 fTr.C00() -= K00 * F00 + K01 * F01;
294
295 fTr.C10() -= K10 * F00 + K11 * F01;
296 fTr.C11() -= K10 * F10 + K11 * F11;
297
298 fTr.C20() -= K20 * F00 + K21 * F01;
299 fTr.C21() -= K20 * F10 + K21 * F11;
300 fTr.C22() -= K20 * F20 + K21 * F21;
301
302 fTr.C30() -= K30 * F00 + K31 * F01;
303 fTr.C31() -= K30 * F10 + K31 * F11;
304 fTr.C32() -= K30 * F20 + K31 * F21;
305 fTr.C33() -= K30 * F30 + K31 * F31;
306
307 fTr.C40() -= K40 * F00 + K41 * F01;
308 fTr.C41() -= K40 * F10 + K41 * F11;
309 fTr.C42() -= K40 * F20 + K41 * F21;
310 fTr.C43() -= K40 * F30 + K41 * F31;
311 fTr.C44() -= K40 * F40 + K41 * F41;
312
313 fTr.C50() -= K50 * F00 + K51 * F01;
314 fTr.C51() -= K50 * F10 + K51 * F11;
315 fTr.C52() -= K50 * F20 + K51 * F21;
316 fTr.C53() -= K50 * F30 + K51 * F31;
317 fTr.C54() -= K50 * F40 + K51 * F41;
318 fTr.C55() -= K50 * F50 + K51 * F51;
319
320 fTr.C60() -= K60 * F00 + K61 * F01;
321 fTr.C61() -= K60 * F10 + K61 * F11;
322 fTr.C62() -= K60 * F20 + K61 * F21;
323 fTr.C63() -= K60 * F30 + K61 * F31;
324 fTr.C64() -= K60 * F40 + K61 * F41;
325 fTr.C65() -= K60 * F50 + K61 * F51;
326 fTr.C66() -= K60 * F60 + K61 * F61;
327 }
328
329 template<typename DataT>
330 void TrackKalmanFilter<DataT>::FilterExtrapolatedXY(const kf::MeasurementXy<DataT>& m, DataT extrX, DataT extrY,
331 const std::array<DataT, kf::TrackParam<DataT>::kNtrackParam>& Jx,
332 const std::array<DataT, kf::TrackParam<DataT>::kNtrackParam>& Jy)
333 {
334 // add a 2-D measurenent (x,y) at some z, that differs from fTr.GetZ()
335 // extrX, extrY are extrapolated track parameters at z, Jx, Jy are derivatives of the extrapolation
336
337 // ! it is assumed that in the track covariance matrix all non-diagonal covariances are 0
338 // ! except of C10
339
340 auto& T = fTr;
341
342 //zeta0 = T.x + Jx[2]*T.Tx() + Jx[3]*T.Ty() + Jx[4]*T.qp - x;
343 //zeta1 = T.y + Jy[2]*T.Tx() + Jy[3]*T.Ty() + Jy[4]*T.qp - y;
344
345 DataT zeta0 = extrX - m.X();
346 DataT zeta1 = extrY - m.Y();
347
348 // H = 1 0 Jx[2] Jx[3] Jx[4] 0
349 // 0 1 Jy[2] Jy[3] Jy[4] 0
350
351 // F = CH'
352 DataT F00 = T.C00();
353 DataT F01 = T.C10();
354 DataT F10 = T.C10();
355 DataT F11 = T.C11();
356 DataT F20 = Jx[2] * T.C22();
357 DataT F21 = Jy[2] * T.C22();
358 DataT F30 = Jx[3] * T.C33();
359 DataT F31 = Jy[3] * T.C33();
360 DataT F40 = Jx[4] * T.C44();
361 DataT F41 = Jy[4] * T.C44();
362
363 // Jx[5,6] and Jy[5,6] are 0.
364
365 DataT S00 = m.Dx2() + F00 + Jx[2] * F20 + Jx[3] * F30 + Jx[4] * F40;
366 DataT S10 = m.Dxy() + F10 + Jy[2] * F20 + Jy[3] * F30 + Jy[4] * F40;
367 DataT S11 = m.Dy2() + F11 + Jy[2] * F21 + Jy[3] * F31 + Jy[4] * F41;
368
369 DataT si = DataT(1.) / (S00 * S11 - S10 * S10);
370
371 DataT S00tmp = S00;
372 S00 = si * S11;
373 S10 = -si * S10;
374 S11 = si * S00tmp;
375
376 T.ChiSq() += zeta0 * zeta0 * S00 + DataT(2.) * zeta0 * zeta1 * S10 + zeta1 * zeta1 * S11;
377 T.Ndf() += m.NdfX() + m.NdfY();
378
379 DataT K00 = F00 * S00 + F01 * S10;
380 DataT K01 = F00 * S10 + F01 * S11;
381 DataT K10 = F10 * S00 + F11 * S10;
382 DataT K11 = F10 * S10 + F11 * S11;
383 DataT K20 = F20 * S00 + F21 * S10;
384 DataT K21 = F20 * S10 + F21 * S11;
385 DataT K30 = F30 * S00 + F31 * S10;
386 DataT K31 = F30 * S10 + F31 * S11;
387 DataT K40 = F40 * S00 + F41 * S10;
388 DataT K41 = F40 * S10 + F41 * S11;
389
390 T.X() -= K00 * zeta0 + K01 * zeta1;
391 T.Y() -= K10 * zeta0 + K11 * zeta1;
392 T.Tx() -= K20 * zeta0 + K21 * zeta1;
393 T.Ty() -= K30 * zeta0 + K31 * zeta1;
394 T.Qp() -= K40 * zeta0 + K41 * zeta1;
395
396 T.C00() -= (K00 * F00 + K01 * F01);
397 T.C10() -= (K10 * F00 + K11 * F01);
398 T.C11() -= (K10 * F10 + K11 * F11);
399 T.C20() = -(K20 * F00 + K21 * F01);
400 T.C21() = -(K20 * F10 + K21 * F11);
401 T.C22() -= (K20 * F20 + K21 * F21);
402 T.C30() = -(K30 * F00 + K31 * F01);
403 T.C31() = -(K30 * F10 + K31 * F11);
404 T.C32() = -(K30 * F20 + K31 * F21);
405 T.C33() -= (K30 * F30 + K31 * F31);
406 T.C40() = -(K40 * F00 + K41 * F01);
407 T.C41() = -(K40 * F10 + K41 * F11);
408 T.C42() = -(K40 * F20 + K41 * F21);
409 T.C43() = -(K40 * F30 + K41 * F31);
410 T.C44() -= (K40 * F40 + K41 * F41);
411 }
412
413
414 template<typename DataT>
415 void TrackKalmanFilter<DataT>::MeasureVelocityWithQp()
416 {
417 // measure velocity using measured qp
418 // assuming particle mass == fMass;
419
420 const DataT kClightNsInv = kf::defs::SpeedOfLightInv<DataT>;
421
422 DataT zeta, HCH;
423 DataT F0, F1, F2, F3, F4, F5, F6;
424 DataT K1, K2, K3, K4, K5, K6;
425
426 //FilterVi(sqrt(DataT(1.) + fMass2 * fQp0 * fQp0) * kClightNsInv);
427 //return;
428
429 DataT vi0 = sqrt(DataT(1.) + fMass2 * fQp0 * fQp0) * kClightNsInv;
430
431 DataT h = fMass2 * fQp0 / sqrt(DataT(1.) + fMass2 * fQp0 * fQp0) * kClightNsInv;
432
433 zeta = vi0 + h * (fTr.Qp() - fQp0) - fTr.Vi();
434
435 fTr.Vi() = vi0;
436
437 // H = (0,0,0,0, h,0, -1)
438
439 // F = CH'
440
441 F0 = h * fTr.C40() - fTr.C60();
442 F1 = h * fTr.C41() - fTr.C61();
443 F2 = h * fTr.C42() - fTr.C62();
444 F3 = h * fTr.C43() - fTr.C63();
445 F4 = h * fTr.C44() - fTr.C64();
446 F5 = h * fTr.C54() - fTr.C65();
447 F6 = h * fTr.C64() - fTr.C66();
448
449 HCH = F4 * h - F6;
450
451 DataT wi = kf::utils::iif(fMask, DataT(1.) / HCH, DataT(0.));
452 DataT zetawi = kf::utils::iif(fMask, zeta / HCH, DataT(0.));
453 fTr.ChiSqTime() += kf::utils::iif(fMask, zeta * zeta * wi, DataT(0.));
454 fTr.NdfTime() += kf::utils::iif(fMask, DataT(1.), DataT(0.));
455
456 K1 = F1 * wi;
457 K2 = F2 * wi;
458 K3 = F3 * wi;
459 K4 = F4 * wi;
460 K5 = F5 * wi;
461 K6 = F6 * wi;
462
463 fTr.X() -= F0 * zetawi;
464 fTr.Y() -= F1 * zetawi;
465 fTr.Tx() -= F2 * zetawi;
466 fTr.Ty() -= F3 * zetawi;
467 fTr.Qp() -= F4 * zetawi;
468 fTr.Time() -= F5 * zetawi;
469 fTr.Vi() -= F6 * zetawi;
470
471 fTr.C00() -= F0 * F0 * wi;
472
473 fTr.C10() -= K1 * F0;
474 fTr.C11() -= K1 * F1;
475
476 fTr.C20() -= K2 * F0;
477 fTr.C21() -= K2 * F1;
478 fTr.C22() -= K2 * F2;
479
480 fTr.C30() -= K3 * F0;
481 fTr.C31() -= K3 * F1;
482 fTr.C32() -= K3 * F2;
483 fTr.C33() -= K3 * F3;
484
485 fTr.C40() -= K4 * F0;
486 fTr.C41() -= K4 * F1;
487 fTr.C42() -= K4 * F2;
488 fTr.C43() -= K4 * F3;
489 fTr.C44() -= K4 * F4;
490
491 fTr.C50() -= K5 * F0;
492 fTr.C51() -= K5 * F1;
493 fTr.C52() -= K5 * F2;
494 fTr.C53() -= K5 * F3;
495 fTr.C54() -= K5 * F4;
496 fTr.C55() -= K5 * F5;
497
498 fTr.C60() -= K6 * F0;
499 fTr.C61() -= K6 * F1;
500 fTr.C62() -= K6 * F2;
501 fTr.C63() -= K6 * F3;
502 fTr.C64() -= K6 * F4;
503 fTr.C65() -= K6 * F5;
504 fTr.C66() -= K6 * F6;
505
506 // fTr.Vi()( fTr.Vi() < DataT(TrackParamV::kClightNsInv) ) = DataT(TrackParamV::kClightNsInv);
507 }
508
509 template<typename DataT>
510 void TrackKalmanFilter<DataT>::FilterVi(DataT vi)
511 {
512 // set inverse velocity to vi
513
514 DataT zeta, HCH;
515 DataT F0, F1, F2, F3, F4, F5, F6;
516 DataT K1, K2, K3, K4, K5; //, K6;
517
518 zeta = fTr.Vi() - vi;
519
520 // H = (0,0,0,0, 0, 0, 1)
521
522 // F = CH'
523
524 F0 = fTr.C60();
525 F1 = fTr.C61();
526 F2 = fTr.C62();
527 F3 = fTr.C63();
528 F4 = fTr.C64();
529 F5 = fTr.C65();
530 F6 = fTr.C66();
531
532 HCH = F6;
533
534 DataT wi = kf::utils::iif(fMask, DataT(1.) / HCH, DataT(0.));
535 DataT zetawi = kf::utils::iif(fMask, zeta / HCH, DataT(0.));
536 fTr.ChiSqTime() += kf::utils::iif(fMask, zeta * zeta * wi, DataT(0.));
537 fTr.NdfTime() += kf::utils::iif(fMask, DataT(1.), DataT(0.));
538
539 K1 = F1 * wi;
540 K2 = F2 * wi;
541 K3 = F3 * wi;
542 K4 = F4 * wi;
543 K5 = F5 * wi;
544 // K6 = F6 * wi;
545
546 fTr.X() -= F0 * zetawi;
547 fTr.Y() -= F1 * zetawi;
548 fTr.Tx() -= F2 * zetawi;
549 fTr.Ty() -= F3 * zetawi;
550 fTr.Qp() -= F4 * zetawi;
551 fTr.Time() -= F5 * zetawi;
552 // fTr.Vi() -= F6 * zetawi;
553 fTr.Vi() = vi;
554
555 fTr.C00() -= F0 * F0 * wi;
556
557 fTr.C10() -= K1 * F0;
558 fTr.C11() -= K1 * F1;
559
560 fTr.C20() -= K2 * F0;
561 fTr.C21() -= K2 * F1;
562 fTr.C22() -= K2 * F2;
563
564 fTr.C30() -= K3 * F0;
565 fTr.C31() -= K3 * F1;
566 fTr.C32() -= K3 * F2;
567 fTr.C33() -= K3 * F3;
568
569 fTr.C40() -= K4 * F0;
570 fTr.C41() -= K4 * F1;
571 fTr.C42() -= K4 * F2;
572 fTr.C43() -= K4 * F3;
573 fTr.C44() -= K4 * F4;
574
575 fTr.C50() -= K5 * F0;
576 fTr.C51() -= K5 * F1;
577 fTr.C52() -= K5 * F2;
578 fTr.C53() -= K5 * F3;
579 fTr.C54() -= K5 * F4;
580 fTr.C55() -= K5 * F5;
581
582 //fTr.C60() -= K6 * F0;
583 //fTr.C61() -= K6 * F1;
584 //fTr.C62() -= K6 * F2;
585 //fTr.C63() -= K6 * F3;
586 //fTr.C64() -= K6 * F4;
587 //fTr.C65() -= K6 * F5;
588 //fTr.C66() -= K6 * F6;
589 fTr.C60() = DataT(0.);
590 fTr.C61() = DataT(0.);
591 fTr.C62() = DataT(0.);
592 fTr.C63() = DataT(0.);
593 fTr.C64() = DataT(0.);
594 fTr.C65() = DataT(0.);
595 fTr.C66() = DataT(1.e-8); // just for a case..
596 }
597
598 template<typename DataT>
599 void TrackKalmanFilter<
600 DataT>::Extrapolate // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
601 (DataT z_out, // extrapolate to this z position
602 const kf::FieldRegion<DataT>& F)
603 {
604 // use Q/p linearisation at fQp0
605
606 if (F.GetFieldType() == kf::EFieldType::Null) {
607 ExtrapolateLineNoField(z_out);
608 }
609 else {
610 DataT sgn = kf::utils::iif(fTr.GetZ() < z_out, DataT(1.), DataT(-1.));
611 while (
612 !kf::utils::isFull(kf::utils::iif(fMask, kf::utils::fabs(z_out - fTr.GetZ()), DataT(0.)) <= DataT(1.e-6))) {
613 DataT zNew = fTr.GetZ() + sgn * fMaxExtraplationStep; // max. 50 cm step
614 zNew = kf::utils::iif(sgn * (z_out - zNew) <= DataT(0.), z_out, zNew);
615 ExtrapolateStep(zNew, F);
616 }
617 }
618 }
619
620 template<typename DataT>
621 void TrackKalmanFilter<
622 DataT>::ExtrapolateStep // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
623 (DataT zOut, // extrapolate to this z position
624 const kf::FieldRegion<DataT>& Field)
625 {
626 // use Q/p linearisation at fQp0
627 // implementation of the Runge-Kutta method without optimization
628 //
629
630 //
631 // Forth-order Runge-Kutta method for solution of the equation
632 // of motion of a particle with parameter qp = Q /P
633 // in the magnetic field B()
634 //
635 // ( x ) tx
636 // ( y ) ty
637 // ( tx) c_light * qp * L * ( tx*ty * Bx - (1+tx*tx) * By + ty * Bz )
638 // d ( ty) / dz = c_light * qp * L * ( (1+ty*ty) * Bx - tx*ty * By - tx * Bz ) ,
639 // ( qp) 0.
640 // ( t ) L * vi
641 // ( vi) 0.
642 //
643 // where L = sqrt ( 1 + tx*tx + ty*ty ) .
644 // c_light = 0.000299792458 [(GeV/c)/kG/cm]
645 // c_light_ns = 29.9792458 [cm/ns]
646 //
647 // In the following for RK step :
648 // r[7] = {x, y, tx, ty, qp, t, vi}
649 // dr(z)/dz = f(z,r)
650 //
651 //
652 //========================================================================
653 //
654 // NIM A395 (1997) 169-184; NIM A426 (1999) 268-282.
655 //
656 // the routine is based on LHC(b) utility code
657 //
658 //========================================================================
659
660 // TODO: check the units, where does 1.e-5 factor come?
661 // c_light = 0.000299792458; // [(GeV/c)/kG/cm]
662 // c_light_ns = 29.9792458 [cm/ns]
663
664 cnst c_light = DataT(1.e-5 * kf::defs::SpeedOfLight<double>);
665
666 //----------------------------------------------------------------
667
668 cnst zMasked = kf::utils::iif(fMask, zOut, fTr.GetZ());
669
670 cnst h = (zMasked - fTr.GetZ());
671
672 cnst stepDz[5] = {0., 0., h * DataT(0.5), h * DataT(0.5), h};
673
674 DataT f[5][7] = {{DataT(0.)}}; // ( d*/dz ) [step]
675 DataT F[5][7][7] = {{{DataT(0.)}}}; // ( d *new [step] / d *old )
676
677 // Runge-Kutta steps
678 //
679
680 DataT r0[7] = {fTr.X(), fTr.Y(), fTr.Tx(), fTr.Ty(), fQp0, fTr.Time(), fTr.Vi()};
681 DataT R0[7][7] = {{DataT(0.)}};
682 for (int i = 0; i < 7; ++i) {
683 R0[i][i] = 1.;
684 }
685
686 for (int step = 1; step <= 4; ++step) {
687
688 DataT rstep[7] = {DataT(0.)};
689 for (int i = 0; i < 7; ++i) {
690 rstep[i] = r0[i] + stepDz[step] * f[step - 1][i];
691 }
692 DataT z = fTr.GetZ() + stepDz[step];
693
694 kf::FieldValue B = Field.Get(rstep[0], rstep[1], z);
695 DataT Bx = B.GetBx();
696 DataT By = B.GetBy();
697 DataT Bz = B.GetBz();
698 // NOTE: SZh 13.08.2024: The rstep[0] and rstep[1] make no effect on the B value, if Field is an approximation
699
700 DataT tx = rstep[2];
701 DataT ty = rstep[3];
702 DataT tx2 = tx * tx;
703 DataT ty2 = ty * ty;
704 DataT txty = tx * ty;
705 DataT L2 = DataT(1.) + tx2 + ty2;
706 DataT L2i = DataT(1.) / L2;
707 DataT L = sqrt(L2);
708 DataT cL = c_light * L;
709 DataT cLqp0 = cL * fQp0;
710
711 f[step][0] = tx;
712 F[step][0][2] = 1.;
713
714 f[step][1] = ty;
715 F[step][1][3] = 1.;
716
717 DataT f2tmp = txty * Bx - (DataT(1.) + tx2) * By + ty * Bz;
718 f[step][2] = cLqp0 * f2tmp;
719
720 F[step][2][2] = cLqp0 * (tx * f2tmp * L2i + ty * Bx - DataT(2.) * tx * By);
721 F[step][2][3] = cLqp0 * (ty * f2tmp * L2i + tx * Bx + Bz);
722 F[step][2][4] = cL * f2tmp;
723
724 DataT f3tmp = -txty * By - tx * Bz + (DataT(1.) + ty2) * Bx;
725 f[step][3] = cLqp0 * f3tmp;
726 F[step][3][2] = cLqp0 * (tx * f3tmp * L2i - ty * By - Bz);
727 F[step][3][3] = cLqp0 * (ty * f3tmp * L2i + DataT(2.) * ty * Bx - tx * By);
728 F[step][3][4] = cL * f3tmp;
729
730 f[step][4] = 0.;
731
732 if (fDoFitVelocity) {
733 DataT vi = rstep[6];
734 f[step][5] = vi * L;
735 F[step][5][2] = vi * tx / L;
736 F[step][5][3] = vi * ty / L;
737 F[step][5][4] = 0.;
738 F[step][5][5] = 0.;
739 F[step][5][6] = L;
740 }
741 else {
742 DataT vi = sqrt(DataT(1.) + fMass2 * fQp0 * fQp0) * kf::defs::SpeedOfLightInv<DataT>;
743 f[step][5] = vi * L;
744 F[step][5][2] = vi * tx / L;
745 F[step][5][3] = vi * ty / L;
746 F[step][5][4] = fMass2 * fQp0 * L / sqrt(DataT(1.) + fMass2 * fQp0 * fQp0) * kf::defs::SpeedOfLightInv<DataT>;
747 F[step][5][5] = 0.;
748 F[step][5][6] = 0.;
749 }
750
751 f[step][6] = 0.;
752
753 } // end of Runge-Kutta step
754
755 DataT r[7] = {DataT(0.)}; // extrapolated parameters
756 DataT R[7][7] = {{DataT(0.)}}; // Jacobian of the extrapolation
757
758 cnst stepW[5] = {0., h / DataT(6.), h / DataT(3.), h / DataT(3.), h / DataT(6.)};
759
760 DataT k[5][7][7] = {{{DataT(0.)}}};
761 for (int step = 1; step <= 4; ++step) {
762 for (int i = 0; i < 7; i++) {
763 for (int j = 0; j < 7; j++) {
764 k[step][i][j] = F[step][i][j];
765 for (int m = 0; m < 7; m++) {
766 k[step][i][j] += stepDz[step] * F[step][i][m] * k[step - 1][m][j];
767 }
768 }
769 }
770 }
771
772 for (int i = 0; i < 7; i++) {
773 for (int j = 0; j < 7; j++) {
774 R[i][j] = R0[i][j];
775 for (int step = 1; step <= 4; step++) {
776 R[i][j] += stepW[step] * k[step][i][j];
777 }
778 }
779 }
780
781 DataT dqp = fTr.Qp() - fQp0;
782
783 for (int i = 0; i < 7; i++) {
784 r[i] = r0[i];
785 for (int step = 1; step <= 4; step++) {
786 r[i] += stepW[step] * f[step][i];
787 }
788 // take into account linearisation at fQp0
789 r[i] += R[i][4] * dqp;
790 }
791
792 // update parameters
793
794 fTr.X() = r[0];
795 fTr.Y() = r[1];
796 fTr.Tx() = r[2];
797 fTr.Ty() = r[3];
798 fTr.Qp() = r[4];
799 fTr.Time() = r[5];
800 fTr.Vi() = r[6];
801
802 //fTr.Vi()( fTr.Vi() < DataT(TrackParamV::kClightNsInv) ) = DataT(TrackParamV::kClightNsInv);
803 fTr.Z() = zMasked;
804
805 // covariance matrix transport
806
807 DataT C[7][7];
808 for (int i = 0; i < 7; i++) {
809 for (int j = 0; j < 7; j++) {
810 C[i][j] = fTr.C(i, j);
811 }
812 }
813
814 DataT RC[7][7];
815 for (int i = 0; i < 7; i++) {
816 for (int j = 0; j < 7; j++) {
817 RC[i][j] = 0.;
818 for (int m = 0; m < 7; m++) {
819 RC[i][j] += R[i][m] * C[m][j];
820 }
821 }
822 }
823 for (int i = 0; i < 7; i++) {
824 for (int j = 0; j < 7; j++) {
825 DataT Cij = 0.;
826 for (int m = 0; m < 7; m++) {
827 Cij += RC[i][m] * R[j][m];
828 }
829 fTr.C(i, j) = Cij;
830 }
831 }
832 }
833
834
835 template<typename DataT>
836 void TrackKalmanFilter<DataT>::ExtrapolateLine(DataT z_out, const kf::FieldRegion<DataT>& F)
837 {
838 // extrapolate the track assuming fQp0 == 0
839 // TODO: write special simplified procedure
840 //
841 auto qp0 = fQp0;
842 fQp0 = DataT(0.);
843 Extrapolate(z_out, F);
844 fQp0 = qp0;
845 }
846
847 template<typename DataT>
848 void TrackKalmanFilter<DataT>::ExtrapolateLineNoField(DataT zOut)
849 {
850 // extrapolate the track assuming no field
851
852 // x += dz * tx
853 // y += dz * ty
854 // t += dz * sqrt ( 1 + tx*tx + ty*ty ) * vi
855
856 if (0) { // debug: full Runge-Kutta extrapolation with zero field values
857 kf::FieldRegion<DataT> F;
858 ExtrapolateStep(zOut, F);
859 return;
860 }
861
862 auto& t = fTr; // use reference to shorten the text
863
864 cnst zMasked = kf::utils::iif(fMask, zOut, t.GetZ());
865
866 cnst dz = (zMasked - t.GetZ());
867
868 DataT tx = t.GetTx();
869 DataT ty = t.GetTy();
870 DataT vi = t.GetVi();
871
872 DataT L = sqrt(DataT(1.) + tx * tx + ty * ty);
873
874 DataT j52 = dz * tx * vi / L;
875 DataT j53 = dz * ty * vi / L;
876 DataT j56 = dz * L;
877
878 // transport parameters
879
880 t.X() += tx * dz;
881 t.Y() += ty * dz;
882 t.Time() += L * vi * dz;
883 t.Z() = zMasked;
884
885 // transport covariance matrix
886
887 // ( 1 0 dz 0 0 0 0 )
888 // ( 0 1 0 dz 0 0 0 )
889 // ( 0 0 1 0 0 0 0 )
890 // J = ( 0 0 0 1 0 0 0 )
891 // ( 0 0 0 0 1 0 0 )
892 // ( 0 0 j52 j53 0 1 j56 )
893 // ( 0 0 0 0 0 0 1 )
894
895
896 // JC = J * C
897
898 DataT jc00 = t.C00() + dz * t.C20();
899 //DataT jc01 = t.C01() + dz * t.C21();
900 DataT jc02 = t.C02() + dz * t.C22();
901 //DataT jc03 = t.C03() + dz * t.C23();
902 //DataT jc04 = t.C04() + dz * t.C24();
903 //DataT jc05 = t.C05() + dz * t.C25();
904 //DataT jc06 = t.C06() + dz * t.C26();
905
906
907 DataT jc10 = t.C10() + dz * t.C30();
908 DataT jc11 = t.C11() + dz * t.C31();
909 DataT jc12 = t.C12() + dz * t.C32();
910 DataT jc13 = t.C13() + dz * t.C33();
911 //DataT jc14 = t.C14() + dz * t.C34();
912 //DataT jc15 = t.C15() + dz * t.C35();
913 //DataT jc16 = t.C16() + dz * t.C36();
914
915 // jc2? = t.C2?
916 // jc3? = t.C3?
917 // jc4? = t.C4?
918
919 DataT jc50 = t.C50() + j52 * t.C20() + j53 * t.C30() + j56 * t.C60();
920 DataT jc51 = t.C51() + j52 * t.C21() + j53 * t.C31() + j56 * t.C61();
921 DataT jc52 = t.C52() + j52 * t.C22() + j53 * t.C32() + j56 * t.C62();
922 DataT jc53 = t.C53() + j52 * t.C23() + j53 * t.C33() + j56 * t.C63();
923 DataT jc54 = t.C54() + j52 * t.C24() + j53 * t.C34() + j56 * t.C64();
924 DataT jc55 = t.C55() + j52 * t.C25() + j53 * t.C35() + j56 * t.C65();
925 DataT jc56 = t.C56() + j52 * t.C26() + j53 * t.C36() + j56 * t.C66();
926
927 // jc6? = t.C6?
928
929 // transpose J
930 //
931 // ( 1 0 0 0 0 0 0 )
932 // ( 0 1 0 0 0 0 0 )
933 // ( dz 0 1 0 0 j52 0 )
934 // J' = ( 0 dz 0 1 0 j53 0 )
935 // ( 0 0 0 0 1 0 0 )
936 // ( 0 0 0 0 0 1 0 )
937 // ( 0 0 0 0 0 j56 1 )
938
939
940 // C = JC * J'
941
942 t.C00() = jc00 + jc02 * dz;
943 t.C10() = jc10 + jc12 * dz;
944 t.C20() = t.C20() + t.C22() * dz;
945 t.C30() = t.C30() + t.C32() * dz;
946 t.C40() = t.C40() + t.C42() * dz;
947 t.C50() = jc50 + jc52 * dz;
948 t.C60() = t.C60() + t.C62() * dz;
949
950 t.C11() = jc11 + jc13 * dz;
951 t.C21() = t.C21() + t.C23() * dz;
952 t.C31() = t.C31() + t.C33() * dz;
953 t.C41() = t.C41() + t.C43() * dz;
954 t.C51() = jc51 + jc53 * dz;
955 t.C61() = t.C61() + t.C63() * dz;
956
957 // t.C22 = jc22 == t.C22 -> unchanged
958 // t.C32 = jc32 == t.C32 -> unchanged
959 // t.C42 = jc42 == t.C42 -> unchanged
960 t.C52() = jc52;
961 // t.C62 = jc62 == t.C62 -> unchanged
962
963 // t.C33 = jc33 == t.C33 -> unchanged
964 // t.C43 = jc43 == t.C43 -> unchanged
965 t.C53() = jc53;
966 // t.C63 = jc63 == t.C63 -> unchanged
967
968 // t.C44 = jc44 == t.C44 -> unchanged
969 t.C54() = jc54;
970 // t.C64 = jc64 == t.C64 -> unchanged
971
972 t.C55() = jc55 + jc52 * j52 + jc53 * j53 + jc56 * j56;
973 t.C65() = t.C65() + t.C62() * j52 + t.C63() * j53 + t.C66() * j56;
974
975 // t.C66 = jc66 = t.C66 -> unchanged
976 }
977
978
979 template<typename DataT>
980 void TrackKalmanFilter<DataT>::MultipleScattering(DataT radThick, DataT tx, DataT ty, DataT qp)
981 {
982 // 1974 - Highland (PDG) correction for multiple scattering
983 // In the formula there is a replacement:
984 // log(X/X0) == log( RadThick * sqrt(1+h) ) -> log(RadThick) + 0.5 * log(1+h);
985 // then we use an approximation:
986 // log(1+h) -> h - h^2/2 + h^3/3 - h^4/4
987
988 DataT txtx = tx * tx;
989 DataT tyty = ty * ty;
990 DataT txtx1 = DataT(1.) + txtx;
991 DataT tyty1 = DataT(1.) + tyty;
992 DataT t = sqrt(txtx1 + tyty);
993 // DataT qpt = qp * t;
994
995 DataT lg = DataT(.0136) * (DataT(1.) + DataT(0.038) * log(radThick * t));
996 lg = kf::utils::iif(lg > DataT(0.), lg, DataT(0.));
997
998 DataT s0 = lg * qp * t;
999 DataT a = (DataT(1.) + fMass2 * qp * qp) * s0 * s0 * t * radThick;
1000
1001 // Approximate formula
1002
1003 // DataT h = txtx + tyty;
1004 // DataT h2 = h * h;
1005 // cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f, c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
1006 // DataT s0 = (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
1007 // DataT a = ( (kONE+mass2*qp0*qp0t)*radThick*s0*s0 );
1008 // DataT a = ((t + fMass2 * qp0 * qp0t) * radThick * s0 * s0);
1009
1010 fTr.C22() += kf::utils::iif(fMask, txtx1 * a, DataT(0.));
1011 fTr.C32() += kf::utils::iif(fMask, tx * ty * a, DataT(0.));
1012 fTr.C33() += kf::utils::iif(fMask, tyty1 * a, DataT(0.));
1013 }
1014
1015 template<typename DataT>
1016 void TrackKalmanFilter<DataT>::MultipleScatteringInThickMaterial(DataT radThick, DataT thickness, bool fDownstream)
1017 {
1018 cnst kONE = 1.;
1019
1020 DataT tx = fTr.Tx();
1021 DataT ty = fTr.Ty();
1022 DataT txtx = tx * tx;
1023 DataT tyty = ty * ty;
1024 DataT txtx1 = txtx + kONE;
1025 DataT h = txtx + tyty;
1026 DataT t = sqrt(txtx1 + tyty);
1027 DataT h2 = h * h;
1028 DataT qp0t = fQp0 * t;
1029
1030 cnst c1(0.0136), c2 = c1 * DataT(0.038), c3 = c2 * DataT(0.5), c4 = -c3 / DataT(2.0), c5 = c3 / DataT(3.0),
1031 c6 = -c3 / DataT(4.0);
1032
1033 DataT s0 = (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
1034 //DataT a = ( (ONE+mass2*qp0*qp0t)*radThick*s0*s0 );
1035 DataT a = ((t + fMass2 * fQp0 * qp0t) * radThick * s0 * s0);
1036
1037 DataT D = (fDownstream) ? DataT(1.) : DataT(-1.);
1038 DataT T23 = (thickness * thickness) / DataT(3.0);
1039 DataT T2 = thickness / DataT(2.0);
1040
1041 fTr.C00() += kf::utils::iif(fMask, txtx1 * a * T23, DataT(0.));
1042 fTr.C10() += kf::utils::iif(fMask, tx * ty * a * T23, DataT(0.));
1043 fTr.C20() += kf::utils::iif(fMask, txtx1 * a * D * T2, DataT(0.));
1044 fTr.C30() += kf::utils::iif(fMask, tx * ty * a * D * T2, DataT(0.));
1045
1046 fTr.C11() += kf::utils::iif(fMask, (kONE + tyty) * a * T23, DataT(0.));
1047 fTr.C21() += kf::utils::iif(fMask, tx * ty * a * D * T2, DataT(0.));
1048 fTr.C31() += kf::utils::iif(fMask, (kONE + tyty) * a * D * T2, DataT(0.));
1049
1050 fTr.C22() += kf::utils::iif(fMask, txtx1 * a, DataT(0.));
1051 fTr.C32() += kf::utils::iif(fMask, tx * ty * a, DataT(0.));
1052 fTr.C33() += kf::utils::iif(fMask, (kONE + tyty) * a, DataT(0.));
1053 }
1054
1055
1056 template<typename DataT>
1057 void TrackKalmanFilter<DataT>::EnergyLossCorrection(DataT radThick, FitDirection direction)
1058 {
1059 cnst qp2cut(1. / (10. * 10.)); // 10 GeV cut
1060 cnst qp02 = kf::utils::max(fQp0 * fQp0, qp2cut);
1061 cnst p2 = DataT(1.) / qp02;
1062 cnst E2 = fMass2 + p2;
1063
1064 cnst bethe = ApproximateBetheBloch(p2 / fMass2);
1065
1066 DataT tr = sqrt(DataT(1.f) + fTr.Tx() * fTr.Tx() + fTr.Ty() * fTr.Ty());
1067
1068 DataT dE = bethe * radThick * tr * 2.33f * 9.34961f;
1069
1070 if (direction == FitDirection::kDownstream) dE = -dE;
1071
1072 cnst ECorrected = sqrt(E2) + dE;
1073 cnst E2Corrected = ECorrected * ECorrected;
1074
1075 DataT corr = sqrt(p2 / (E2Corrected - fMass2));
1076 DataTmask ok = (corr == corr) && fMask;
1077 corr = kf::utils::iif(ok, corr, DataT(1.));
1078
1079 fQp0 *= corr;
1080 fTr.Qp() *= corr;
1081 fTr.C40() *= corr;
1082 fTr.C41() *= corr;
1083 fTr.C42() *= corr;
1084 fTr.C43() *= corr;
1085 fTr.C44() *= corr * corr;
1086 fTr.C54() *= corr;
1087 }
1088
1089
1090 template<typename DataT>
1091 void TrackKalmanFilter<DataT>::EnergyLossCorrection(int atomicZ, DataTscal atomicA, DataTscal rho, DataTscal radLen,
1092 DataT radThick, FitDirection direction)
1093 {
1094 cnst qp2cut(1. / (10. * 10.)); // 10 GeV cut
1095 cnst qp02 = kf::utils::max(fQp0 * fQp0, qp2cut);
1096 cnst p2 = DataT(1.) / qp02;
1097 cnst E2 = fMass2 + p2;
1098
1099 DataT i;
1100 if (atomicZ < 13) {
1101 i = (12. * atomicZ + 7.) * 1.e-9;
1102 }
1103 else {
1104 i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
1105 }
1106
1107 cnst bethe = ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
1108
1109 DataT tr = sqrt(DataT(1.f) + fTr.Tx() * fTr.Tx() + fTr.Ty() * fTr.Ty());
1110
1111 DataT dE = bethe * radThick * tr * radLen * rho;
1112
1113 if (direction == FitDirection::kDownstream) dE = -dE;
1114
1115 cnst ECorrected = (sqrt(E2) + dE);
1116 cnst E2Corrected = ECorrected * ECorrected;
1117
1118 DataT corr = sqrt(p2 / (E2Corrected - fMass2));
1119 DataTmask ok = (corr == corr) && fMask;
1120 corr = kf::utils::iif(ok, corr, DataT(1.));
1121
1122 fQp0 *= corr;
1123 fTr.Qp() *= corr;
1124
1125 DataT P(sqrt(p2)); // GeV
1126
1127 DataT Z(atomicZ);
1128 DataT A(atomicA);
1129 DataT RHO(rho);
1130
1131 DataT STEP = radThick * tr * radLen;
1132 DataT EMASS(0.511 * 1e-3); // GeV
1133
1134 DataT BETA = P / ECorrected;
1135 DataT GAMMA = ECorrected / fMass;
1136
1137 // Calculate xi factor (KeV).
1138 DataT XI = (DataT(153.5) * Z * STEP * RHO) / (A * BETA * BETA);
1139
1140 // Maximum energy transfer to atomic electron (KeV).
1141 DataT ETA = BETA * GAMMA;
1142 DataT ETASQ = ETA * ETA;
1143 DataT RATIO = EMASS / fMass;
1144 DataT F1 = DataT(2.) * EMASS * ETASQ;
1145 DataT F2 = DataT(1.) + DataT(2.) * RATIO * GAMMA + RATIO * RATIO;
1146 DataT EMAX = DataT(1e6) * F1 / F2;
1147
1148 DataT DEDX2 = XI * EMAX * (DataT(1.) - (BETA * BETA / DataT(2.))) * DataT(1e-12);
1149
1150 DataT P2 = P * P;
1151 DataT SDEDX = (E2 * DEDX2) / (P2 * P2 * P2);
1152
1153 // T.fTr.C40() *= corr;
1154 // T.fTr.C41() *= corr;
1155 // T.fTr.C42() *= corr;
1156 // T.fTr.C43() *= corr;
1157 // T.fTr.C44() *= corr*corr;
1158 fTr.C44() += kf::utils::fabs(SDEDX);
1159 }
1160
1161
1162 template<typename DataT>
1163 void TrackKalmanFilter<DataT>::GetExtrapolatedXYline(DataT z, const kf::FieldRegion<DataT>& F, DataT& extrX,
1164 DataT& extrY,
1165 std::array<DataT, kf::TrackParam<DataT>::kNtrackParam>& Jx,
1166 std::array<DataT, kf::TrackParam<DataT>::kNtrackParam>& Jy) const
1167 {
1168 // extrapolate track assuming it is straight (qp==0)
1169 // return the extrapolated X, Y and the derivatives of the extrapolated X and Y
1170
1171 cnst c_light(0.000299792458);
1172
1173 cnst tx = fTr.GetTx();
1174 cnst ty = fTr.GetTy();
1175 DataT dz = z - fTr.GetZ();
1176
1177 DataT xx = tx * tx;
1178 DataT yy = ty * ty;
1179 DataT xy = tx * ty;
1180
1181 DataT Ay = -xx - DataT(1.);
1182 DataT Bx = yy + DataT(1.);
1183
1184 DataT ct = c_light * sqrt(DataT(1.) + xx + yy);
1185
1186
1187 DataT Sx, Sy, Sz;
1188 std::tie(Sx, Sy, Sz) = F.GetDoubleIntegrals(fTr.GetX(), fTr.GetY(), fTr.GetZ(), //
1189 fTr.GetX() + dz * tx, fTr.GetY() + dz * ty, z);
1190
1191
1192 extrX = fTr.GetX() + tx * dz;
1193 extrY = fTr.GetY() + ty * dz;
1194
1195 Jx.fill(DataT(0.)); // for a case
1196 Jy.fill(DataT(0.));
1197
1198 Jx[0] = DataT(1.);
1199 Jx[1] = DataT(0.);
1200 Jx[2] = dz;
1201 Jx[3] = DataT(0.);
1202 Jx[4] = ct * (Sx * xy + Sy * Ay + Sz * ty);
1203 Jx[5] = DataT(0.);
1204 Jx[6] = DataT(0.);
1205
1206 Jy[0] = DataT(0.);
1207 Jy[1] = DataT(1.);
1208 Jy[2] = DataT(0.);
1209 Jy[3] = dz;
1210 Jy[4] = ct * (Sx * Bx - Sy * xy - Sz * tx);
1211 Jy[5] = DataT(0.);
1212 Jx[6] = DataT(0.);
1213 }
1214
1215
1216 template<typename DataT>
1217 void TrackKalmanFilter<DataT>::FilterWithTargetAtLine(DataT targZ, const kf::MeasurementXy<DataT>& targXY,
1218 const kf::FieldRegion<DataT>& F)
1219 {
1220 // Add the target constraint to a straight line track
1221
1222 DataT eX, eY;
1223 std::array<DataT, kf::TrackParamV::kNtrackParam> Jx, Jy;
1224 GetExtrapolatedXYline(targZ, F, eX, eY, Jx, Jy);
1225 FilterExtrapolatedXY(targXY, eX, eY, Jx, Jy);
1226 }
1227
1228 template<typename DataT>
1229 DataT TrackKalmanFilter<DataT>::ApproximateBetheBloch(DataT bg2)
1230 {
1231 //
1232 // This is the parameterization of the Bethe-Bloch formula inspired by Geant.
1233 //
1234 // bg2 - (beta*gamma)^2
1235 // kp0 - density [g/cm^3]
1236 // kp1 - density effect first junction point
1237 // kp2 - density effect second junction point
1238 // kp3 - mean excitation energy [GeV]
1239 // kp4 - mean Z/A
1240 //
1241 // The default values for the kp* parameters are for silicon.
1242 // The returned value is in [GeV/(g/cm^2)].
1243 //
1244
1245 cnst kp0 = 2.33f;
1246 cnst kp1 = 0.20f;
1247 cnst kp2 = 3.00f;
1248 cnst kp3 = 173e-9f;
1249 cnst kp4 = 0.49848f;
1250
1251 constexpr DataTscal mK = 0.307075e-3f; // [GeV*cm^2/g]
1252 constexpr DataTscal _2me = 1.022e-3f; // [GeV/c^2]
1253 cnst rho = kp0;
1254 cnst x0 = kp1 * 2.303f;
1255 cnst x1 = kp2 * 2.303f;
1256 cnst mI = kp3;
1257 cnst mZA = kp4;
1258 cnst maxT = _2me * bg2; // neglecting the electron mass
1259
1260 //*** Density effect
1261 DataT d2(0.f);
1262 cnst x = 0.5f * log(bg2);
1263 cnst lhwI = log(28.816f * 1e-9f * sqrt(rho * mZA) / mI);
1264
1265 DataTmask init = x > x1;
1266 d2 = kf::utils::iif(init, lhwI + x - 0.5f, DataT(0.));
1267 cnst r = (x1 - x) / (x1 - x0);
1268 init = (x > x0) & (x1 > x);
1269 d2 = kf::utils::iif(init, lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r, d2);
1270
1271 return mK * mZA * (DataT(1.f) + bg2) / bg2
1272 * (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (DataT(1.f) + bg2) - d2);
1273 }
1274
1275 template<typename DataT>
1276 DataT TrackKalmanFilter<DataT>::ApproximateBetheBloch(DataT bg2, DataT kp0, DataT kp1, DataT kp2, DataT kp3,
1277 DataT kp4)
1278 {
1279 //
1280 // This is the parameterization of the Bethe-Bloch formula inspired by Geant.
1281 //
1282 // bg2 - (beta*gamma)^2
1283 // kp0 - density [g/cm^3]
1284 // kp1 - density effect first junction point
1285 // kp2 - density effect second junction point
1286 // kp3 - mean excitation energy [GeV]
1287 // kp4 - mean Z/A
1288 //
1289 // The default values for the kp* parameters are for silicon.
1290 // The returned value is in [GeV/(g/cm^2)].
1291 //
1292
1293 // cnst &kp0 = 2.33f;
1294 // cnst &kp1 = 0.20f;
1295 // cnst &kp2 = 3.00f;
1296 // cnst &kp3 = 173e-9f;
1297 // cnst &kp4 = 0.49848f;
1298
1299 constexpr DataTscal mK = 0.307075e-3f; // [GeV*cm^2/g]
1300 constexpr DataTscal _2me = 1.022e-3f; // [GeV/c^2]
1301 DataT rho = kp0;
1302 cnst x0 = kp1 * 2.303f;
1303 cnst x1 = kp2 * 2.303f;
1304 DataT mI = kp3;
1305 DataT mZA = kp4;
1306 cnst maxT = _2me * bg2; // neglecting the electron mass
1307
1308 //*** Density effect
1309 DataT d2(0.f);
1310 cnst x = 0.5f * log(bg2);
1311 cnst lhwI = log(28.816f * 1e-9f * sqrt(rho * mZA) / mI);
1312
1313 DataTmask init = x > x1;
1314 d2 = kf::utils::iif(init, lhwI + x - 0.5f, DataT(0.));
1315 cnst r = (x1 - x) / (x1 - x0);
1316 init = (x > x0) & (x1 > x);
1317 d2 = kf::utils::iif(init, lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r, d2);
1318
1319 return mK * mZA * (DataT(1.f) + bg2) / bg2
1320 * (0.5f * log(_2me * bg2 * maxT / (mI * mI)) - bg2 / (DataT(1.f) + bg2) - d2);
1321 }
1322
1323
1324 template<typename DataT>
1325 std::tuple<DataT, DataT> TrackKalmanFilter<DataT>::GetChi2XChi2U(kf::MeasurementXy<DataT> m, DataT x, DataT y,
1326 DataT C00, DataT C10, DataT C11)
1327 {
1328
1329 DataT chi2x{0.};
1330
1331 { // filter X measurement
1332 DataT zeta = x - m.X();
1333
1334 // F = CH'
1335 DataT F0 = C00;
1336 DataT F1 = C10;
1337
1338 DataT HCH = F0;
1339
1340 DataT wi = DataT(1.) / (m.Dx2() + HCH);
1341 DataT zetawi = zeta * wi;
1342 chi2x = m.NdfX() * zeta * zetawi;
1343
1344 DataT K1 = F1 * wi;
1345
1346 x -= F0 * zetawi;
1347 y -= F1 * zetawi;
1348
1349 C00 -= F0 * F0 * wi;
1350 C10 -= K1 * F0;
1351 C11 -= K1 * F1;
1352 }
1353
1354 DataT chi2u{0.};
1355
1356 { // filter U measurement, we need only chi2 here
1357 DataT cosPhi = -m.Dxy() / m.Dx2();
1358 DataT u = cosPhi * m.X() + m.Y();
1359 DataT du2 = m.Dy2() + cosPhi * m.Dxy();
1360
1361 DataT zeta = cosPhi * x + y - u;
1362
1363 // F = CH'
1364 DataT F0 = cosPhi * C00 + C10;
1365 DataT F1 = cosPhi * C10 + C11;
1366
1367 DataT HCH = (F0 * cosPhi + F1);
1368
1369 chi2u += m.NdfY() * zeta * zeta / (du2 + HCH);
1370 }
1371
1372 return std::tuple<DataT, DataT>(chi2x, chi2u);
1373 }
1374
1375
1376 template<typename DataT>
1377 void TrackKalmanFilter<DataT>::GuessTrack(const DataT& trackZ, const DataT hitX[], const DataT hitY[],
1378 const DataT hitZ[], const DataT hitT[], const DataT By[],
1379 const DataTmask hitW[], const DataTmask hitWtime[], int NHits)
1380 {
1381 // gives nice initial approximation for x,y,tx,ty - almost same as KF fit. qp - is shifted by 4%, resid_ual - ~3.5% (KF fit resid_ual - 1%).
1382
1383 const DataT c_light(0.000299792458), c_light_i(DataT(1.) / c_light);
1384
1385 DataT A0 = 0., A1 = 0., A2 = 0., A3 = 0., A4 = 0., A5 = 0.;
1386 DataT a0 = 0., a1 = 0., a2 = 0.;
1387 DataT b0 = 0., b1 = 0., b2 = 0.;
1388
1389 DataT time = 0.;
1390 DataTmask isTimeSet = DataTmask(false);
1391
1392 DataT prevZ = 0.;
1393 DataT sy = 0., Sy = 0.; // field integrals
1394
1395 for (int i = 0; i < NHits; i++) {
1396
1397 DataT w = kf::utils::iif(hitW[i], DataT(1.), DataT(0.));
1398
1399 DataTmask setTime = (!isTimeSet) && hitWtime[i] && hitW[i];
1400 time = kf::utils::iif(setTime, hitT[i], time);
1401 isTimeSet = isTimeSet || setTime;
1402
1403 DataT x = hitX[i];
1404 DataT y = hitY[i];
1405 DataT z = hitZ[i] - trackZ;
1406
1407 {
1408 DataT dZ = z - prevZ;
1409 Sy += w * dZ * sy + DataT(0.5) * dZ * dZ * By[i];
1410 sy += w * dZ * By[i];
1411 prevZ = kf::utils::iif(hitW[i], z, prevZ);
1412 }
1413
1414 DataT S = Sy;
1415
1416 DataT wz = w * z;
1417 DataT wS = w * S;
1418
1419 A0 += w;
1420 A1 += wz;
1421 A2 += wz * z;
1422 A3 += wS;
1423 A4 += wS * z;
1424 A5 += wS * S;
1425
1426 a0 += w * x;
1427 a1 += wz * x;
1428 a2 += wS * x;
1429
1430 b0 += w * y;
1431 b1 += wz * y;
1432 b2 += wS * y;
1433 }
1434
1435 DataT m00 = A0;
1436 DataT m01 = A1;
1437 DataT m02 = A3;
1438
1439 DataT m11 = A2;
1440 DataT m12 = A4;
1441
1442 DataT m22 = A5;
1443
1444 DataT m21 = m12;
1445
1446 // { m00 m01 m02 } ( a0 )
1447 // { m01 m11 m12 } = x * ( a1 )
1448 // { m02 m21 m22 } ( a2 )
1449
1450 { // triangulation row 0
1451 m11 = m00 * m11 - m01 * m01;
1452 m12 = m00 * m12 - m01 * m02;
1453 a1 = m00 * a1 - m01 * a0;
1454
1455 m21 = m00 * m21 - m02 * m01;
1456 m22 = m00 * m22 - m02 * m02;
1457 a2 = m00 * a2 - m02 * a0;
1458 }
1459
1460 { // triangulation step row 1
1461 m22 = m11 * m22 - m21 * m12;
1462 a2 = m11 * a2 - m21 * a1;
1463 }
1464
1465 DataT L = 0.;
1466 { // diagonalization row 2
1467 L = kf::utils::iif((kf::utils::fabs(m22) > DataT(1.e-4)), a2 / m22, DataT(0.));
1468 a1 -= L * m12;
1469 a0 -= L * m02;
1470 }
1471
1472 { // diagonalization row 1
1473 fTr.Tx() = a1 / m11;
1474 a0 -= fTr.Tx() * m01;
1475 }
1476
1477 { // diagonalization row 0
1478 fTr.X() = a0 / m00;
1479 }
1480
1481 DataT txtx1 = DataT(1.) + fTr.Tx() * fTr.Tx();
1482 L = L / txtx1;
1483 DataT L1 = L * fTr.Tx();
1484
1485 A1 = A1 + A3 * L1;
1486 A2 = A2 + (A4 + A4 + A5 * L1) * L1;
1487 b1 += b2 * L1;
1488
1489 // { A0 A1 } = x * ( b0 )
1490 // { A1 A2 } ( b1 )
1491
1492 A2 = A0 * A2 - A1 * A1;
1493 b1 = A0 * b1 - A1 * b0;
1494
1495 fTr.Ty() = b1 / A2;
1496 fTr.Y() = (b0 - A1 * fTr.Ty()) / A0;
1497
1498 fTr.Qp() = -L * c_light_i / sqrt(txtx1 + fTr.Ty() * fTr.Ty());
1499 fTr.Time() = time;
1500 fTr.Z() = trackZ;
1501 fTr.Vi() = kf::defs::SpeedOfLightInv<DataT>;
1502 fQp0 = fTr.Qp();
1503 }
1504
1505 template class TrackKalmanFilter<kf::fvec>;
1506 template class TrackKalmanFilter<float>;
1507 template class TrackKalmanFilter<double>;
1508
1509} // namespace cbm::algo::kf
KfMeasurementXy.h
Definition of the KfMeasurementXy class.
sqrt
friend fvec sqrt(const fvec &a)
Definition KfSimdPseudo.h:118
log
friend fvec log(const fvec &a)
Definition KfSimdPseudo.h:122
KfTrackKalmanFilter.h
Track fit utilities for the CA tracking based on the Kalman filter.
h
Data class with information on a STS local track.
cbm::algo::kf::FieldRegion
Magnetic field region, corresponding to a hit triplet.
Definition KfFieldRegion.h:257
cbm::algo::kf::FieldRegion::GetDoubleIntegrals
std::tuple< T, T, T > GetDoubleIntegrals(const T &x1, const T &y1, const T &z1, const T &x2, const T &y2, const T &z2) const
Gets the double integrals of the field along the track.
Definition KfFieldRegion.h:337
cbm::algo::kf::FieldValue
Magnetic flux density vector.
Definition KfFieldValue.h:31
cbm::algo::kf::FieldValue::GetBz
T GetBz() const
Gets magnetic flux density z-component [kG].
Definition KfFieldValue.h:108
cbm::algo::kf::FieldValue::GetBx
T GetBx() const
Gets magnetic flux x-component [kG].
Definition KfFieldValue.h:102
cbm::algo::kf::FieldValue::GetBy
T GetBy() const
Gets magnetic flux density y-component [kG].
Definition KfFieldValue.h:105
cbm::algo::kf::Field
Magnetic field manager class.
Definition KfField.h:188
cbm::algo::kf::MeasurementU
The class describes a 1D - measurement U in XY coordinate system.
Definition KfMeasurementU.h:35
cbm::algo::kf::MeasurementU::SetSinPhi
void SetSinPhi(DataT sinPhi)
Definition KfMeasurementU.h:64
cbm::algo::kf::MeasurementU::SetCosPhi
void SetCosPhi(DataT cosPhi)
Definition KfMeasurementU.h:63
cbm::algo::kf::MeasurementXy
The class describes a 2D - measurement (x, y) in XY coordinate system.
Definition KfMeasurementXy.h:28
cbm::algo::kf::TrackKalmanFilter::Filter1d
void Filter1d(const kf::MeasurementU< DataT > &m)
filter the track with the 1d measurement
Definition KfTrackKalmanFilter.cxx:15
cbm::algo::kf::TrackKalmanFilter::FilterWithTargetAtLine
void FilterWithTargetAtLine(DataT targZ, const kf::MeasurementXy< DataT > &targXYInfo, const kf::FieldRegion< DataT > &F)
add target measuremet to the track using linearisation at a straight line
Definition KfTrackKalmanFilter.cxx:1217
cbm::algo::kf::TrackKalmanFilter::MeasureVelocityWithQp
void MeasureVelocityWithQp()
measure the track velocity with the track Qp and the mass
Definition KfTrackKalmanFilter.cxx:415
cbm::algo::kf::TrackKalmanFilter::FilterExtrapolatedXY
void FilterExtrapolatedXY(const kf::MeasurementXy< DataT > &m, DataT extrX, DataT extrY, const std::array< DataT, kf::TrackParam< DataT >::kNtrackParam > &Jx, const std::array< DataT, kf::TrackParam< DataT >::kNtrackParam > &Jy)
Definition KfTrackKalmanFilter.cxx:330
cbm::algo::kf::TrackKalmanFilter::GetChi2XChi2U
static std::tuple< DataT, DataT > GetChi2XChi2U(kf::MeasurementXy< DataT > m, DataT x, DataT y, DataT C00, DataT C10, DataT C11)
git two chi^2 components of the track fit to measurement
Definition KfTrackKalmanFilter.cxx:1325
cbm::algo::kf::TrackKalmanFilter::FilterXY
void FilterXY(const kf::MeasurementXy< DataT > &m, bool skipUnmeasuredCoordinates=false)
filter the track with the XY measurement
Definition KfTrackKalmanFilter.cxx:189
cbm::algo::kf::TrackKalmanFilter::ExtrapolateLineNoField
void ExtrapolateLineNoField(DataT z_out)
extrapolate the track to the given Z assuming no magnetic field
Definition KfTrackKalmanFilter.cxx:848
cbm::algo::kf::TrackKalmanFilter::MultipleScatteringInThickMaterial
void MultipleScatteringInThickMaterial(DataT radThick, DataT thickness, bool fDownstream)
apply multiple scattering correction in thick material to the track
Definition KfTrackKalmanFilter.cxx:1016
cbm::algo::kf::TrackKalmanFilter::FilterVi
void FilterVi(DataT vi)
filter the inverse speed
Definition KfTrackKalmanFilter.cxx:510
cbm::algo::kf::TrackKalmanFilter::ApproximateBetheBloch
static DataT ApproximateBetheBloch(DataT bg2)
Approximate mean energy loss with Bethe-Bloch formula.
Definition KfTrackKalmanFilter.cxx:1229
cbm::algo::kf::TrackKalmanFilter::MultipleScattering
void MultipleScattering(DataT radThick, DataT tx0, DataT ty0, DataT qp0)
apply multiple scattering correction to the track with the given Qp0
Definition KfTrackKalmanFilter.cxx:980
cbm::algo::kf::TrackKalmanFilter::GetExtrapolatedXYline
void GetExtrapolatedXYline(DataT z, const kf::FieldRegion< DataT > &F, DataT &extrX, DataT &extrY, std::array< DataT, kf::TrackParam< DataT >::kNtrackParam > &Jx, std::array< DataT, kf::TrackParam< DataT >::kNtrackParam > &Jy) const
extrapolate track as a line, return the extrapolated X, Y and the Jacobians
Definition KfTrackKalmanFilter.cxx:1163
cbm::algo::kf::TrackKalmanFilter::GuessTrack
void GuessTrack(const DataT &trackZ, const DataT hitX[], const DataT hitY[], const DataT hitZ[], const DataT hitT[], const DataT By[], const DataTmask hitW[], const DataTmask hitWtime[], int NHits)
fast guess of track parameterts based on its hits
Definition KfTrackKalmanFilter.cxx:1377
cbm::algo::kf::TrackKalmanFilter::DataTscal
kf::utils::scaltype< DataT > DataTscal
Definition KfTrackKalmanFilter.h:44
cbm::algo::kf::TrackKalmanFilter::EnergyLossCorrection
void EnergyLossCorrection(DataT radThick, FitDirection direction)
Definition KfTrackKalmanFilter.cxx:1057
cbm::algo::kf::TrackKalmanFilter::FilterTime
void FilterTime(DataT t, DataT dt2, const DataTmask &m)
filter the track with the time measurement
Definition KfTrackKalmanFilter.cxx:105
cbm::algo::kf::TrackKalmanFilter::DataTmask
kf::utils::masktype< DataT > DataTmask
Definition KfTrackKalmanFilter.h:45
cbm::algo::kf::TrackKalmanFilter::ExtrapolateLine
void ExtrapolateLine(DataT z_out, const kf::FieldRegion< DataT > &F)
extrapolate the track to the given Z using linearization at the straight line
Definition KfTrackKalmanFilter.cxx:836
cbm::algo::kf::TrackParam
TrackParam classes of different types.
cbm::algo::kf::defs::SpeedOfLightInv
constexpr auto SpeedOfLightInv
Inverse speed of light [ns/cm].
Definition KfDefs.h:107
cbm::algo::kf::defs::SpeedOfLight
constexpr auto SpeedOfLight
Speed of light [cm/ns].
Definition KfDefs.h:104
cbm::algo::kf::utils::max
T max(const T &a, const T &b)
Definition KfUtils.h:56
cbm::algo::kf::utils::iif
fvec iif(const fmask &m, const fvec &t, const fvec &f)
Definition KfUtils.h:29
cbm::algo::kf::utils::fabs
fvec fabs(const fvec &v)
Definition KfUtils.h:30
cbm::algo::kf::utils::isFull
bool isFull(const T &m)
Definition KfUtils.h:45
cbm::algo::kf::EFieldType::Null
@ Null
No magnetic field.
cbm::algo::kf::sgn
fscal sgn(fscal x)
Definition KfSimdPseudo.h:234