CbmRoot
LitFiltration.h
Go to the documentation of this file.
1 
10 #ifndef LITFILTRATION_H_
11 #define LITFILTRATION_H_
12 
13 #include "LitPixelHit.h"
14 #include "LitScalPixelHit.h"
15 #include "LitStripHit.h"
16 #include "LitTrackParam.h"
17 
18 namespace lit {
19  namespace parallel {
20 
32  template<class T>
33  inline void
34  LitFiltration(LitTrackParam<T>& par, const LitPixelHit<T>& hit, T& chiSq) {
35  static const T ONE = 1., TWO = 2.;
36 
37  T dxx = hit.Dx * hit.Dx;
38  T dxy = hit.Dxy;
39  T dyy = hit.Dy * hit.Dy;
40 
41  // calculate residuals
42  T dx = hit.X - par.X;
43  T dy = hit.Y - par.Y;
44 
45  // Calculate and inverse residual covariance matrix
46  T t = ONE
47  / (dxx * dyy + dxx * par.C5 + dyy * par.C0 + par.C0 * par.C5
48  - dxy * dxy - TWO * dxy * par.C1 - par.C1 * par.C1);
49  T R00 = (dyy + par.C5) * t;
50  T R01 = -(dxy + par.C1) * t;
51  T R11 = (dxx + par.C0) * t;
52 
53  // Calculate Kalman gain matrix
54  T K00 = par.C0 * R00 + par.C1 * R01;
55  T K01 = par.C0 * R01 + par.C1 * R11;
56  T K10 = par.C1 * R00 + par.C5 * R01;
57  T K11 = par.C1 * R01 + par.C5 * R11;
58  T K20 = par.C2 * R00 + par.C6 * R01;
59  T K21 = par.C2 * R01 + par.C6 * R11;
60  T K30 = par.C3 * R00 + par.C7 * R01;
61  T K31 = par.C3 * R01 + par.C7 * R11;
62  T K40 = par.C4 * R00 + par.C8 * R01;
63  T K41 = par.C4 * R01 + par.C8 * R11;
64 
65  // Calculate filtered state vector
66  par.X += K00 * dx + K01 * dy;
67  par.Y += K10 * dx + K11 * dy;
68  par.Tx += K20 * dx + K21 * dy;
69  par.Ty += K30 * dx + K31 * dy;
70  par.Qp += K40 * dx + K41 * dy;
71 
72  // Calculate filtered covariance matrix
73  T cIn[15] = {par.C0,
74  par.C1,
75  par.C2,
76  par.C3,
77  par.C4,
78  par.C5,
79  par.C6,
80  par.C7,
81  par.C8,
82  par.C9,
83  par.C10,
84  par.C11,
85  par.C12,
86  par.C13,
87  par.C14};
88 
89  par.C0 += -K00 * cIn[0] - K01 * cIn[1];
90  par.C1 += -K00 * cIn[1] - K01 * cIn[5];
91  par.C2 += -K00 * cIn[2] - K01 * cIn[6];
92  par.C3 += -K00 * cIn[3] - K01 * cIn[7];
93  par.C4 += -K00 * cIn[4] - K01 * cIn[8];
94 
95  par.C5 += -K11 * cIn[5] - K10 * cIn[1];
96  par.C6 += -K11 * cIn[6] - K10 * cIn[2];
97  par.C7 += -K11 * cIn[7] - K10 * cIn[3];
98  par.C8 += -K11 * cIn[8] - K10 * cIn[4];
99 
100  par.C9 += -K20 * cIn[2] - K21 * cIn[6];
101  par.C10 += -K20 * cIn[3] - K21 * cIn[7];
102  par.C11 += -K20 * cIn[4] - K21 * cIn[8];
103 
104  par.C12 += -K30 * cIn[3] - K31 * cIn[7];
105  par.C13 += -K30 * cIn[4] - K31 * cIn[8];
106 
107  par.C14 += -K40 * cIn[4] - K41 * cIn[8];
108 
109  // Calculate chi-square
110  T xmx = hit.X - par.X;
111  T ymy = hit.Y - par.Y;
112  T norm = dxx * dyy - dxx * par.C5 - dyy * par.C0 + par.C0 * par.C5
113  - dxy * dxy + TWO * dxy * par.C1 - par.C1 * par.C1;
114  chiSq = ((xmx * (dyy - par.C5) - ymy * (dxy - par.C1)) * xmx
115  + (-xmx * (dxy - par.C1) + ymy * (dxx - par.C0)) * ymy)
116  / norm;
117  }
118 
130  template<class T>
131  inline void
132  LitFiltration(LitTrackParam<T>& par, const LitStripHit<T>& hit, T& chiSq) {
133  static const T ONE = 1., TWO = 2.;
134 
135  T duu = hit.Du * hit.Du;
136  T phiCosSq = hit.phiCos * hit.phiCos;
137  T phiSinSq = hit.phiSin * hit.phiSin;
138  T phi2SinCos = TWO * hit.phiCos * hit.phiSin;
139 
140  // residual
141  T r = hit.U - par.C0 * hit.phiCos - par.C1 * hit.phiSin;
142  T norm =
143  duu + par.C0 * phiCosSq + phi2SinCos * par.C1 + par.C5 * phiSinSq;
144  // myf norm = duu + cIn[0] * phiCos + cIn[5] * phiSin;
145  T R = rcp(norm);
146 
147  // Calculate Kalman gain matrix
148  T K0 = par.C0 * hit.phiCos + par.C1 * hit.phiSin;
149  T K1 = par.C1 * hit.phiCos + par.C5 * hit.phiSin;
150  T K2 = par.C2 * hit.phiCos + par.C6 * hit.phiSin;
151  T K3 = par.C3 * hit.phiCos + par.C7 * hit.phiSin;
152  T K4 = par.C4 * hit.phiCos + par.C8 * hit.phiSin;
153 
154  T KR0 = K0 * R;
155  T KR1 = K1 * R;
156  T KR2 = K2 * R;
157  T KR3 = K3 * R;
158  T KR4 = K4 * R;
159 
160  // Calculate filtered state vector
161  par.X += KR0 * r;
162  par.Y += KR1 * r;
163  par.Tx += KR2 * r;
164  par.Ty += KR3 * r;
165  par.Qp += KR4 * r;
166 
167  // Calculate filtered covariance matrix
168  par.C0 -= KR0 * K0;
169  par.C1 -= KR0 * K1;
170  par.C2 -= KR0 * K2;
171  par.C3 -= KR0 * K3;
172  par.C4 -= KR0 * K4;
173 
174  par.C5 -= KR1 * K1;
175  par.C6 -= KR1 * K2;
176  par.C7 -= KR1 * K3;
177  par.C8 -= KR1 * K4;
178 
179  par.C9 -= KR2 * K2;
180  par.C10 -= KR2 * K3;
181  par.C11 -= KR2 * K4;
182 
183  par.C12 -= KR3 * K3;
184  par.C13 -= KR3 * K4;
185 
186  par.C14 -= KR4 * K4;
187 
188  // Calculate chi-square
189  T ru = hit.U - par.X * hit.phiCos - par.Y * hit.phiSin;
190  chiSq =
191  (ru * ru)
192  / (duu - phiCosSq * par.C0 - phi2SinCos * par.C1 - phiSinSq * par.C5);
193  }
194 
203  const LitScalPixelHit& hit,
204  fscal& chiSq) {
205  LitPixelHitScal pixelHit;
206  pixelHit.X = hit.X;
207  pixelHit.Y = hit.Y;
208  pixelHit.Z = hit.Z;
209  pixelHit.Dx = hit.Dx;
210  pixelHit.Dy = hit.Dy;
211  pixelHit.Dxy = hit.Dxy;
212  LitFiltration<fscal>(par, pixelHit, chiSq);
213  }
214 
215  } // namespace parallel
216 } // namespace lit
217 #endif /* LITFILTRATION_H_ */
lit::parallel::LitTrackParam::C2
T C2
Definition: LitTrackParam.h:100
fscal
float fscal
Definition: L1/vectors/P4_F32vec4.h:250
lit::parallel::LitScalPixelHit
Base class for scalar pixel hits.
Definition: LitScalPixelHit.h:31
lit::parallel::LitTrackParam::C10
T C10
Definition: LitTrackParam.h:100
lit::parallel::LitFiltration
void LitFiltration(LitTrackParam< T > &par, const LitPixelHit< T > &hit, T &chiSq)
Function implements Kalman filter update step for pixel hit.
Definition: LitFiltration.h:34
lit::parallel::LitTrackParam::C13
T C13
Definition: LitTrackParam.h:100
lit::parallel::LitPixelHit::X
T X
Definition: LitPixelHit.h:62
lit::parallel::LitTrackParam::X
T X
Definition: LitTrackParam.h:92
lit::parallel::LitTrackParam::C8
T C8
Definition: LitTrackParam.h:100
LitPixelHit.h
Base class for pixel hits.
lit::parallel::LitTrackParam::Y
T Y
Definition: LitTrackParam.h:93
lit::parallel::LitScalPixelHit::Dxy
fscal Dxy
Definition: LitScalPixelHit.h:75
lit::parallel::LitTrackParam::C0
T C0
Definition: LitTrackParam.h:100
LitTrackParam.h
Track parameters data class.
lit::parallel::LitScalPixelHit::Z
fscal Z
Definition: LitScalPixelHit.h:78
lit::parallel::LitTrackParam::C9
T C9
Definition: LitTrackParam.h:100
lit::parallel::LitPixelHit::Dx
T Dx
Definition: LitPixelHit.h:63
lit::parallel::LitTrackParam::C3
T C3
Definition: LitTrackParam.h:100
LitStripHit::Du
T Du
Definition: LitStripHit.h:64
lit::parallel::LitTrackParam::C14
T C14
Definition: LitTrackParam.h:100
lit::parallel::LitTrackParam::C4
T C4
Definition: LitTrackParam.h:100
lit::parallel::LitTrackParam
Track parameters data class.
Definition: LitTrackParam.h:34
lit::parallel::LitTrackParam::C12
T C12
Definition: LitTrackParam.h:100
lit::parallel::LitTrackParam::Ty
T Ty
Definition: LitTrackParam.h:96
lit::parallel::LitScalPixelHit::Dx
fscal Dx
Definition: LitScalPixelHit.h:74
lit::parallel::LitTrackParam::Tx
T Tx
Definition: LitTrackParam.h:95
lit::parallel::LitScalPixelHit::Y
fscal Y
Definition: LitScalPixelHit.h:73
lit::parallel::LitScalPixelHit::Dy
fscal Dy
Definition: LitScalPixelHit.h:74
lit::parallel::LitScalPixelHit::X
fscal X
Definition: LitScalPixelHit.h:73
lit::parallel::LitPixelHit::Dy
T Dy
Definition: LitPixelHit.h:63
lit::parallel::LitPixelHit::Y
T Y
Definition: LitPixelHit.h:62
LitScalPixelHit.h
Base class for scalar pixel hits.
LitStripHit.h
Base class for strip hits.
LitStripHit
Base class for strip hits.
Definition: LitStripHit.h:28
lit::parallel::LitPixelHit::Dxy
T Dxy
Definition: LitPixelHit.h:64
lit::parallel::LitTrackParam::C5
T C5
Definition: LitTrackParam.h:100
lit::parallel::LitTrackParam::C1
T C1
Definition: LitTrackParam.h:100
LitStripHit::phiSin
T phiSin
Definition: LitStripHit.h:62
lit::parallel::LitPixelHit::Z
T Z
Definition: LitPixelHit.h:62
lit::parallel::LitTrackParam::C11
T C11
Definition: LitTrackParam.h:100
lit::parallel::LitTrackParam::C7
T C7
Definition: LitTrackParam.h:100
LitStripHit::phiCos
T phiCos
Definition: LitStripHit.h:61
lit::parallel::LitTrackParam::C6
T C6
Definition: LitTrackParam.h:100
lit::parallel::LitTrackParam::Qp
T Qp
Definition: LitTrackParam.h:97
lit::parallel::LitPixelHit
Base class for pixel hits.
Definition: LitPixelHit.h:29
LitStripHit::U
T U
Definition: LitStripHit.h:63
lit::parallel::rcp
fscal rcp(const fscal &a)
Returns reciprocal.
Definition: LitMath.h:26
NS_L1TrackFitter::ONE
const fvec ONE
Definition: L1TrackFitter.cxx:33
lit
Definition: LitTrackFinderNNVecElectron.h:19