FEDRA emulsion software from the OPERA Collaboration
List of all members
VERTEX::Kalman Class Reference

#include <VtKalman.hh>

Collaboration diagram for VERTEX::Kalman:
Collaboration graph
[legend]

Public Member Functions

— Constructors —
 Kalman (const Relation *const relation)
 
 ~Kalman ()
 
— Access methods —
double chi2 () const
 $\chi^2$ contribution of track More...
 
double chi2s () const
 smoothed $\chi^2$ More...
 
double tx () const
 refitted $t_x$ More...
 
double ty () const
 refitted $t_y$ More...
 
double p () const
 refitted $p$ More...
 
double px () const
 refitted $p_x$ More...
 
double py () const
 refitted $p_y$ More...
 
double pz () const
 refitted $p_z$ More...
 
double ex () const
 refitted $e_x$ More...
 
double ey () const
 refitted $e_y$ More...
 
double ez () const
 refitted $e_z$ More...
 
double E (double rm=0.) const
 refitted Energy $E = \sqrt{p^2 + m^2}$ More...
 
MATRIX::VtVector evec () const
 $\vec{v} = (e_x,e_y,e_z)$ unit vector along refitted track More...
 
MATRIX::VtVector tvec () const
 $\vec{v} = (t_x,t_y,1.)$ refitted slope vector More...
 
MATRIX::VtVector pvec () const
 $\vec{v} = (p_x,p_y,p_z)$ refitted mom. vector More...
 
— Cov. matrices for internal use —
const MATRIX::VtSymMatrix & W () const
 
const MATRIX::VtSymMatrix & C () const
 
const MATRIX::VtSymMatrix & CINV () const
 
const MATRIX::VtSqMatrix & F () const
 
const MATRIX::VtSqMatrix & ES () const
 
const MATRIX::VtSymMatrix & DS () const
 
const MATRIX::VtVector & xv () const
 
const MATRIX::VtVector & xnk () const
 
const MATRIX::VtVector & qvs () const
 $\vec{v} = (t_x,t_y,p)$ More...
 
— Expert functions —
void use_momentum (const bool use)
 set if track momentum should be used in fit More...
 
bool use_momentum () const
 
double set_chi2 (const double chi2)
 
void init ()
 initialize the kalman structure More...
 
std::ostream & print (std::ostream &os) const
 called by operator<<() More...
 
bool filter (double z, const MATRIX::VtSymMatrix &prCINV, const MATRIX::VtVector &prkal_xv)
 
double filter_chi2 (double z, double prChi2, const MATRIX::VtSymMatrix &prCINV, const MATRIX::VtVector &prkal_xv)
 
bool inverse_filter (double z, const MATRIX::VtSymMatrix &CINVn, const MATRIX::VtSymMatrix &prCINV, const MATRIX::VtVector &kal_xvn)
 
void smooth (double z, const MATRIX::VtVector &xvs, const MATRIX::VtSymMatrix &Cn)
 
const MATRIX::VtVector calc_dp (double z, const MATRIX::VtVector &xk, const MATRIX::VtVector &qk) const
 
const MATRIX::VtVector calc_pcAx (const MATRIX::VtVector &xk) const
 
const MATRIX::VtVector calc_AGpc (void) const
 
const MATRIX::VtVector calc_qk (const MATRIX::VtVector &xk) const
 
void calc_qvs (const MATRIX::VtVector &xvs)
 
double calc_dchi2 (double z, const MATRIX::VtSymMatrix &prCINV, const MATRIX::VtVector &xk, const MATRIX::VtVector &prxk, const MATRIX::VtVector &qk) const
 
void calc_pc (double z)
 

— Mass constraint fit functions —

const Relation * rel
 
MATRIX::CMatrix G
 
bool k_use_momentum
 
double k_tx
 
double k_ty
 
MATRIX::VtVector k_qvs
 
MATRIX::VtVector k_xv
 
MATRIX::VtVector k_qv
 
MATRIX::VtVector k_pc
 
MATRIX::VtVector k_xnk
 
MATRIX::VtSymMatrix k_W
 
MATRIX::VtMatrix k_GB
 
MATRIX::VtMatrix k_WBG
 
MATRIX::CMatrix k_Gb
 
MATRIX::VtSymMatrix k_C
 
MATRIX::VtSymMatrix k_CINV
 
MATRIX::VtSqMatrix k_F
 
MATRIX::VtSqMatrix k_ES
 
MATRIX::VtSymMatrix k_DS
 
MATRIX::VtVector k_alpc
 
MATRIX::VtVector k_alp
 
MATRIX::VtVector k_nalpc
 
double k_erg
 
double k_chi2
 
double k_chi2s
 
void alpc_init (void)
 fill $\vec{\alpha}_c^{(0)}$ More...
 
void alp_init (void)
 fill $\vec{\alpha}^{(0)}$ More...
 
void calc_ealpc (void)
 calculate unit vector $\vec{\alpha}_c/|\vec{\alpha}_c|$, energy from alpc More...
 
