MultiverseCorrelation.cxx
Go to the documentation of this file.
2 #include "CAFAna/Core/Spectrum.h"
3 
4 #include <iostream>
5 
6 using namespace std;
7 
8 namespace ana
9 {
10 
11  ////////
12  MultiverseCorrelation::MultiverseCorrelation(std::vector<TH1D*> hist){
13  vhistIn.push_back(hist);
14  }
15  ////////
16  MultiverseCorrelation::MultiverseCorrelation(std::vector<ana::Spectrum*> spectra){
17  vhistIn.push_back(SpectraToHist(spectra));
18  }
19  ////////
20  MultiverseCorrelation::MultiverseCorrelation(std::vector<const ana::Spectrum*> spectra){
21  std::vector<ana::Spectrum*> nonconst_spectra;
22  for(unsigned int i = 0; i < spectra.size(); i++)
23  nonconst_spectra.push_back(const_cast<Spectrum*>(spectra[i]));
24  vhistIn.push_back(SpectraToHist(nonconst_spectra));
25  }
26  ////////
27  MultiverseCorrelation::MultiverseCorrelation(std::vector< std::vector<TH1D*> > vhist){
28  for(unsigned int i=0;i<vhist.size();i++)vhistIn.push_back(vhist[i]);
29  }
30  ////////
31  MultiverseCorrelation::MultiverseCorrelation(std::vector< std::vector<ana::Spectrum*> > vspectra){
32  for(unsigned int i=0;i<vspectra.size();i++)vhistIn.push_back(SpectraToHist(vspectra[i]));
33  }
34  ////////
35  MultiverseCorrelation::~MultiverseCorrelation(){}
36 
37  ////////
38  std::vector<TH1D*> MultiverseCorrelation::SpectraToHist(std::vector<ana::Spectrum*> spectra){
39  std::vector<TH1D*> tmp_val;
40  for(unsigned int i = 0; i < spectra.size(); i++)
41  tmp_val.push_back(spectra[i]->ToTH1(spectra[i]->POT()));
42  return tmp_val;
43  }
44  ////////
45  void MultiverseCorrelation::Unifyinghists(){
46  //for now, assuming that all histograms have the same binning:
47  int nbin = vhistIn[0][0]->GetXaxis()->GetNbins();
48  double fval = vhistIn[0][0]->GetXaxis()->GetBinLowEdge(1);
49  double lval = vhistIn[0][0]->GetXaxis()->GetBinUpEdge(nbin);
50  unsigned int nuniv = vhistIn[0].size();
51  unsigned int nvars = vhistIn.size();
52  int nunibin = nvars * nbin;
53  double unifval = fval;
54  double unilval = vhistIn.size() * lval;
55  for(unsigned i=0;i<nuniv;i++){
56  UnifiedhistIn.push_back( new TH1D(Form("uni%d",i),"",nunibin,unifval,unilval));
57  for(unsigned int j=0;j<nvars;j++){
58  for(int k=1;k<=nbin;k++){
59  int current_bin = k + nbin * j;
60  UnifiedhistIn[i]->SetBinContent(current_bin,vhistIn[j][i]->GetBinContent(k));
61  }
62  }
63  }
64 
65  fcov_hist = new TH2D("Cov" ,"",nunibin,unifval,unilval,nunibin,unifval,unilval);
66  fcorr_hist = new TH2D("Corr","",nunibin,unifval,unilval,nunibin,unifval,unilval);
67 
68  }
69  ////////
70  void MultiverseCorrelation::calculate_correlation_matrix(){
71  if(UnifiedhistIn.size()==0)MultiverseCorrelation::calculate_covariance_matrix();
72 
73  for(int i=1;i<=fcorr_hist->GetXaxis()->GetNbins();i++){
74  for(int j=1;j<=fcorr_hist->GetYaxis()->GetNbins();j++){
75  double cov_i_j = fcov_hist->GetBinContent(i,j);
76  double cov_i_i = fcov_hist->GetBinContent(i,i);
77  double cov_j_j = fcov_hist->GetBinContent(j,j);
78  double den = sqrt(cov_i_i*cov_j_j);
79  double valcor = 0;
80  if(den>0)valcor = cov_i_j/den;
81  fcorr_hist->SetBinContent(i,j,valcor);
82  }
83  }
84  }
85 
86  ////////
87  void MultiverseCorrelation::calculate_covariance_matrix(){
88  // Some safeguard.
89  if(!vhistIn[0].size()){
90  cerr << "Number of universes is zero. No correlation is calculated." << endl;
91  return;
92  }
93 
94  MultiverseCorrelation::Unifyinghists();
95 
96  // first calculate the mean
97  TH1D* htemp = (TH1D*)UnifiedhistIn[0]->Clone("htmp");
98  TH1D* hmean = (TH1D*)UnifiedhistIn[0]->Clone("htmp");
99  htemp->Reset();
100  hmean->Reset();
101 
102  for(unsigned int iu=0; iu<UnifiedhistIn.size(); iu++){
103  hmean->Add(UnifiedhistIn[iu]);
104  }
105  hmean->Scale(1.0/double(UnifiedhistIn.size()));
106 
107  //now calculate the covariance matrix
108  TMatrixD fcov_matrix(htemp->GetNbinsX(), htemp->GetNbinsX());
109  for(unsigned int iu=0; iu<UnifiedhistIn.size(); iu++){
110  for(int i=0; i<htemp->GetNbinsX(); i++){
111  double xi=UnifiedhistIn[iu]->GetBinContent(i+1);
112  double ximean=hmean->GetBinContent(i+1);
113  for(int j=i; j<htemp->GetNbinsX(); j++){
114  double xj=UnifiedhistIn[iu]->GetBinContent(j+1);
115  double xjmean=hmean->GetBinContent(j+1);
116  fcov_matrix[i][j]+=(xi-ximean)*(xj-xjmean) ;
117  }
118  }
119  }
120 
121  const double N=double(UnifiedhistIn.size());
122  for(int i=0; i<htemp->GetNbinsX(); i++){
123  for(int j=i; j<htemp->GetNbinsX(); j++){
124  fcov_matrix[i][j]/=(N-1.0);
125  fcov_matrix[j][i]=fcov_matrix[i][j];
126  }
127  }
128 
129  delete hmean;
130  delete htemp;
131 
132  for(int i=1;i<=fcov_hist->GetXaxis()->GetNbins();i++){
133  for(int j=1;j<=fcov_hist->GetYaxis()->GetNbins();j++){
134  fcov_hist->SetBinContent(i,j,fcov_matrix[i-1][j-1]);
135  }
136  }
137 
138  }
139 
140  //////
141  TH2D* MultiverseCorrelation::GetCorrelationMatrix(){
142  calculate_correlation_matrix();
143  return fcorr_hist;
144  }
145 
146  TH2D* MultiverseCorrelation::GetCovarianceMatrix(){
147  calculate_covariance_matrix();
148  return fcov_hist;
149  }
150 
151 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Float_t den
Definition: plot.C:36
T sqrt(T number)
Definition: d0nt_math.hpp:156
OStream cerr
Definition: OStream.cxx:7
const double j
Definition: BetheBloch.cxx:29
unsigned int nuniv
Definition: PlotUnfolding.C:14
TH1F * hmean
Definition: plots_total.C:40
const int nvars