DecorrelateMatrix.h
Go to the documentation of this file.
2 
4 #include "CAFAna/Core/Sample.h"
5 #include "CAFAna/Core/Registry.h"
7 #include "CAFAna/Cuts/QuantileCuts.h"
13 
16 
17 #include "NuXAna/Cuts/NusCuts18.h"
18 #include "NuXAna/Cuts/NusCuts20.h"
22 
23 #include "TFile.h"
24 #include "TKey.h"
25 
26 
27 std::vector<double> DecorrelateND(TMatrixD* inmat, int ndbins, int fdbins){
28 
29  TMatrixD matrix = *inmat;
30 
31  // Create and populate submatrices
32  TMatrixD s11(ndbins, ndbins);
33  TMatrixD s12(fdbins, ndbins);
34  TMatrixD s21(ndbins, fdbins);
35  TMatrixD s22(fdbins, fdbins);
36 
37  for (int i = 0; i < ndbins; ++i) {
38  for (int j = 0; j < ndbins; ++j) {
39  s11(i,j) = matrix(i,j);
40  }
41  }
42 
43  for (int i = 0; i < ndbins; ++i) {
44  for (int j = 0; j < fdbins; ++j) {
45  s12(j,i) = matrix(i,ndbins+j);
46  s21(i,j) = matrix(i,ndbins+j);
47  }
48  }
49 
50  for (int i = 0; i < fdbins; ++i) {
51  for (int j = 0; j < fdbins; ++j) {
52  s22(i,j) = matrix(ndbins+i, ndbins+j);
53  }
54  }
55 
56  // Decorrelate
57  TMatrixD tmp(ndbins, ndbins);
58  s22.Invert();
59  s21 *= s22;
60  tmp.Mult(s21, s12);
61  s11 -= tmp;
62 
63  std::vector<double> ret(ndbins);
64  for (int i = 0; i < ndbins; ++i) {
65  ret[i] = s11(i, i);
66  //std::cout << "nd error bin " << i << ret[i] << std::endl;
67  }
68 
69  return ret;
70 }
71 
72 std::vector<double> DecorrelateFD(TMatrixD* inmat, int ndbins, int fdbins){
73 
74  TMatrixD matrix = *inmat;
75 
76  TMatrixD s11(fdbins, fdbins);
77  TMatrixD s12(ndbins, fdbins);
78  TMatrixD s21(fdbins, ndbins);
79  TMatrixD s22(ndbins, ndbins);
80 
81  for (int i = 0; i < ndbins; ++i) {
82  for (int j = 0; j < ndbins; ++j) {
83  s22(i,j) = matrix(i,j);
84  }
85  }
86 
87  for (int i = 0; i < ndbins; ++i) {
88  for (int j = 0; j < fdbins; ++j) {
89  s12(i,j) = matrix(i,ndbins+j);
90  s21(j,i) = matrix(i,ndbins+j);
91  }
92  }
93 
94  for (int i = 0; i < fdbins; ++i) {
95  for (int j = 0; j < fdbins; ++j) {
96  s11(i,j) = matrix(ndbins+i, ndbins+j);
97  }
98  }
99 
100  TMatrixD tmp(fdbins, fdbins);
101  s22.Invert();
102  s21 *= s22;
103  tmp.Mult(s21, s12);
104  s11 -= tmp;
105 
106  std::vector<double> ret(fdbins);
107  for (int i = 0; i < fdbins; ++i) {
108  ret[i] = s11(i, i);
109  //std::cout << "fd error bin " << i << ret[i] << std::endl;
110  }
111 
112  return ret;
113 
114 }
115 
116 std::vector<double> GetDecorrelatedUncertainty(TMatrixD* inmat,
117  int ndbins, int fdbins){
118 
119  std::vector<double> decorrnd = DecorrelateND(inmat, ndbins, fdbins);
120  std::vector<double> decorrfd = DecorrelateFD(inmat, ndbins, fdbins);
121 
122  std::vector<double> decorrerr;
123  for (size_t i = 0; i < decorrnd.size(); ++i){
124  decorrerr.push_back(std::sqrt(decorrnd.at(i)));
125  }
126  for (size_t i = 0; i < decorrfd.size(); ++i){
127  decorrerr.push_back(std::sqrt(decorrfd.at(i)));
128  }
129 
130  return decorrerr;
131 }
132 
133 
T sqrt(T number)
Definition: d0nt_math.hpp:156
Float_t tmp
Definition: plot.C:36
std::vector< double > DecorrelateND(TMatrixD *inmat, int ndbins, int fdbins)
const double j
Definition: BetheBloch.cxx:29
std::vector< double > DecorrelateFD(TMatrixD *inmat, int ndbins, int fdbins)
std::vector< double > GetDecorrelatedUncertainty(TMatrixD *inmat, int ndbins, int fdbins)