const MATRIX::VtVector & alpc () const
 state vector $\vec{\alpha_c}=(t_x,t_y,p)$ More...
 
const MATRIX::VtVector & alp () const
 state vector $\vec{\alpha}=(t_x,t_y,p)$ More...
 
double xn () const
 $n_x = t_x * e_z$ More...
 
double yn () const
 $n_y = t_y * e_z$ More...
 
double zn () const
 $n_z = 1/\sqrt{1 + t_x^2 + t_y^2}$ More...
 
double erg () const
 $E_i=\sqrt{m^2 + p^2}$, $m$ = track rest-mass More...
 
const MATRIX::VtVector & nalpc () const
 $\vec{n} = (n_x,n_y,n_z)$ More...
 
MATRIX::VtVector & qvs_nc ()
 return non-const reference to qvs More...
 

Detailed Description

This class contains the Vt Kalman filter matrices and vectors. Most of the member functions are for internal use only. @memo Kalman filter class

Constructor & Destructor Documentation

◆ Kalman()

VERTEX::Kalman::Kalman ( const Relation *const  relation)
40 :
41 rel(relation),
42 G(rel->track.G()),
43 k_use_momentum(true),
44 k_tx (0.),
45 k_ty (0.),
46 k_qvs (MATRIX::VtVector(3)),
47 k_xv (MATRIX::VtVector(3)),
48 k_qv (MATRIX::VtVector(3)),
49 k_pc (MATRIX::VtVector(5)),
50 k_xnk (MATRIX::VtVector(3)),
51 k_W (MATRIX::VtSymMatrix(3)),
52 k_GB (MATRIX::VtMatrix(5,3)),
53 k_WBG (MATRIX::VtMatrix(3,5)),
54 k_C (MATRIX::VtSymMatrix(3)),
55 k_CINV (MATRIX::VtSymMatrix(3)),
56 k_F (MATRIX::VtSqMatrix(3)),
57 k_ES (MATRIX::VtSqMatrix(3)),
58 k_DS (MATRIX::VtSymMatrix(3)),
59 k_alpc (MATRIX::VtVector(3)),
60 k_alp (MATRIX::VtVector(3)),
61 k_nalpc (MATRIX::VtVector(3)),
62 k_erg (0.),
63 k_chi2 (0.),
64 k_chi2s (0.)
65 {}
MATRIX::VtMatrix
Definition: VtMatrix.hh:49
MATRIX::VtSqMatrix
Definition: VtSqMatrix.hh:50
MATRIX::VtSymMatrix
Definition: VtSymMatrix.hh:49
MATRIX::VtVector
Definition: VtVector.hh:45
VERTEX::Kalman::rel
const Relation * rel
Definition: VtKalman.hh:218
VERTEX::Kalman::k_GB
MATRIX::VtMatrix k_GB
Definition: VtKalman.hh:232
VERTEX::Kalman::G
MATRIX::CMatrix G
Definition: VtKalman.hh:219
VERTEX::Kalman::k_alpc
MATRIX::VtVector k_alpc
Definition: VtKalman.hh:240
VERTEX::Kalman::k_xv
MATRIX::VtVector k_xv
Definition: VtKalman.hh:224
VERTEX::Kalman::k_nalpc
MATRIX::VtVector k_nalpc
Definition: VtKalman.hh:242
VERTEX::Kalman::k_chi2s
double k_chi2s
Definition: VtKalman.hh:245
VERTEX::Kalman::k_alp
MATRIX::VtVector k_alp
Definition: VtKalman.hh:241
VERTEX::Kalman::k_erg
double k_erg
Definition: VtKalman.hh:243
VERTEX::Kalman::k_chi2
double k_chi2
Definition: VtKalman.hh:244
VERTEX::Kalman::k_pc
MATRIX::VtVector k_pc
Definition: VtKalman.hh:229
VERTEX::Kalman::k_DS
MATRIX::VtSymMatrix k_DS
Definition: VtKalman.hh:239
VERTEX::Kalman::k_tx
double k_tx
Definition: VtKalman.hh:221
VERTEX::Kalman::k_WBG
MATRIX::VtMatrix k_WBG
Definition: VtKalman.hh:233
VERTEX::Kalman::k_qv
MATRIX::VtVector k_qv
Definition: VtKalman.hh:228
VERTEX::Kalman::k_ES
MATRIX::VtSqMatrix k_ES
Definition: VtKalman.hh:238
VERTEX::Kalman::k_xnk
MATRIX::VtVector k_xnk
Definition: VtKalman.hh:230
VERTEX::Kalman::k_ty
double k_ty
Definition: VtKalman.hh:222
VERTEX::Kalman::k_qvs
MATRIX::VtVector k_qvs
Definition: VtKalman.hh:223
VERTEX::Kalman::k_F
MATRIX::VtSqMatrix k_F
Definition: VtKalman.hh:237
VERTEX::Kalman::k_CINV
MATRIX::VtSymMatrix k_CINV
Definition: VtKalman.hh:236
VERTEX::Kalman::k_W
MATRIX::VtSymMatrix k_W
Definition: VtKalman.hh:231
VERTEX::Kalman::k_use_momentum
bool k_use_momentum
Definition: VtKalman.hh:220
VERTEX::Kalman::k_C
MATRIX::VtSymMatrix k_C
Definition: VtKalman.hh:235
VERTEX::Relation::track
Track & track
Definition: VtRelation.hh:72
VERTEX::Track::G
const MATRIX::CMatrix & G() const
inverse $5\times5$ covariance matrix
Definition: VtTrack.C:206

