ppfx_diag_utils.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include "TH1.h"
4 #include "TMatrixD.h"
5 #include "TMatrixDSymEigen.h"
6 #include "TMatrixDSym.h"
7 #include "TVectorD.h"
8 
9 #include <fstream>
10 #include <iostream>
11 #include <cmath>
12 #include <algorithm>
13 #include <memory>
14 
15 TH1D* CombineHistograms(std::vector<TH1D*> hists)
16 {
17  static int uniq = 0;
18  std::string hist_name = "hist_combined"+std::to_string(uniq);
19  int nbins = 0;
20  std::for_each(hists.begin(), hists.end(), [&](TH1D* hist){ nbins = nbins + hist->GetNbinsX(); });
21  TH1D* ret = new TH1D(hist_name.c_str(), "hist_combined", nbins, 0, nbins);
22 
23  int bin = 1;
24  for(auto hist: hists){
25  for(int binIdx = 1; binIdx <= hist->GetNbinsX(); binIdx++){
26  ret->SetBinContent(bin, hist->GetBinContent(binIdx));
27  bin++;
28  }
29  }
30  uniq++;
31  return ret;
32 
33 }
34 
35 std::vector<TH1D*> SplitHistograms(TH1D* hist, int nhists)
36 {
37  static int uniq2 = 0;
38  int ret_bins = hist->GetNbinsX()/nhists;
39  std::vector<TH1D*> hists;
40 
41  for(int i = 0; i < nhists; i++){
42  std::string hist_name = "hist_split"+std::to_string(uniq2);
43  TH1D* ret = new TH1D(hist_name.c_str(), "hist_split", ret_bins, 0, ret_bins);
44 
45  for(int bin = 1; bin <= ret_bins; bin++)
46  ret->SetBinContent(bin, hist->GetBinContent(ret_bins*i + bin));
47 
48  hists.push_back(ret);
49  uniq2++;
50  }
51 
52  return hists;
53 }
54 
55 TH1D* AddHistograms(std::vector<TH1D*> hists)
56 {
57  TH1D* ret = (TH1D*)hists[0]->Clone("");
58  for(int i = 1; i < (int)hists.size(); i++)
59  ret->Add(hists[i]);
60 
61  return ret;
62 }
63 
64 std::unique_ptr<TMatrixDSym> GetCovMx(const std::vector<TH1D*> &hists, TH1D* hist_cv)
65 {
66  if (hists.size() < 2)
67  return std::unique_ptr<TMatrixDSym>(nullptr);
68 
69  int nBins = hists[0]->GetNbinsX(); // 1 and nBins are inclusive
70 
71  // std::vector<double> binMeans(nBins);
72  // std::for_each(hists.begin(), hists.end(),
73  // [&](TH1D* hist){
74  // for(int bin = 1; bin <= nBins; bin++)
75  // binMeans[bin-1] += hist->GetBinContent(bin)/hists.size();
76  // });
77  //
78  // auto covMx = std::make_unique<TMatrixDSym>(nBins);
79  // std::for_each(hists.begin(), hists.end(),
80  // [&](TH1D* hist){
81  // for(int bin1 = 1; bin1 <= nBins; bin1++){
82  // double x1 = hist->GetBinContent(bin1);
83  // for(int bin2 = bin1; bin2 <= nBins; bin2++){
84  // double x2 = hist->GetBinContent(bin2);
85  // (*covMx)(bin1-1, bin2-1) += (x1-binMeans[bin1-1])*(x2-binMeans[bin2-1])/(hists.size()-1);
86  // if(bin1 != bin2)
87  // (*covMx)(bin2-1,bin1-1) = (*covMx)(bin1-1,bin2-1);
88  // }
89  // }
90  // });
91  auto covMx = std::make_unique<TMatrixDSym>(nBins);
92  std::for_each(hists.begin(), hists.end(),
93  [&](TH1D* hist){
94  for(int bin1 = 1; bin1 <= nBins; bin1++){
95  double x1 = hist->GetBinContent(bin1);
96  for(int bin2 = bin1; bin2 <= nBins; bin2++){
97  double x2 = hist->GetBinContent(bin2);
98  (*covMx)(bin1-1, bin2-1) += (x1-hist_cv->GetBinContent(bin1))*(x2-hist_cv->GetBinContent(bin2))/(hists.size()-1);
99  if(bin1 != bin2)
100  (*covMx)(bin2-1,bin1-1) = (*covMx)(bin1-1,bin2-1);
101  }
102  }
103  });
104 
105 
106 
107  // Check if its symmetric
108  if(covMx->IsSymmetric()) std::cout << "Covariance Matrix is symmetric!" << std::endl;
109  assert(covMx->IsSymmetric());
110 
111  return covMx;
112 }
113 
114 void CrossCheckDiag(TMatrixDSym covMx, TMatrixD V,
115  TMatrixD eigenvectors, TVectorD eigenvalues, double tolerance = 1e-14)
116 {
117  TMatrixD C = eigenvectors*V;
118  TMatrixD I(TMatrixD::kUnit, V);
119  std::cout << "Checking for Orthogonality.." << std::endl;
120  VerifyMatrixIdentity(C, I, 1, tolerance); // orthogonality
121 
122  std::cout << "Verifying Eigenvectors.." << std::endl;
123  for(int i = 0; i < eigenvalues.GetNrows(); i++){
124  TVectorD left = eigenvectors[i];
125  left *= covMx;
126  TVectorD right = eigenvectors[i];
127  right *= eigenvalues[i];
128  std::cout << "Component " << i << std::endl;
129  VerifyVectorIdentity(left, right, 1, tolerance); // eigenvector relations
130  }
131 
132  std::cout << "Checking if Covariance Matrix Trace = Sum of Eigenvalues.." << std::endl;
133  double trace = 0;
134  double eigensum = 0;
135  for(int i = 0; i < (covMx).GetNrows(); i++) trace += (covMx)(i,i);
136  for(int i = 0; i < eigenvalues.GetNrows(); i++) eigensum += eigenvalues[i];
137  std::cout << "covariance trace : " << sqrt(trace) << " eigenvalue sum : " << sqrt(eigensum) << std::endl;
138  assert(std::abs(trace-eigensum) < tolerance); // coverage of variances
139  std::cout << "Continuing.." << std::endl;
140 
141 }
142 
int nBins
Definition: plotROC.py:16
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:112
Float_t x1[n_points_granero]
Definition: compare.C:5
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH1D * CombineHistograms(std::vector< TH1D * > hists)
std::unique_ptr< TMatrixDSym > GetCovMx(const std::vector< TH1D * > &hists, TH1D *hist_cv)
TString hists[nhists]
Definition: bdt_com.C:3
void CrossCheckDiag(TMatrixDSym covMx, TMatrixD V, TMatrixD eigenvectors, TVectorD eigenvalues, double tolerance=1e-14)
float abs(float number)
Definition: d0nt_math.hpp:39
const double C
const int nbins
Definition: cellShifts.C:15
std::vector< TH1D * > SplitHistograms(TH1D *hist, int nhists)
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:104
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T trace(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
Definition: trace.hpp:19
assert(nhit_max >=nhit_nbins)
TH1D * AddHistograms(std::vector< TH1D * > hists)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
Float_t e
Definition: plot.C:35