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  }
t
TTree * t
Definition: check_shower.C:4
Track
Definition: Track.h:10

◆ 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: