Functions | Variables
genie_diag_utils.h File Reference
#include "3FlavorAna/Systs/JointAna2018Systs.h"
#include "TH1.h"
#include "TMatrixD.h"
#include "TMatrixDSymEigen.h"
#include "TMatrixDSym.h"
#include "TVectorD.h"
#include "Math/ProbFunc.h"
#include <fstream>
#include <iostream>
#include <cmath>
#include <algorithm>

Go to the source code of this file.

Functions

TH1D * ErrorBand (TH1D *hNom, std::vector< TH1D * > hUnivs, int nsigmas)
 
TH1D * CombineHistograms (std::vector< TH1D * > hists)
 
std::vector< TH1D * > SplitHistograms (TH1D *hist, int nhists)
 
std::unique_ptr< TMatrixDSym > GetCovMx (const std::vector< TH1D * > &hists)
 
void CrossCheckDiag (TMatrixDSym covMx, TMatrixD V, TMatrixD eigenvectors, TVectorD eigenvalues, double tolerance=1e-14)
 

Variables

double suppress = 1000000
 

Function Documentation

TH1D* CombineHistograms ( std::vector< TH1D * >  hists)

Definition at line 50 of file genie_diag_utils.h.

References bin, analysePickle::hist, MECModelEnuComparisons::i, nbins, runNovaSAM::ret, and art::to_string().

Referenced by ConvertNormalBasis(), genie_syst_pca(), genie_syst_shifts(), GetFNRatio(), ppfx_syst_pca_fn(), and SavePCAShifts().

50  {
51  static int uniq = 0;
52  std::string hist_name = "hist_combined"+std::to_string(uniq);
53 
54  int nbins = 0;
55  for(unsigned int i = 0;i < hists.size();++i) nbins += hists[i]->GetNbinsX();
56 
57  TH1D* ret = new TH1D(hist_name.c_str(),"",nbins,0,nbins);
58 
59  int bin = 1;
60  for(auto hist: hists){
61  for(int i = 1;i <= hist->GetNbinsX();++i){
62  ret->SetBinContent(bin,hist->GetBinContent(i));
63  ++bin;
64  }
65  }
66 
67  ++uniq;
68  return ret;
69 }
TString hists[nhists]
Definition: bdt_com.C:3
const int nbins
Definition: cellShifts.C:15
float bin[41]
Definition: plottest35.C:14
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void CrossCheckDiag ( TMatrixDSym  covMx,
TMatrixD  V,
TMatrixD  eigenvectors,
TVectorD  eigenvalues,
double  tolerance = 1e-14 
)

Definition at line 130 of file genie_diag_utils.h.

References std::abs(), C, om::cout, allTimeWatchdog::endl, MECModelEnuComparisons::i, art::left(), art::right(), std::sqrt(), and stan::math::trace().

Referenced by genie_syst_pca(), and ppfx_syst_pca_fn().