◆ ~Kalman()

VERTEX::Kalman::~Kalman ( )
67{}

Member Function Documentation

◆ alp()

const MATRIX::VtVector & VERTEX::Kalman::alp ( ) const

state vector $\vec{\alpha}=(t_x,t_y,p)$

◆ alp_init()

void VERTEX::Kalman::alp_init ( void  )

fill $\vec{\alpha}^{(0)}$

◆ alpc()

const MATRIX::VtVector & VERTEX::Kalman::alpc ( ) const

state vector $\vec{\alpha_c}=(t_x,t_y,p)$

◆ alpc_init()

void VERTEX::Kalman::alpc_init ( void  )

fill $\vec{\alpha}_c^{(0)}$

◆ C()

const MATRIX::VtSymMatrix & VERTEX::Kalman::C ( ) const

◆ calc_AGpc()

const VtVector VERTEX::Kalman::calc_AGpc ( void  ) const
111 {
112 VtVector AGpc(3);
113
114 // row 1 + 2
115 for(unsigned int i=0; i<2; ++i) {
116 for(unsigned int j=0; j < k_Gb.ncol(); ++j) {
117 AGpc[i] += k_Gb(i,j) * k_pc[j];
118 }
119 }
120 // row 3
121 for(unsigned int j=0; j < k_Gb.ncol(); ++j) {
122 AGpc[2] += (-k_tx * k_Gb(0,j) - k_ty * k_Gb(1,j)) * k_pc[j];
123 }
124
125 return AGpc;
126 }
MATRIX::VtMatrix::ncol
unsigned int ncol() const
no of columns $m$
Definition: VtMatrix.hh:63
VERTEX::Kalman::k_Gb
MATRIX::CMatrix k_Gb
Definition: VtKalman.hh:234

◆ calc_dchi2()

double VERTEX::Kalman::calc_dchi2 ( double  z,
const MATRIX::VtSymMatrix &  prCINV,
const MATRIX::VtVector &  xk,
const MATRIX::VtVector &  prxk,
const MATRIX::VtVector &  qk 
) const
240 {
241
242 // compute vector x_k - x_k-1
243 const VtVector dx = xk - prxk;
244
245#ifdef VtDEBUG
246 cout << "Vector prxk: " << prxk << endl;
247 cout << "Vector dx: " << dx << endl;
248#endif
249
250 // compute vector p_k - c_k^(0) - A_k*x_k - B_k*q_k (eq. 14)
251 VtVector dp = calc_dp(z, xk, qk);
252
253 unsigned int dim = 4;
254 if(k_use_momentum) {
255 dim = 5;
256 }
257
258 // calculate dp^T * G_k * dp and add it to chi2 (in Vt: chi2f)
259 double chi2 = G.product(dp, dim);
260
261#ifdef VtDEBUG
262 cout << "chi2 1st term: " << chi2 << endl;
263#endif
264
265 // compute term (x_k-x_k-1)^T * C_k-1^-1 * (x_k-x_k-1)
266 double chi2b = prCINV.product(dx);
267
268#ifdef VtDEBUG
269 cout << "chi2 2nd term: " << chi2b << endl;
270#endif
271
272 return chi2 + chi2b;
273 } // calc_dchi2
MATRIX::VtSymMatrix::product
double product(const VtVector &rhs, unsigned int dim=0) const
compute $v^t*A*v$
Definition: VtSymMatrix.C:351
VERTEX::Kalman::chi2
double chi2() const
$\chi^2$ contribution of track
VERTEX::Kalman::calc_dp
const MATRIX::VtVector calc_dp(double z, const MATRIX::VtVector &xk, const MATRIX::VtVector &qk) const
Definition: VtKalman.C:131

◆ calc_dp()

const VtVector VERTEX::Kalman::calc_dp ( double  z,
const MATRIX::VtVector &  xk,
const MATRIX::VtVector &  qk 
) const
133 {
134
135 // compute p_k - c_k^(0) - A_k*x_k - B_k * q_k (in Vt: dp)
136 VtVector dp = calc_pcAx(xk);
137 dp[0] += z * qk[0];
138 dp[1] += z * qk[1];
139 dp[2] -= qk[0];
140 dp[3] -= qk[1];
141
142 if(k_use_momentum) {
143 dp[4] -= qk[2];
144 }
145
146#ifdef VtDEBUG
147 cout << "Vector dp: " << dp << endl;
148 cout << " z: " << z << endl;
149 cout << "Vector qk: " << qk << endl;
150#endif
151
152 return dp;
153 }
VERTEX::Kalman::calc_pcAx
const MATRIX::VtVector calc_pcAx(const MATRIX::VtVector &xk) const
Definition: VtKalman.C:158

