ROOT   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
13
14#if defined(DEBUG) || defined(WARNINGMSG)
15#include "Minuit2/MnPrint.h"
16#endif
17
18
19namespace 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
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
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
const double * Data() const
Definition: LASymMatrix.h:233
unsigned int Nrow() const
Definition: LASymMatrix.h:239
unsigned int size() const
Definition: LASymMatrix.h:237
MinimumError keeps the inv.
Definition: MinimumError.h:26
MnAlgebraicSymMatrix Hessian() const
Definition: MinimumError.h:62
MnUserCovariance operator()(const MnUserCovariance &, unsigned int) const
Class containing the covariance matrix data represented as a vector of size n*(n+1)/2 Used to hide in...
const Int_t n
Definition: legend1.C:16
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:22
VSD Structures.
Definition: StringConv.hxx:21
auto * l
Definition: textangle.C:4