132 {
133  // Number of cross-checks after diagonalizing
134  TMatrixD C = eigenvectors*V;
135  TMatrixD I(TMatrixD::kUnit, V);
136  std::cout << "Checking for Orthogonality.." << std::endl;
137  VerifyMatrixIdentity(C, I, 1, tolerance); // orthogonality
138 
139  std::cout << "Verifying Eigenvectors.." << std::endl;
140  for(int i = 0; i < eigenvalues.GetNrows(); i++){
141  TVectorD left = eigenvectors[i];
142  left *= covMx;
143  TVectorD right = eigenvectors[i];
144  right *= eigenvalues[i];
145  std::cout << "Component " << i << std::endl;
146  VerifyVectorIdentity(left, right, 1, tolerance); // eigenvector relations
147  }
148 
149  std::cout << "Checking if Covariance Matrix Trace = Sum of Eigenvalues.." << std::endl;
150  double trace = 0;
151  double eigensum = 0;
152  for(int i = 0; i < (covMx).GetNrows(); i++) trace += (covMx)(i,i);
153  for(int i = 0; i < eigenvalues.GetNrows(); i++) eigensum += eigenvalues[i];
154  std::cout << "covariance trace : " << sqrt(trace) << " eigenvalue sum : " << sqrt(eigensum) << std::endl;
155  std::cout << std::abs(trace-eigensum) << std::endl;
156  //assert(std::abs(trace-eigensum) < tolerance); // coverage of variances
157  std::cout << "Continuing.." << std::endl;
158 
159 }
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:112
T sqrt(T number)
Definition: d0nt_math.hpp:156
float abs(float number)
Definition: d0nt_math.hpp:39
const double C
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
T trace(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
Definition: trace.hpp:19
TH1D* ErrorBand ( TH1D *  hNom,
std::vector< TH1D * >  hUnivs,
int  nsigmas 
)

Definition at line 19 of file genie_diag_utils.h.

References events(), MECModelEnuComparisons::i, allTimeWatchdog::index, calib::j, cet::sqlite::max(), min(), stan::math::normal_cdf(), elec2geo::pos, runNovaSAM::ret, and art::to_string().

Referenced by FillHists().

19  {
20  TH1D* ret = (TH1D*)hNom->Clone(((std::string)hNom->GetName()+"_"+std::to_string(nsigmas)).c_str());
21 
22  for(int i = 1;i <= hNom->GetNbinsX();++i){
23  double nom = hNom->GetBinContent(i);
24 
25  std::vector<double> events;
26  for(unsigned int j = 0;j < hUnivs.size();++j)
27  events.push_back(hUnivs[j]->GetBinContent(i));
28  sort(events.begin(),events.end());
29 
30  double pos=0;
31  for(unsigned int j = 0;j < events.size()-1;++j)
32  if(nom >= events[j] && nom <= events[j+1]){
33  pos = j+0.5;
34  break;
35  }
36 
37  double count_fraction = ROOT::Math::normal_cdf(nsigmas)-ROOT::Math::normal_cdf(0);
38 
39  int index = 0;
40  int boundIdx = pos + count_fraction*(double)events.size();
41  if(nsigmas >= 0) index = min(boundIdx, (int)events.size()-1);
42  else index = max(boundIdx, 0);
43 
44  ret->SetBinContent(i,events[index]);
45  }
46 
47  return ret;
48 }
const double j
Definition: BetheBloch.cxx:29
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
return_type< T_y, T_loc, T_scale >::type normal_cdf(const T_y &y, const T_loc &mu, const T_scale &sigma)
Definition: normal_cdf.hpp:39
void events(int which)
Definition: Cana.C:52
std::unique_ptr<TMatrixDSym> GetCovMx ( const std::vector< TH1D * > &  hists)

Definition at line 92 of file genie_diag_utils.h.

References ana::assert(), bin, om::cout, allTimeWatchdog::endl, analysePickle::hist, plotROC::nBins, x1, and submit_syst::x2.

Referenced by genie_syst_pca(), and ppfx_syst_pca_fn().

93 {
94  // a bit different from the one in Utilities because I use TMatrixDSym
95  if (hists.size() < 2)
96  return std::unique_ptr<TMatrixDSym>(nullptr);
97 
98  int nBins = hists[0]->GetNbinsX(); // 1 and nBins are inclusive
99 
100  std::vector<double> binMeans(nBins);
101  std::for_each(hists.begin(), hists.end(),
102  [&](TH1D* hist){
103  for(int bin = 1; bin <= nBins; bin++)
104  binMeans[bin-1] += hist->GetBinContent(bin)/hists.size();
105  });
106 
107  auto covMx = std::make_unique<TMatrixDSym>(nBins);
108  std::for_each(hists.begin(), hists.end(),
109  [&](TH1D* hist){
110  for(int bin1 = 1; bin1 <= nBins; bin1++){
111  double x1 = hist->GetBinContent(bin1);
112  for(int bin2 = bin1; bin2 <= nBins; bin2++){
113  double x2 = hist->GetBinContent(bin2);
114  (*covMx)(bin1-1, bin2-1) += (x1-binMeans[bin1-1])*(x2-binMeans[bin2-1])/(hists.size()-1);
115  if(bin1 != bin2)
116  (*covMx)(bin2-1,bin1-1) = (*covMx)(bin1-1,bin2-1);
117  }
118  }
119  });
120 
121 
122 
123  // Check if its symmetric
124  if(covMx->IsSymmetric()) std::cout << "Covariance Matrix is symmetric!" << std::endl;
125  assert(covMx->IsSymmetric());
126 
127  return covMx;
128 }
int nBins
Definition: plotROC.py:16
Float_t x1[n_points_granero]
Definition: compare.C:5
TString hists[nhists]
Definition: bdt_com.C:3
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
assert(nhit_max >=nhit_nbins)
std::vector<TH1D*> SplitHistograms ( TH1D *  hist,
int  nhists 
)

Definition at line 71 of file genie_diag_utils.h.

References bin, hists, MECModelEnuComparisons::i, runNovaSAM::ret, and art::to_string().

Referenced by ConvertNormalBasis(), draw_fn_coverage(), genie_check_var(), genie_syst_shifts(), GetFNRatio(), SavePCAForCAFAna(), and SavePCAShifts().

72 {
73  // get back the flavor spectra
74  static int uniq2 = 0;
75  int ret_bins = hist->GetNbinsX()/nhists;
76  std::vector<TH1D*> hists;
77 
78  for(int i = 0; i < nhists; i++){
79  std::string hist_name = "hist_split"+std::to_string(uniq2);
80  TH1D* ret = new TH1D(hist_name.c_str(), "", ret_bins, 0, ret_bins);
81 
82  for(int bin = 1; bin <= ret_bins; bin++)
83  ret->SetBinContent(bin, hist->GetBinContent(ret_bins*i + bin));
84 
85  hists.push_back(ret);
86  uniq2++;
87  }
88 
89  return hists;
90 }
TString hists[nhists]
Definition: bdt_com.C:3
float bin[41]
Definition: plottest35.C:14
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32

Variable Documentation

double suppress = 1000000

Definition at line 17 of file genie_diag_utils.h.

Referenced by genie_syst_pca(), and genie_syst_shifts().