◆ calc_ealpc()

void VERTEX::Kalman::calc_ealpc ( void  )

calculate unit vector $\vec{\alpha}_c/|\vec{\alpha}_c|$, energy from alpc

574 {
575
576 k_nalpc[2] = 1./sqrt(1. + sqr(k_alpc[0]) + sqr(k_alpc[1]));
577 k_nalpc[0] = k_alpc[0] * k_nalpc[2];
578 k_nalpc[1] = k_alpc[1] * k_nalpc[2];
579
580 k_erg = sqrt(sqr(rel->track.rm()) + sqr(k_alpc[2]));
581 return;
582 }
VERTEX::Track::rm
double rm() const
rest-mass
VERTEX::sqr
const T sqr(const T &x)
compute the square of a number: $x*x$
Definition: VtUtil.hh:37

◆ calc_pc()

void VERTEX::Kalman::calc_pc ( double  z)
213 {
214
215 const Track& t = rel->track;
216
217 // compute vector p_k - c_k^(0) (in Vt: pc)
218 k_pc[0] = t.x() - z * k_tx;
219 k_pc[1] = t.y() - z * k_ty;
220 k_pc[2] = t.tx();
221 k_pc[3] = t.ty();
222
223 if(k_use_momentum)
224 k_pc[4] = t.p();
225
226#ifdef VtDEBUG
227 cout << "Vector k_pc: " << k_pc << endl;
228#endif
229
230 return;
231 }
Track
Definition: Track.h:10
t
TTree * t
Definition: check_shower.C:4

◆ calc_pcAx()

const VtVector VERTEX::Kalman::calc_pcAx ( const MATRIX::VtVector &  xk) const
158 {
159
160 VtVector pcAx(5);
161 pcAx[0] = k_pc[0] - xk[0] + k_tx * xk[2];
162 pcAx[1] = k_pc[1] - xk[1] + k_ty * xk[2];
163 pcAx[2] = k_pc[2];
164 pcAx[3] = k_pc[3];
165
166 if(k_use_momentum)
167 pcAx[4] = k_pc[4];
168
169#ifdef VtDEBUG
170 cout << " xvs[1]: " << xk[0]
171 << " xvs[2]: " << xk[1]
172 << " xvs[3]: " << xk[2]
173 << " a33: " << k_tx
174 << " a34: " << k_ty
175 << endl;
176 cout << "Vector pc: " << k_pc << endl;
177 cout << "Vector pcAx: " << pcAx << endl;
178#endif
179 return pcAx;
180 }

◆ calc_qk()

const VtVector VERTEX::Kalman::calc_qk ( const MATRIX::VtVector &  xk) const
185 {
186
187 // compute p_k-c_k^(0) - A_k*x_k
188 VtVector pcAx = calc_pcAx(xk);
189
190 // compute q_k^n = W_k*B_k^T*G_k * (p_k - c_k^(0) - A_k*x_n) (in Vt: qvs)
191 const VtVector qk = k_WBG * pcAx;
192
193#ifdef VtDEBUG
194 cout << "Vector qk: " << qk << endl;
195#endif
196
197 return qk;
198 }

◆ calc_qvs()

void VERTEX::Kalman::calc_qvs ( const MATRIX::VtVector &  xvs)
203 {
204
205 // compute q_k^n = W_k*B_k^T*G_k * (p_k - c_k^(0) - A_k*x_n) (in Vt: qvs)
206 k_qvs = calc_qk(xvs);
207 return;
208 }
VERTEX::Kalman::calc_qk
const MATRIX::VtVector calc_qk(const MATRIX::VtVector &xk) const
Definition: VtKalman.C:185

◆ chi2()

double VERTEX::Kalman::chi2 ( ) const

$\chi^2$ contribution of track

◆ chi2s()

double VERTEX::Kalman::chi2s ( ) const

smoothed $\chi^2$

◆ CINV()

const MATRIX::VtSymMatrix & VERTEX::Kalman::CINV ( ) const

◆ DS()

const MATRIX::VtSymMatrix & VERTEX::Kalman::DS ( ) const

◆ E()

double VERTEX::Kalman::E ( double  rm = 0.) const

refitted Energy $E = \sqrt{p^2 + m^2}$

◆ erg()

double VERTEX::Kalman::erg ( ) const

$E_i=\sqrt{m^2 + p^2}$, $m$ = track rest-mass

◆ ES()

const MATRIX::VtSqMatrix & VERTEX::Kalman::ES ( ) const

◆ evec()

MATRIX::VtVector VERTEX::Kalman::evec ( ) const

$\vec{v} = (e_x,e_y,e_z)$ unit vector along refitted track

◆ ex()

double VERTEX::Kalman::ex ( ) const

refitted $e_x$

◆ ey()

double VERTEX::Kalman::ey ( ) const

