7 #include "TMatrixDSymEigen.h" 8 #include "TMatrixDSym.h" 10 #include "Math/ProbFunc.h" 19 TH1D*
ErrorBand(TH1D* hNom,std::vector<TH1D*> hUnivs,
int nsigmas){
22 for(
int i = 1;
i <= hNom->GetNbinsX();++
i){
23 double nom = hNom->GetBinContent(
i);
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());
31 for(
unsigned int j = 0;
j < events.size()-1;++
j)
32 if(nom >= events[
j] && nom <= events[
j+1]){
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);
44 ret->SetBinContent(
i,events[index]);
55 for(
unsigned int i = 0;
i < hists.size();++
i) nbins += hists[
i]->GetNbinsX();
60 for(
auto hist: hists){
61 for(
int i = 1;
i <=
hist->GetNbinsX();++
i){
62 ret->SetBinContent(bin,
hist->GetBinContent(
i));
75 int ret_bins = hist->GetNbinsX()/nhists;
76 std::vector<TH1D*>
hists;
78 for(
int i = 0;
i < nhists;
i++){
80 TH1D*
ret =
new TH1D(hist_name.c_str(),
"", ret_bins, 0, ret_bins);
83 ret->SetBinContent(
bin, hist->GetBinContent(ret_bins*
i +
bin));
96 return std::unique_ptr<TMatrixDSym>(
nullptr);
98 int nBins = hists[0]->GetNbinsX();
100 std::vector<double> binMeans(nBins);
101 std::for_each(hists.begin(), hists.end(),
104 binMeans[
bin-1] +=
hist->GetBinContent(
bin)/hists.size();
107 auto covMx = std::make_unique<TMatrixDSym>(
nBins);
108 std::for_each(hists.begin(), hists.end(),
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);
116 (*covMx)(bin2-1,bin1-1) = (*covMx)(bin1-1,bin2-1);
125 assert(covMx->IsSymmetric());
137 VerifyMatrixIdentity(C, I, 1, tolerance);
140 for(
int i = 0;
i < eigenvalues.GetNrows();
i++){
144 right *= eigenvalues[
i];
146 VerifyVectorIdentity(left, right, 1, tolerance);
149 std::cout <<
"Checking if Covariance Matrix Trace = Sum of Eigenvalues.." <<
std::endl;
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];
std::vector< TH1D * > SplitHistograms(TH1D *hist, int nhists)
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Float_t x1[n_points_granero]
TH1D * CombineHistograms(std::vector< TH1D * > hists)
static float min(const float a, const float b, const float c)
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
std::unique_ptr< TMatrixDSym > GetCovMx(const std::vector< TH1D * > &hists)
T trace(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m)
assert(nhit_max >=nhit_nbins)
std::string to_string(ModuleType mt)
void CrossCheckDiag(TMatrixDSym covMx, TMatrixD V, TMatrixD eigenvectors, TVectorD eigenvalues, double tolerance=1e-14)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
TH1D * ErrorBand(TH1D *hNom, std::vector< TH1D * > hUnivs, int nsigmas)
return_type< T_y, T_loc, T_scale >::type normal_cdf(const T_y &y, const T_loc &mu, const T_scale &sigma)