Logo ROOT   6.13/01
Reference Guide
MnCovarianceSqueeze.cxx
Go to the documentation of this file.
1 // @(#)root/minuit2:$Id$
2 // Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT *
7  * *
8  **********************************************************************/
9 
12 #include "Minuit2/MinimumError.h"
13 
14 #if defined(DEBUG) || defined(WARNINGMSG)
15 #include "Minuit2/MnPrint.h"
16 #endif
17 
18 
19 namespace ROOT {
20 
21  namespace Minuit2 {
22 
23 
25  // squeeze MnUserCovariance class
26  // MnUserCovariance contasins the error matrix. Need to invert first to get the hessian, then
27  // after having squuezed the hessian, need to invert again to get the new error matrix
28  assert(cov.Nrow() > 0);
29  assert(n < cov.Nrow());
30 
31  MnAlgebraicSymMatrix hess(cov.Nrow());
32  for(unsigned int i = 0; i < cov.Nrow(); i++) {
33  for(unsigned int j = i; j < cov.Nrow(); j++) {
34  hess(i,j) = cov(i,j);
35  }
36  }
37 
38  int ifail = Invert(hess);
39 
40  if(ifail != 0) {
41 #ifdef WARNINGMSG
42  MN_INFO_MSG("MnUserCovariance inversion failed; return diagonal matrix;");
43 #endif
44  MnUserCovariance result(cov.Nrow() - 1);
45  for(unsigned int i = 0, j =0; i < cov.Nrow(); i++) {
46  if(i == n) continue;
47  result(j,j) = cov(i,i);
48  j++;
49  }
50  return result;
51  }
52 
53  MnAlgebraicSymMatrix squeezed = (*this)(hess, n);
54 
55  ifail = Invert(squeezed);
56  if(ifail != 0) {
57 #ifdef WARNINGMSG
58  MN_INFO_MSG("MnUserCovariance back-inversion failed; return diagonal matrix;");
59 #endif
60  MnUserCovariance result(squeezed.Nrow());
61  for(unsigned int i = 0; i < squeezed.Nrow(); i++) {
62  result(i,i) = 1./squeezed(i,i);
63  }
64  return result;
65  }
66 
67  return MnUserCovariance(std::vector<double>(squeezed.Data(), squeezed.Data() + squeezed.size()), squeezed.Nrow());
68 }
69 
70 MinimumError MnCovarianceSqueeze::operator()(const MinimumError& err, unsigned int n) const {
71  // squueze the minimum error class
72  // Remove index-row on the Hessian matrix and the get the new correct error matrix
73  // (inverse of new Hessian)
74  MnAlgebraicSymMatrix hess = err.Hessian();
75  MnAlgebraicSymMatrix squeezed = (*this)(hess, n);
76  int ifail = Invert(squeezed);
77  if(ifail != 0) {
78 #ifdef WARNINGMSG
79  MN_INFO_MSG("MnCovarianceSqueeze: MinimumError inversion fails; return diagonal matrix.");
80 #endif
81  MnAlgebraicSymMatrix tmp(squeezed.Nrow());
82  for(unsigned int i = 0; i < squeezed.Nrow(); i++) {
83  tmp(i,i) = 1./squeezed(i,i);
84  }
86  }
87 
88  return MinimumError(squeezed, err.Dcovar());
89 }
90 
92  // squueze a symmetrix matrix (remove entire row and column n)
93  assert(hess.Nrow() > 0);
94  assert(n < hess.Nrow());
95 
96  MnAlgebraicSymMatrix hs(hess.Nrow() - 1);
97  for(unsigned int i = 0, j = 0; i < hess.Nrow(); i++) {
98  if(i == n) continue;
99  for(unsigned int k = i, l = j; k < hess.Nrow(); k++) {
100  if(k == n) continue;
101  hs(j,l) = hess(i,k);
102  l++;
103  }
104  j++;
105  }
106 
107  return hs;
108 }
109 
110  } // namespace Minuit2
111 
112 } // namespace ROOT
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:22
Namespace for new ROOT classes and functions.
Definition: TFoamSampler.h:19
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
const double * Data() const
Definition: LASymMatrix.h:233
MnAlgebraicSymMatrix Hessian() const
Definition: MinimumError.h:62
unsigned int Nrow() const
Definition: LASymMatrix.h:239
MnUserCovariance operator()(const MnUserCovariance &, unsigned int) const
MinimumError keeps the inv.
Definition: MinimumError.h:26
unsigned int size() const
Definition: LASymMatrix.h:237
Class containing the covariance matrix data represented as a vector of size n*(n+1)/2 Used to hide in...