refitted $e_y$

◆ ez()

double VERTEX::Kalman::ez ( ) const

refitted $e_z$

◆ F()

const MATRIX::VtSqMatrix & VERTEX::Kalman::F ( ) const

◆ filter()

bool VERTEX::Kalman::filter ( double  z,
const MATRIX::VtSymMatrix &  prCINV,
const MATRIX::VtVector &  prkal_xv 
)
280 {
281 // task: propagate inverse cov. Matrix to estimated z of Vertex
282 // and project it into (tx,ty,p) space. Finally invert the matrix.
283 // z = -bk13 in Vt package
284
285 const double z2 = z*z;
286#ifdef VtDEBUG
287 cout << " z: " << z << endl;
288#endif
289
290 // save tx,ty at begin of filter step (needed in case
291 // of refit by track removal from vertex)
292 k_tx = k_qvs[0];
293 k_ty = k_qvs[1];
294
295#ifdef VtDEBUG
296 cout << "k_tx: " << k_tx << " k_ty: " << k_ty << endl;
297#endif
298
299 // compute Matrix W = (B^T*G*B)^-1
300 k_W(0,0) = z2 * G(0,0) - 2. * z * G(0,2) + G(2,2);
301 k_W(1,1) = z2 * G(1,1) - 2. * z * G(1,3) + G(3,3);
302 k_W(0,1) =
303 k_W(1,0) = z2 * G(1,0) - z * ( G(1,2) + G(3,0) ) + G(3,2);
304
305 if(k_use_momentum) {
306 k_W(0,2) = k_W(2,0) = G(4,2) - z * G(4,0);
307 k_W(1,2) = k_W(2,1) = G(4,3) - z * G(4,1);
308 k_W(2,2) = G(4,4);
309 }
310
311#ifdef VtDEBUG
312 cout << "Matrix G: " << G << endl;
313 cout << "Matrix k_W: " << k_W << endl;
314#endif
315
316 if(k_W.invert(k_use_momentum) == false)
317 return false;
318
319#ifdef VtDEBUG
320 cout << "Matrix inverted k_W: " << k_W << endl;
321#endif
322
323 // compute Matrix G_k*B_k = B_k^T*G_k (in Vt: AGB)
324 for(unsigned int i=0; i < G.nrow(); ++i) {
325 k_GB(i,0) = G(i,2) - z * G(i,0);
326 k_GB(i,1) = G(i,3) - z * G(i,1);
327 }
328 if(k_use_momentum) {
329 for(unsigned int i=0; i < G.nrow(); ++i)
330 k_GB(i,2) = G(i,4);
331 }
332
333#ifdef VtDEBUG
334 cout << "Matrix k_GB: " << k_GB << endl;
335#endif
336
337 // compute Matrix W*(GB)^T = W*B^T*G (in Vt: WBG)
338 k_WBG = k_W * k_GB.T();
339
340#ifdef VtDEBUG
341 cout << "Matrix k_WBG: " << k_WBG << endl;
342#endif
343
344 // compute Matrix (G*B) * (W*B^T*G) (in Vt: GBWBG)
345 const VtSymMatrix GBWBG = k_GB * k_WBG;
346
347#ifdef VtDEBUG
348 cout << "Matrix GBWBG: " << GBWBG << endl;
349#endif
350
351 // compute Matrix G^B = G - GBWB^TG (in Vt: GB, eq. 10)
352 k_Gb = G - GBWBG;
353
354#ifdef VtDEBUG
355 cout << "Matrix G: " << G;
356 cout << "Matrix k_Gb: " << k_Gb;
357#endif
358
359 // compute Matrix A*G^B*A (in Vt: AGA)
360 VtSymMatrix AGA(3);
361 AGA(0,0) = k_Gb(0,0);
362 AGA(0,1) = AGA(1,0) = k_Gb(0,1);
363 AGA(1,1) = k_Gb(1,1);
364 AGA(0,2) = AGA(2,0) = -k_tx*k_Gb(0,0) - k_ty*k_Gb(1,0);
365 AGA(1,2) = AGA(2,1) = -k_tx*k_Gb(0,1) - k_ty*k_Gb(1,1);
366 AGA(2,2) = k_tx*k_tx*k_Gb(0,0) + 2.*k_tx*k_ty*k_Gb(0,1) +
367 k_ty*k_ty*k_Gb(1,1);
368
369#ifdef VtDEBUG
370 cout << "Matrix AGA: " << AGA << endl;
371 cout << "Matrix prCINV: " << prCINV << endl;
372#endif
373
374 // compute k_C, in Vt: C_k = (C_k-1^-1 + A^T_k G^B_k A_k)^-1 (eq. 10)
375 k_CINV = prCINV + AGA;
376
377#ifdef VtDEBUG
378 cout << "Matrix k_CINV: " << k_CINV << endl;
379#endif
380
381 k_C = k_CINV;
382 if( k_C.VtDsinv() == false )
383 return false; // track does not fit onto vertex
384
385#ifdef VtDEBUG
386 cout << "Matrix k_C: " << k_C << endl;
387 cout << "k_C*k_CINV: " << k_C*k_CINV << endl;
388#endif
389
390 // --------> Determine filtered state vector <-----------
391
392 // compute equation 12 (see Vt docu)
393 // compute vector p_k - c_k^(0) (in Vt: pc)
394 calc_pc(z);
395
396 // compute Vector (A^T*G^B) * pc (in Vt: AGpc)
397 const VtVector AGpc = calc_AGpc();
398
399#ifdef VtDEBUG
400 cout << "Vector AGpc: " << AGpc << endl;
401#endif
402
403
404 // compute x_k = C_k * [C_k-1^-1*x_k-1 + A^TG^B (p_k - c_k^(0))] (eq. 12)
405 // Cx = expression in brackets [...]
406 const VtVector Cx = AGpc + prCINV * prkal_xv;
407
408#ifdef VtDEBUG
409 cout << "prkal_xv: " << prkal_xv << endl;
410 cout << "Vector Cx: " << Cx << endl;
411#endif
412
413 // compute k_xv (in Vt: x_k): filtered vtx position
414 k_xv = k_C * Cx;
415
416#ifdef VtDEBUG
417 cout << "Vector k_xv: " << k_xv << endl;
418#endif
419
420 // compute equation 13
421 // compute q_k = W_kB^T_kG_k * [p_k - c_k^(0) - A_k*x_k] (in Vt: qv)
422 k_qv = calc_qk(k_xv);
423
424#ifdef VtDEBUG
425 cout << "Vector k_qv: " << k_qv << endl;
426#endif
427
428 return true;
429 } // filter
MATRIX::VtMatrix::nrow
unsigned int nrow() const
no of rows $n$
Definition: VtMatrix.hh:61
MATRIX::VtMatrix::T
const VtMatrix T(void) const
return transpose
Definition: VtMatrix.C:349
MATRIX::VtSymMatrix::invert
bool invert(const bool use_momentum)
calc inverse using momentum or not
Definition: VtSymMatrix.C:265
MATRIX::VtSymMatrix::VtDsinv
bool VtDsinv(int dim=0)
transform to inverse
Definition: VtSymMatrix.C:292
VERTEX::Kalman::calc_AGpc
const MATRIX::VtVector calc_AGpc(void) const
Definition: VtKalman.C:111
VERTEX::Kalman::calc_pc
void calc_pc(double z)
Definition: VtKalman.C:213

◆ filter_chi2()

double VERTEX::Kalman::filter_chi2 ( double  z,
double  prChi2,
const MATRIX::VtSymMatrix &  prCINV,
const MATRIX::VtVector &  prkal_xv 
)
437 {
438
439 k_chi2 = prChi2 + calc_dchi2(z, prCINV, k_xv, prkal_xv, k_qv);
440
441#ifdef VtDEBUG
442 cout << "k_chi2: " << k_chi2 << endl;
443#endif
444 return k_chi2;
445 } // filter_chi2
VERTEX::Kalman::calc_dchi2
double calc_dchi2(double z, const MATRIX::VtSymMatrix &prCINV, const MATRIX::VtVector &xk, const MATRIX::VtVector &prxk, const MATRIX::VtVector &qk) const
Definition: VtKalman.C:236

◆ init()

void VERTEX::Kalman::init ( )

initialize the kalman structure

80 {
81 k_tx = rel->track.tx();
82 k_ty = rel->track.ty();
83 k_qvs = VtVector(rel->track.tx(),rel->track.ty(),rel->track.p());
84 k_qv = k_qvs;
85 k_xv = rel->vertex.xvs();
86 k_C = rel->vertex.CS();
87 k_CINV = k_C.dsinv();
88 k_chi2 = k_chi2s = 0.;
89
90 k_GB.clear();
91
92 return;
93 }
MATRIX::VtMatrix::clear
void clear(void)
set matrix elements to 0
Definition: VtMatrix.C:394
MATRIX::VtSymMatrix::dsinv
const VtSymMatrix dsinv(int dim=0) const
return inverse
Definition: VtSymMatrix.C:275
VERTEX::Relation::vertex
Vertex & vertex
Definition: VtRelation.hh:73
VERTEX::Track::ty
float ty() const
slope (in Vt: p(2,..))
Definition: VtTrack.C:154
VERTEX::Track::p
float p() const
momentum (in Vt: p(5,..))
Definition: VtTrack.C:155
VERTEX::Track::tx
float tx() const
slope (in Vt: p(1,..))
Definition: VtTrack.C:153
VERTEX::Vertex::xvs
const MATRIX::VtVector & xvs() const
vector of vertex-positions
VERTEX::Vertex::CS
const MATRIX::VtSymMatrix & CS() const
$3\times3$ Vertex covariance matrix

◆ inverse_filter()

bool VERTEX::Kalman::inverse_filter ( double  z,
const MATRIX::VtSymMatrix &  CINVn,
const MATRIX::VtSymMatrix &  prCINV,
const MATRIX::VtVector &  kal_xvn 
)
512 {
513
514 // compute C_k^n* (eq. 20)
515 // it is: C_k = (C_k-1^-1 + A_k^TG_k^BA_k)^-1
516 // => A_kG_k^BA_k = C_k^-1 - C_k-1^-1
517 // => C_k^n* = (C_n^-1 - (C_k^-1 - C_k-1^-1))^-1
518 const VtSymMatrix CnkINV = CINVn - (k_CINV - prCINV);
519
520#ifdef VtDEBUG
521 cout << "CnkINV: " << CnkINV << endl;
522#endif
523
524 VtSymMatrix Cnk = CnkINV;
525 if(Cnk.VtDsinv() == false)
526 return false;
527
528 // compute vector p_k - c_k^(0) (in Vt: pc)
529 calc_pc(z);
530
531 // compute Vector (A^T*G^B) * pc (in Vt: AGpc)
532 const VtVector AGpc = calc_AGpc();
533
534 // compute x_k^n* = C_k^n* * [C_n^-1*x_n - A_k^TG_k^B (p_k - c_k^(0))] (eq. 21)
535 // Cx = expression in brackets [...]
536 const VtVector Cx = CINVn * k_xvn - AGpc;
537
538#ifdef VtDEBUG
539 cout << "Vector Cx: " << Cx << endl;
540#endif
541
542 // in Vt: xnk, (eq. 21)
543 k_xnk = Cnk * Cx;
544
545
546 // ------------> inverse filter of vertex-chi2
547 k_chi2s = calc_dchi2(z, prCINV, k_xvn, k_xnk, k_qvs);
548
549
550 // compute x_n - x_k^n* (in Vt: dx)
551 VtVector dx = k_xvn - k_xnk;
552
553 // compute p_k - c_k^(0) - A_k*x_n - B_k * q_k^n (in Vt: dp)
554 // q_k^n are the smoothed track parameters
555 VtVector dp = calc_dp(z, k_xvn, k_qvs);
556
557 unsigned int dim = 4;
558 if(k_use_momentum) {
559 dim = 5;
560 }
561
562 // compute dp^T*G_k*dp
563 double chi2sr = G.product(dp,dim);
564
565 // compute chi2_k,S = dp^T*G_k*dp + dx^T*(C_k^n*)^-1*dx (eq. 23)
566 k_chi2s = CnkINV.product(dx) + chi2sr;
567
568 return true;
569 }

◆ nalpc()

const MATRIX::VtVector & VERTEX::Kalman::nalpc ( ) const

$\vec{n} = (n_x,n_y,n_z)$

◆ p()

double VERTEX::Kalman::p ( ) const

refitted $p$

◆ print()

std::ostream & VERTEX::Kalman::print ( std::ostream &  os) const

called by operator<<()

72 {
73 return os << k_C;
74 }

◆ pvec()

MATRIX::VtVector VERTEX::Kalman::pvec ( ) const

$\vec{v} = (p_x,p_y,p_z)$ refitted mom. vector

◆ px()

double VERTEX::Kalman::px ( ) const

refitted $p_x$

◆ py()

double VERTEX::Kalman::py ( ) const

refitted $p_y$

◆ pz()

double VERTEX::Kalman::pz ( ) const

refitted $p_z$

◆ qvs()

const MATRIX::VtVector & VERTEX::Kalman::qvs ( ) const

$\vec{v} = (t_x,t_y,p)$

◆ qvs_nc()

MATRIX::VtVector & VERTEX::Kalman::qvs_nc ( )

return non-const reference to qvs

◆ set_chi2()

double VERTEX::Kalman::set_chi2 ( const double  chi2)

◆ smooth()

void VERTEX::Kalman::smooth ( double  z,
const MATRIX::VtVector &  xvs,
const MATRIX::VtSymMatrix &  Cn 
)
452 {
453
454 // compute q_k^n = W_k*B_k^T*G_k * (p_k - c_k^(0) - A_k*x_n) (in Vt: qvs)
455 calc_qvs(xvs);
456
457#ifdef VtDEBUG
458 cout << "Vector k_qvs: " << k_qvs << endl;
459 cout << "smooth: Matrix WBG:" << k_WBG << endl;
460#endif
461
462 // compute Matrix B_k^T*G*A_k (eq. 11)
463 VtSqMatrix BGA(3);
464
465 BGA(0,0) = G(2,0) - z * G(0,0);
466 BGA(0,1) = G(2,1) - z * G(0,1);
467 BGA(1,0) = G(3,0) - z * G(1,0);
468 BGA(1,1) = G(3,1) - z * G(1,1);
469
470 BGA(0,2) = -k_tx*BGA(0,0) - k_ty*BGA(0,1);
471 BGA(1,2) = -k_tx*BGA(1,0) - k_ty*BGA(1,1);
472
473 if(k_use_momentum) {
474 BGA(2,0) = G(4,0);
475 BGA(2,1) = G(4,1);
476 BGA(2,2) = -k_tx*BGA(2,0) - k_ty*BGA(2,1);
477 }
478
479#ifdef VtDEBUG
480 cout << "Matrix BGA: " << BGA << endl;
481#endif
482
483 // compute Matrix F_k = W_kB_k^TG_kA_k (eq. 11) (in Vt: F)
484 k_F = k_W * BGA;
485
486#ifdef VtDEBUG
487 cout << "Matrix k_F: " << k_F << endl;
488#endif
489
490 // compute Matrix E_k^n = -F_kC_n (eq. 18) (in Vt: ES)
491 k_ES = -k_F * Cn;
492
493#ifdef VtDEBUG
494 cout << "Marix k_ES: " << k_ES << endl;
495#endif
496
497 // compute Matrix D_k^n = W_k - E_k^n*F_k^T (eq. 18) (in Vt: DS)
498 k_DS = k_W - k_ES * k_F.T();
499
500#ifdef VtDEBUG
501 cout << "Matrix k_DS: " << k_DS << endl;
502#endif
503 return;
504 } // smooth
VERTEX::Kalman::calc_qvs
void calc_qvs(const MATRIX::VtVector &xvs)
Definition: VtKalman.C:203

