GENFIT  Rev:NoNumberAvailable
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
Tools.cc
Go to the documentation of this file.
1 /* Copyright 2008-2010, Technische Universitaet Muenchen,
2  Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch
3 
4  This file is part of GENFIT.
5 
6  GENFIT is free software: you can redistribute it and/or modify
7  it under the terms of the GNU Lesser General Public License as published
8  by the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  GENFIT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public License
17  along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #include "Tools.h"
21 
22 #include <cmath>
23 #include <memory>
24 #include <typeinfo>
25 #include <cassert>
26 
27 #include <TDecompChol.h>
28 #include <TMatrixDSymEigen.h>
29 #include <TMatrixTSymCramerInv.h>
30 #include <TMath.h>
31 
32 #include "AbsHMatrix.h"
33 #include "Exception.h"
34 
35 // Use Cramer inversion for small matrices?
36 static const bool useCramer = false;
37 
38 namespace genfit {
39 
40 void tools::invertMatrix(const TMatrixDSym& mat, TMatrixDSym& inv, double* determinant){
41  inv.ResizeTo(mat);
42 
43  // check if numerical limits are reached (i.e at least one entry < 1E-100 and/or at least one entry > 1E100)
44  if (!(mat<1.E100) || !(mat>-1.E100)){
45  Exception e("Tools::invertMatrix() - cannot invert matrix, entries too big (>1e100)",
46  __LINE__,__FILE__);
47  e.setFatal();
48  throw e;
49  }
50  // do the trivial inversions for 1x1 and 2x2 matrices manually
51  if (mat.GetNrows() == 1){
52  if (determinant != NULL) *determinant = mat(0,0);
53  inv(0,0) = 1./mat(0,0);
54  return;
55  }
56 
57  if (mat.GetNrows() == 2){
58  double det = mat(0,0)*mat(1,1) - mat(1,0)*mat(1,0);
59  if (determinant != NULL) *determinant = det;
60  if(fabs(det) < 1E-50){
61  Exception e("Tools::invertMatrix() - cannot invert matrix , determinant = 0",
62  __LINE__,__FILE__);
63  e.setFatal();
64  throw e;
65  }
66  det = 1./det;
67  inv(0,0) = det * mat(1,1);
68  inv(0,1) = inv(1,0) = -det * mat(1,0);
69  inv(1,1) = det * mat(0,0);
70  return;
71  }
72 
73 
74  if (useCramer && mat.GetNrows() <= 6){
75  Bool_t (*inversion)(TMatrixDSym&, Double_t*) = 0;
76  inv.ResizeTo(mat);
77  inv = mat;
78  switch (mat.GetNrows()) {
79  case 3:
80  inversion = TMatrixTSymCramerInv::Inv3x3; break;
81  case 4:
82  inversion = TMatrixTSymCramerInv::Inv4x4; break;
83  case 5:
84  inversion = TMatrixTSymCramerInv::Inv5x5; break;
85  case 6:
86  inversion = TMatrixTSymCramerInv::Inv6x6; break;
87  }
88 
89  Bool_t success = inversion(inv, determinant);
90  if (!success){
91  Exception e("Tools::invertMatrix() - cannot invert matrix, determinant = 0",
92  __LINE__,__FILE__);
93  e.setFatal();
94  throw e;
95  }
96  return;
97  }
98 
99  // else use TDecompChol
100  bool status = 0;
101  TDecompChol invertAlgo(mat, 1E-50);
102 
103  status = invertAlgo.Invert(inv);
104  if(status == 0){
105  Exception e("Tools::invertMatrix() - cannot invert matrix, status = 0",
106  __LINE__,__FILE__);
107  e.setFatal();
108  throw e;
109  }
110 
111  if (determinant != NULL) {
112  double d1, d2;
113  invertAlgo.Det(d1, d2);
114  *determinant = ldexp(d1, d2);
115  }
116 }
117 
118 void tools::invertMatrix(TMatrixDSym& mat, double* determinant){
119  // check if numerical limits are reached (i.e at least one entry < 1E-100 and/or at least one entry > 1E100)
120  if (!(mat<1.E100) || !(mat>-1.E100)){
121  Exception e("Tools::invertMatrix() - cannot invert matrix, entries too big (>1e100)",
122  __LINE__,__FILE__);
123  e.setFatal();
124  throw e;
125  }
126  // do the trivial inversions for 1x1 and 2x2 matrices manually
127  if (mat.GetNrows() == 1){
128  if (determinant != NULL) *determinant = mat(0,0);
129  mat(0,0) = 1./mat(0,0);
130  return;
131  }
132 
133  if (mat.GetNrows() == 2){
134  double *arr = mat.GetMatrixArray();
135  double det = arr[0]*arr[3] - arr[1]*arr[1];
136  if (determinant != NULL) *determinant = det;
137  if(fabs(det) < 1E-50){
138  Exception e("Tools::invertMatrix() - cannot invert matrix, determinant = 0",
139  __LINE__,__FILE__);
140  e.setFatal();
141  throw e;
142  }
143  det = 1./det;
144  double temp[3];
145  temp[0] = det * arr[3];
146  temp[1] = -det * arr[1];
147  temp[2] = det * arr[0];
148  //double *arr = mat.GetMatrixArray();
149  arr[0] = temp[0];
150  arr[1] = arr[2] = temp[1];
151  arr[3] = temp[2];
152  return;
153  }
154 
155  if (useCramer && mat.GetNrows() <= 6){
156  Bool_t (*inversion)(TMatrixDSym&, Double_t*) = 0;
157  switch (mat.GetNrows()) {
158  case 3:
159  inversion = TMatrixTSymCramerInv::Inv3x3; break;
160  case 4:
161  inversion = TMatrixTSymCramerInv::Inv4x4; break;
162  case 5:
163  inversion = TMatrixTSymCramerInv::Inv5x5; break;
164  case 6:
165  inversion = TMatrixTSymCramerInv::Inv6x6; break;
166  }
167 
168  Bool_t success = inversion(mat, determinant);
169  if (!success){
170  Exception e("Tools::invertMatrix() - cannot invert matrix, determinant = 0",
171  __LINE__,__FILE__);
172  e.setFatal();
173  throw e;
174  }
175  return;
176  }
177 
178  // else use TDecompChol
179  bool status = 0;
180  TDecompChol invertAlgo(mat, 1E-50);
181 
182  status = invertAlgo.Invert(mat);
183  if(status == 0){
184  Exception e("Tools::invertMatrix() - cannot invert matrix, status = 0",
185  __LINE__,__FILE__);
186  e.setFatal();
187  throw e;
188  }
189 
190  if (determinant != NULL) {
191  double d1, d2;
192  invertAlgo.Det(d1, d2);
193  *determinant = ldexp(d1, d2);
194  }
195 }
196 
197 
198 // Solves R^T x = b, replaces b with the result x. R is assumed
199 // to be upper-diagonal. This is forward substitution, but with
200 // indices flipped.
201 bool tools::transposedForwardSubstitution(const TMatrixD& R, TVectorD& b)
202 {
203  size_t n = R.GetNrows();
204  double *const bk = b.GetMatrixArray();
205  const double *const Rk = R.GetMatrixArray();
206  for (unsigned int i = 0; i < n; ++i) {
207  double sum = bk[i];
208  for (unsigned int j = 0; j < i; ++j) {
209  sum -= bk[j]*Rk[j*n + i]; // already replaced previous elements in b.
210  }
211  if (Rk[i*n+i] == 0)
212  return false;
213  bk[i] = sum / Rk[i*n + i];
214  }
215  return true;
216 }
217 
218 
219 // Same, but for one column of the matrix b. Used by transposedInvert below
220 // assumes b(i,j) == (i == j)
221 bool tools::transposedForwardSubstitution(const TMatrixD& R, TMatrixD& b, int nCol)
222 {
223  size_t n = R.GetNrows();
224  double *const bk = b.GetMatrixArray() + nCol;
225  const double *const Rk = R.GetMatrixArray();
226  for (unsigned int i = nCol; i < n; ++i) {
227  double sum = (i == (size_t)nCol);
228  for (unsigned int j = 0; j < i; ++j) {
229  sum -= bk[j*n]*Rk[j*n + i]; // already replaced previous elements in b.
230  }
231  if (Rk[i*n+i] == 0)
232  return false;
233  bk[i*n] = sum / Rk[i*n + i];
234  }
235  return true;
236 }
237 
238 
239 // inv will be the inverse of the transposed of the upper-right matrix R
240 bool tools::transposedInvert(const TMatrixD& R, TMatrixD& inv)
241 {
242  bool result = true;
243 
244  inv.ResizeTo(R);
245  double *const invk = inv.GetMatrixArray();
246  int nRows = inv.GetNrows();
247  for (int i = 0; i < nRows; ++i)
248  for (int j = 0; j < nRows; ++j)
249  invk[i*nRows + j] = (i == j);
250 
251  for (int i = 0; i < inv.GetNcols(); ++i)
252  result = result && transposedForwardSubstitution(R, inv, i);
253 
254  return result;
255 }
256 
257 // This replaces A with an upper right matrix connected to A by a
258 // orthogonal transformation. I.e., it computes the R from a QR
259 // decomposition of A replacing A.
260 void tools::QR(TMatrixD& A)
261 {
262  int nCols = A.GetNcols();
263  int nRows = A.GetNrows();
264  assert(nRows >= nCols);
265  // This uses Businger and Golub's algorithm from Handbook for
266  // Automatical Computation, Vol. 2, Chapter 8, but without
267  // pivoting. I.e., we stop at the middle of page 112. We don't
268  // explicitly calculate the orthogonal matrix.
269 
270  double *const ak = A.GetMatrixArray();
271  // No variable-length arrays in C++, alloca does the exact same thing ...
272  double *const u = (double *)alloca(sizeof(double)*nRows);
273 
274  // Main loop over matrix columns.
275  for (int k = 0; k < nCols; ++k) {
276  double akk = ak[k*nCols + k];
277 
278  double sum = akk*akk;
279  // Put together a housholder transformation.
280  for (int i = k + 1; i < nRows; ++i) {
281  sum += ak[i*nCols + k]*ak[i*nCols + k];
282  u[i] = ak[i*nCols + k];
283  }
284  double sigma = sqrt(sum);
285  double beta = 1/(sum + sigma*fabs(akk));
286  // The algorithm uses only the uk[i] for i >= k.
287  u[k] = copysign(sigma + fabs(akk), akk);
288 
289  // Calculate y (again taking into account zero entries). This
290  // encodes how the (sub)matrix changes by the householder transformation.
291  for (int i = k; i < nCols; ++i) {
292  double y = 0;
293  for (int j = k; j < nRows; ++j)
294  y += u[j]*ak[j*nCols + i];
295  y *= beta;
296  // ... and apply the changes.
297  for (int j = k; j < nRows; ++j)
298  ak[j*nCols + i] -= u[j]*y; //y[j];
299  }
300  }
301 
302  // Zero below diagonal
303  for (int i = 1; i < nCols; ++i)
304  for (int j = 0; j < i; ++j)
305  ak[i*nCols + j] = 0.;
306  for (int i = nCols; i < nRows; ++i)
307  for (int j = 0; j < nCols; ++j)
308  ak[i*nCols + j] = 0.;
309 }
310 
311 // This replaces A with an upper right matrix connected to A by a
312 // orthogonal transformation. I.e., it computes the R from a QR
313 // decomposition of A replacing A. Simultaneously it transforms b by
314 // the inverse orthogonal transformation.
315 //
316 // The purpose is this: the least-squared problem
317 // ||Ax - b|| = min
318 // is equivalent to
319 // ||QRx - b|| = ||Rx - Q'b|| = min
320 // where Q' denotes the transposed (i.e. inverse).
321 void tools::QR(TMatrixD& A, TVectorD& b)
322 {
323  int nCols = A.GetNcols();
324  int nRows = A.GetNrows();
325  assert(nRows >= nCols);
326  assert(b.GetNrows() == nRows);
327  // This uses Businger and Golub's algorithm from Handbook for
328  // Automatic Computation, Vol. 2, Chapter 8, but without pivoting.
329  // I.e., we stop at the middle of page 112. We don't explicitly
330  // calculate the orthogonal matrix, but Q'b which is not done
331  // explicitly in Businger et al.
332  // Also in Numer. Math. 7, 269-276 (1965)
333 
334  double *const ak = A.GetMatrixArray();
335  double *const bk = b.GetMatrixArray();
336  // No variable-length arrays in C++, alloca does the exact same thing ...
337  //double * u = (double *)alloca(sizeof(double)*nRows);
338  double u[500];
339 
340  // Main loop over matrix columns.
341  for (int k = 0; k < nCols; ++k) {
342  double akk = ak[k*nCols + k];
343 
344  double sum = akk*akk;
345  // Put together a housholder transformation.
346  for (int i = k + 1; i < nRows; ++i) {
347  sum += ak[i*nCols + k]*ak[i*nCols + k];
348  u[i] = ak[i*nCols + k];
349  }
350  double sigma = sqrt(sum);
351  double beta = 1/(sum + sigma*fabs(akk));
352  // The algorithm uses only the uk[i] for i >= k.
353  u[k] = copysign(sigma + fabs(akk), akk);
354 
355  // Calculate b (again taking into account zero entries). This
356  // encodes how the (sub)vector changes by the householder transformation.
357  double yb = 0;
358  for (int j = k; j < nRows; ++j)
359  yb += u[j]*bk[j];
360  yb *= beta;
361  // ... and apply the changes.
362  for (int j = k; j < nRows; ++j)
363  bk[j] -= u[j]*yb;
364 
365  // Calculate y (again taking into account zero entries). This
366  // encodes how the (sub)matrix changes by the householder transformation.
367  for (int i = k; i < nCols; ++i) {
368  double y = 0;
369  for (int j = k; j < nRows; ++j)
370  y += u[j]*ak[j*nCols + i];
371  y *= beta;
372  // ... and apply the changes.
373  for (int j = k; j < nRows; ++j)
374  ak[j*nCols + i] -= u[j]*y;
375  }
376  }
377 
378  // Zero below diagonal
379  for (int i = 1; i < nCols; ++i)
380  for (int j = 0; j < i; ++j)
381  ak[i*nCols + j] = 0.;
382  for (int i = nCols; i < nRows; ++i)
383  for (int j = 0; j < nCols; ++j)
384  ak[i*nCols + j] = 0.;
385 }
386 
387 // This averages the covariance matrices C1, C2 in a numerically
388 // stable way by using matrix square roots. No optimizations
389 // performed, so use with care.
390 void tools::safeAverage(const TMatrixDSym& C1, const TMatrixDSym& C2,
391  TMatrixDSym& result)
392 {
393  /*
394  The algorithm proceeds as follows:
395  write C1 = S1 S1' (prime for transpose),
396  C2 = S2 S2'
397  Then the inverse of the average can be written as ("." for matrix
398  multiplication)
399  C^-1 = ((S1'^-1, S2'^-1) . (S1'^-1) )
400  ( (S2'^-1) )
401  Inserting an orthogonal matrix T in the middle:
402  C^-1 = ((S1'^-1, S2'^-1) . T . T' . (S1'^-1) )
403  ( (S2'^-1) )
404  doesn't change this because T.T' = 1.
405  Now choose T s.t. T'.(S1'^-1, S2'^-1)' is an upper right matrix. We
406  use Tools::QR for the purpose, as we don't actually need T.
407 
408  Then the inverse needed to obtain the covariance matrix can be
409  obtained by inverting the upper right matrix, which is squared to
410  obtained the new covariance matrix. */
411  TDecompChol dec1(C1);
412  dec1.Decompose();
413  TDecompChol dec2(C2);
414  dec2.Decompose();
415 
416  const TMatrixD& S1 = dec1.GetU();
417  const TMatrixD& S2 = dec2.GetU();
418 
419  TMatrixD S1inv, S2inv;
420  transposedInvert(S1, S1inv);
421  transposedInvert(S2, S2inv);
422 
423  TMatrixD A(2 * S1.GetNrows(), S1.GetNcols());
424  for (int i = 0; i < S1.GetNrows(); ++i) {
425  for (int j = 0; j < S2.GetNcols(); ++j) {
426  A(i, j) = S1inv(i, j);
427  A(i + S1.GetNrows(), j) = S2inv(i, j);
428  }
429  }
430 
431  QR(A);
432  A.ResizeTo(S1.GetNrows(), S1.GetNrows());
433 
434  TMatrixD inv;
435  transposedInvert(A, inv);
436 
437  result.ResizeTo(inv.GetNcols(), inv.GetNcols());
438  result = TMatrixDSym(TMatrixDSym::kAtA, inv);
439 }
440 
441 
442 void
443 tools::noiseMatrixSqrt(const TMatrixDSym& noise,
444  TMatrixD& noiseSqrt)
445 {
446  // This is the slowest part of the whole Sqrt Kalman. Using an LDLt
447  // transform is probably the easiest way of remedying this.
448  TMatrixDSymEigen eig(noise);
449  noiseSqrt.ResizeTo(noise);
450  noiseSqrt = eig.GetEigenVectors();
451  double* pNoiseSqrt = noiseSqrt.GetMatrixArray();
452  const TVectorD& evs(eig.GetEigenValues());
453  const double* pEvs = evs.GetMatrixArray();
454  // GetEigenVectors is such that noise = noiseSqrt * evs * noiseSqrt'
455  // We're evaluating the first product with the eigenvalues replaced
456  // by their square roots, so we're multiplying with a diagonal
457  // matrix from the right.
458  int iCol = 0;
459  for (; iCol < noiseSqrt.GetNrows(); ++iCol) {
460  double ev = pEvs[iCol] > 0 ? sqrt(pEvs[iCol]) : 0;
461  // if (ev == 0)
462  // break;
463  for (int j = 0; j < noiseSqrt.GetNrows(); ++j) {
464  pNoiseSqrt[j*noiseSqrt.GetNcols() + iCol] *= ev;
465  }
466  }
467  if (iCol < noiseSqrt.GetNcols()) {
468  // Hit zero eigenvalue, resize matrix
469  noiseSqrt.ResizeTo(noiseSqrt.GetNrows(), iCol);
470  }
471 
472  // noiseSqrt * noiseSqrt' = noise
473 }
474 
475 // Transports the state.
476 void
477 tools::kalmanPrediction(const TVectorD& x,
478  const TVectorD& delta, const TMatrixD& F,
479  TVectorD& xNew)
480 {
481  xNew = x;
482  xNew *= F;
483  xNew += delta;
484 }
485 
486 // Transports the square root of the covariance matrix using a
487 // square-root formalism
488 //
489 // With covariance square root S, transport matrix F and noise matrix
490 // square root Q.
491 void
493  const TMatrixD& F, const TMatrixD& Q,
494  TMatrixD& Snew)
495 {
496  Snew.ResizeTo(S.GetNrows() + Q.GetNcols(),
497  S.GetNcols());
498 
499  // This overwrites all elements, no precautions necessary
500  Snew.SetSub(0, 0, TMatrixD(S, TMatrixD::kMultTranspose, F));
501  if (Q.GetNcols() != 0)
502  Snew.SetSub(S.GetNrows(), 0, TMatrixD(TMatrixD::kTransposed, Q));
503 
504  tools::QR(Snew);
505 
506  // The result is in the upper right corner of the matrix.
507  Snew.ResizeTo(S.GetNrows(), S.GetNrows());
508 }
509 
510 
511 // Kalman measurement update (no transport)
512 // x, S : state prediction, covariance square root
513 // res, R, H : residual, measurement covariance square root, H matrix of the measurement
514 // gives the update (new state = x + update) and the updated covariance square root.
515 // S and Snew are allowed to refer to the same object.
516 void
517 tools::kalmanUpdateSqrt(const TMatrixD& S,
518  const TVectorD& res, const TMatrixD& R,
519  const AbsHMatrix* H,
520  TVectorD& update, TMatrixD& SNew)
521 {
522  TMatrixD pre(S.GetNrows() + R.GetNrows(),
523  S.GetNcols() + R.GetNcols());
524  pre.SetSub(0, 0, R); /* Zeros in upper right block */
525  pre.SetSub(R.GetNrows(), 0, H->MHt(S)); pre.SetSub(R.GetNrows(), R.GetNcols(), S);
526 
527  tools::QR(pre);
528  const TMatrixD& r = pre;
529 
530  const TMatrixD& a(r.GetSub(0, R.GetNrows()-1,
531  0, R.GetNcols()-1));
532  TMatrixD K(TMatrixD::kTransposed, r.GetSub(0, R.GetNrows()-1, R.GetNcols(), pre.GetNcols()-1));
533  SNew = r.GetSub(R.GetNrows(), pre.GetNrows()-1, R.GetNcols(), pre.GetNcols()-1);
534 
535  update.ResizeTo(res);
536  update = res;
538  update *= K;
539 }
540 
541 
542 // Kalman transport + measurement update
543 // S : covariance square root (pre-prediction)
544 // transport matrix F and noise matrix square root Q
545 // res, R, H : residual, measurement covariance square root, H matrix of the measurement
546 // The new state is xnew = F*xold + update
547 void
549  const TMatrixD& F, const TMatrixD& Q,
550  const TVectorD& res, const TMatrixD& R,
551  const AbsHMatrix* H,
552  TVectorD& update, TMatrixD& SNew)
553 {
554  TMatrixD pre(S.GetNrows() + Q.GetNrows() + R.GetNrows(),
555  S.GetNcols() + R.GetNcols());
556  TMatrixD SFt(S, TMatrixD::kMultTranspose, F);
557  pre.SetSub( 0, 0, R); /* upper right block is zero */
558  pre.SetSub( R.GetNrows(), 0,H->MHt(SFt)); pre.SetSub(R.GetNrows(), R.GetNcols(),SFt);
559  if (Q.GetNcols() > 0) { // needed to suppress warnings when inserting an empty Q
560  TMatrixD Qt(TMatrixD::kTransposed, Q);
561  pre.SetSub(S.GetNrows()+R.GetNrows(),0,H->MHt(Qt)); pre.SetSub(S.GetNrows()+R.GetNrows(),R.GetNcols(), Qt);
562  }
563 
564  tools::QR(pre);
565  const TMatrixD& r = pre;
566 
567  TMatrixD a(r.GetSub(0, R.GetNrows()-1, 0, R.GetNcols()-1));
568  TMatrixD K(TMatrixD::kTransposed, r.GetSub(0, R.GetNrows()-1, R.GetNcols(), pre.GetNcols()-1));
569  SNew = r.GetSub(R.GetNrows(), R.GetNrows() + S.GetNrows() - 1,
570  R.GetNcols(), pre.GetNcols() - 1);
571 
572  update.ResizeTo(res);
573  update = res;
575  update *= K;
576 }
577 
578 
579 } /* End of namespace genfit */
void kalmanPrediction(const TVectorD &x, const TVectorD &delta, const TMatrixD &F, TVectorD &xNew)
Transport the state.
Definition: Tools.cc:477
void kalmanPredictionUpdateSqrt(const TMatrixD &S, const TMatrixD &F, const TMatrixD &Q, const TVectorD &res, const TMatrixD &R, const AbsHMatrix *H, TVectorD &update, TMatrixD &SNew)
Definition: Tools.cc:548
void noiseMatrixSqrt(const TMatrixDSym &noise, TMatrixD &noiseSqrt)
Calculate a sqrt for the positive semidefinite noise matrix. Rows corresponding to zero eigenvalues a...
Definition: Tools.cc:443
static const bool useCramer
Definition: Tools.cc:36
void QR(TMatrixD &A)
Replaces A with an upper right matrix connected to A by an orthongonal transformation. I.e., it computes R from a QR decomposition of A = QR, replacing A.
Definition: Tools.cc:260
HMatrix for projecting from AbsTrackRep parameters to measured parameters in a DetPlane.
Definition: AbsHMatrix.h:37
bool transposedForwardSubstitution(const TMatrixD &R, TVectorD &b)
Solves R^t x = b, replacing b with the solution for x. R is assumed to be upper diagonal.
Definition: Tools.cc:201
void kalmanUpdateSqrt(const TMatrixD &S, const TVectorD &res, const TMatrixD &R, const AbsHMatrix *H, TVectorD &update, TMatrixD &SNew)
Calculate the Kalman measurement update with no transport. x, S : state prediction, covariance square root res, R, H : residual, measurement covariance square root, H matrix of the measurement.
Definition: Tools.cc:517
virtual TMatrixD MHt(const TMatrixDSym &M) const
M*H^t.
Definition: AbsHMatrix.h:52
void safeAverage(const TMatrixDSym &C1, const TMatrixDSym &C2, TMatrixDSym &result)
This averages the covariance matrices C1, C2 in a numerically stable way by using matrix square roots...
Definition: Tools.cc:390
void setFatal(bool b=true)
Set fatal flag.
Definition: Exception.h:61
void kalmanPredictionCovSqrt(const TMatrixD &S, const TMatrixD &F, const TMatrixD &Q, TMatrixD &Snew)
Calculates the square root of the covariance matrix after the Kalman prediction (i.e. extrapolation) with transport matrix F and the noise square root Q. Gives the new covariance square root.
Definition: Tools.cc:492
bool transposedInvert(const TMatrixD &R, TMatrixD &inv)
Inverts the transpose of the upper right matrix R into inv.
Definition: Tools.cc:240
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: Exception.h:48
void invertMatrix(const TMatrixDSym &mat, TMatrixDSym &inv, double *determinant=NULL)
Invert a matrix, throwing an Exception when inversion fails. Optional calculation of determinant...
Definition: Tools.cc:40
Defines for I/O streams used for error and debug printing.