◆ tvec()

MATRIX::VtVector VERTEX::Kalman::tvec ( ) const

$\vec{v} = (t_x,t_y,1.)$ refitted slope vector

◆ tx()

double VERTEX::Kalman::tx ( ) const

refitted $t_x$

◆ ty()

double VERTEX::Kalman::ty ( ) const

refitted $t_y$

◆ use_momentum() [1/2]

bool VERTEX::Kalman::use_momentum ( ) const

◆ use_momentum() [2/2]

void VERTEX::Kalman::use_momentum ( const bool  use)

set if track momentum should be used in fit

98 {
99 k_use_momentum = use;
100 // select proper inverse cov. matrix
101 if(k_use_momentum)
102 G = rel->track.G();
103 else
104 G = rel->track.GM();
105 return;
106 }
VERTEX::Track::GM
const MATRIX::CMatrix & GM() const
inverse $4\times4$ covariance matrix, without momentum
Definition: VtTrack.C:207

◆ W()

const MATRIX::VtSymMatrix & VERTEX::Kalman::W ( ) const

◆ xn()

double VERTEX::Kalman::xn ( ) const

$n_x = t_x * e_z$

◆ xnk()

const MATRIX::VtVector & VERTEX::Kalman::xnk ( ) const

◆ xv()

const MATRIX::VtVector & VERTEX::Kalman::xv ( ) const

◆ yn()

double VERTEX::Kalman::yn ( ) const

$n_y = t_y * e_z$

◆ zn()

double VERTEX::Kalman::zn ( ) const

$n_z = 1/\sqrt{1 + t_x^2 + t_y^2}$

Member Data Documentation

◆ G

MATRIX::CMatrix VERTEX::Kalman::G
private

◆ k_alp

MATRIX::VtVector VERTEX::Kalman::k_alp
private

◆ k_alpc

MATRIX::VtVector VERTEX::Kalman::k_alpc
private

◆ k_C

MATRIX::VtSymMatrix VERTEX::Kalman::k_C
private

◆ k_chi2

double VERTEX::Kalman::k_chi2
private

◆ k_chi2s

double VERTEX::Kalman::k_chi2s
private

◆ k_CINV

MATRIX::VtSymMatrix VERTEX::Kalman::k_CINV
private

◆ k_DS

MATRIX::VtSymMatrix VERTEX::Kalman::k_DS
private

◆ k_erg

double VERTEX::Kalman::k_erg
private

◆ k_ES

MATRIX::VtSqMatrix VERTEX::Kalman::k_ES
private

◆ k_F

MATRIX::VtSqMatrix VERTEX::Kalman::k_F
private

◆ k_GB

MATRIX::VtMatrix VERTEX::Kalman::k_GB
private

◆ k_Gb

MATRIX::CMatrix VERTEX::Kalman::k_Gb
private

◆ k_nalpc

MATRIX::VtVector VERTEX::Kalman::k_nalpc
private

◆ k_pc

MATRIX::VtVector VERTEX::Kalman::k_pc
private

◆ k_qv

MATRIX::VtVector VERTEX::Kalman::k_qv
private

◆ k_qvs

MATRIX::VtVector VERTEX::Kalman::k_qvs
private

◆ k_tx

double VERTEX::Kalman::k_tx
private

◆ k_ty

double VERTEX::Kalman::k_ty
private

◆ k_use_momentum

bool VERTEX::Kalman::k_use_momentum
private

◆ k_W

MATRIX::VtSymMatrix VERTEX::Kalman::k_W
private

◆ k_WBG

MATRIX::VtMatrix VERTEX::Kalman::k_WBG
private

◆ k_xnk

MATRIX::VtVector VERTEX::Kalman::k_xnk
private

◆ k_xv

MATRIX::VtVector VERTEX::Kalman::k_xv
private

◆ rel

const Relation* VERTEX::Kalman::rel
private
The documentation for this class was generated from the following files: