ppfx_check_var.C
Go to the documentation of this file.
1 #include "ppfx_plot_utils.h"
2 #include "ppfx_diag_utils.h"
3 
4 #include "TStopwatch.h"
5 #include "TCanvas.h"
6 #include "TH2.h"
7 #include "TFile.h"
8 #include "TKey.h"
9 #include "TLegend.h"
10 #include "TLine.h"
11 #include "TText.h"
12 #include "TLatex.h"
13 #include "TMatrixD.h"
14 #include "TStyle.h"
15 #include "TColor.h"
16 #include "TDecompSVD.h"
17 #include "TGraph.h"
18 #include "TRandom3.h"
19 #include "TF1.h"
20 #include "TMath.h"
21 #include "TString.h"
22 #include "TPaveText.h"
23 
24 
25 #include <fstream>
26 #include <iostream>
27 #include <cmath>
28 #include <algorithm>
29 
30 
31 void draw_percentiles(TH1* hist, double maxy = 1)
32 {
33  TPaveText* pt = new TPaveText(0.5, 0.25, 0.7, 0.55, "NB NDC");
34  for(int i = 2; i <= hist->GetNbinsX(); i++){
35  if(hist->GetBinContent(i) > 0.85*maxy && hist->GetBinContent(i-1) < 0.85*maxy) {
36  draw_h_line(0, i, hist->GetBinContent(i), kRed);
37  draw_v_line(i, 0, hist->GetBinContent(i), kRed);
38  std::string text85 = "85% Npc : " + std::to_string(i);
39  pt->AddText(text85.c_str());
40  continue;
41  }
42  if(hist->GetBinContent(i) > 0.9*maxy && hist->GetBinContent(i-1) < 0.9*maxy) {
43  draw_h_line(0, i, hist->GetBinContent(i), kRed);
44  draw_v_line(i, 0, hist->GetBinContent(i), kRed);
45  std::string text90 = "90% Npc : " + std::to_string(i);
46  pt->AddText(text90.c_str());
47  continue;
48  }
49  if(hist->GetBinContent(i) > 0.95*maxy && hist->GetBinContent(i-1) < 0.95*maxy) {
50  draw_h_line(0, i, hist->GetBinContent(i), kRed);
51  draw_v_line(i, 0, hist->GetBinContent(i), kRed);
52  std::string text95 = "95% Npc : " + std::to_string(i);
53  pt->AddText(text95.c_str());
54  continue;
55  }
56  }
57  pt->Draw();
58 
59  gPad->Update();
60 }
61 
62 TH1D* get_rms(std::vector<TH1D*> fn_univs, TH1D* fn_cv)
63 {
64  TH1D* h_rms = (TH1D*)fn_cv->Clone("");
65  for(int j = 1; j <= (int) fn_cv->GetNbinsX(); j++){
66  double stddev = 0;
67  for(int i = 0; i < (int) fn_univs.size(); i++)
68  stddev += std::pow(fn_univs[i]->GetBinContent(j)-fn_cv->GetBinContent(j),2)/(fn_univs.size()-1);
69  stddev = (fn_cv->GetBinContent(j)+sqrt(stddev))/fn_cv->GetBinContent(j);
70  h_rms->SetBinContent(j, stddev);
71  }
72 
73  return h_rms;
74 
75 }
76 
77 TH1D* make_pc_univ_rms(std::vector<TH1D*> fn_ratios, TH1D* fn_cv, int nunivs = 5000)
78 {
79  gRandom->SetSeed(0);
80  TH1D* hunit = (TH1D*)fn_cv->Clone("");
81  hunit->Divide(fn_cv);
82  std::vector<TH1D*> fn_pc_univs;
83  TH1D* hzero = (TH1D*)hunit->Clone("");
84  hzero->Add(hunit, -1.);
85 
86  // for(int i = 0; i < nunivs; i++) {
87  // TH1D* fn_pc_univ = (TH1D*) fn_cv->Clone("");
88  // std::for_each(fn_ratios.begin(), fn_ratios.end(),
89  // [&](TH1D* fn_ratio){
90  // TH1D* fn_pc = (TH1D*)fn_ratio->Clone("");
91  // fn_pc->Divide(fn_cv);
92  // fn_pc->Add(hunit, -1);
93  // fn_pc->Scale(gRandom->Gaus());
94  // fn_pc->Add(hunit, +1);
95  // fn_pc_univ->Multiply(fn_pc);
96  // });
97  // fn_pc_univs.push_back(fn_pc_univ);
98  // }
99  //
100  // TH1D* h_fn_pc_rms = get_rms(fn_pc_univs, fn_cv);
101  // return h_fn_pc_rms;
102 
103  TH1D* fn_pc_univ = (TH1D*) fn_cv->Clone("");
104  std::for_each(fn_ratios.begin(), fn_ratios.end(),
105  [&](TH1D* fn_ratio){
106  TH1D* fn_pc = (TH1D*)fn_ratio->Clone("");
107  fn_pc->Divide(fn_cv);
108  fn_pc->Add(hunit, -1);
109  TH1D* fn_pc_clone = (TH1D*)fn_pc->Clone("");
110  fn_pc->Multiply(fn_pc_clone);
111  hzero->Add(fn_pc, +1);
112  });
113  TH1D* h_fn_pc_rms = (TH1D*)hzero->Clone("");
114  for(int i = 1; i <= hzero->GetNbinsX(); i++)
115  h_fn_pc_rms->SetBinContent(i, std::sqrt(hzero->GetBinContent(i)));
116 
117  h_fn_pc_rms->Add(hunit, +1);
118  return h_fn_pc_rms;
119 
120 }
121 
122 double get_metric(TH1D* pc, TH1D* univ, TH1D* weight)
123 {
124  double metric = 0.;
125  double univ_fnrms = 0.;
126  double pc_fnrms = 0.;
127  for(int i = 1; i <= pc->GetNbinsX(); i++){
128  // univ_fnrms += weight->GetBinContent(i)*std::pow((univ->GetBinContent(i)-1),2);
129  // pc_fnrms += weight->GetBinContent(i)*std::pow((pc->GetBinContent(i)-1),2);
130  if(pc->GetBinContent(i) > 0.99*(univ->GetBinContent(i)))
131  metric += 1./pc->GetNbinsX();
132  }
133  // metric = (1+pc_fnrms)/(1+univ_fnrms);
134  return metric;
135 }
136 
137 void draw_fn_coverage(std::vector<TH1D*> fn_univs, std::vector<TH1D*> fn_ratios, TH1D* fn_cv, int npcs, TH1D* near_cv)
138 {
139 
140  // TFile* saveHists = new TFile("ppfx_coverage.root", "RECREATE");
141 
142  std::vector<std::string> labels = {
143  "#nu_{#mu} F/N RHC", "#nu_{e} F/N RHC", "#bar{#nu_{#mu}} F/N RHC", "#bar{#nu_{e}} F/N RHC",
144  "#nu_{#mu} F/N FHC", "#nu_{e} F/N FHC", "#bar{#nu_{#mu}} F/N FHC", "#bar{#nu_{e}} F/N FHC",
145  };
146  std::vector<std::vector<TH1D*>> ppfx_univ_shifts;
147  for(int i = 0; i < (int)fn_univs.size(); i++){
148  TH1D* h = (TH1D*)fn_univs[i]->Clone("");
149  h->Divide(fn_cv);
150  std::vector<TH1D*> hppfx_univ_comps = SplitHistograms(h, labels.size());
151  ppfx_univ_shifts.push_back(hppfx_univ_comps);
152  }
153  TH1D* hppfx_rms = get_rms(fn_univs, fn_cv);
154  std::vector<TH1D*> hppfx_rms_comps = SplitHistograms(hppfx_rms, labels.size());
155 
156  TVectorD metric(npcs);
157 
158  for(int npc = 1; npc <= (int)npcs; npc++){
159  std::vector<TH1D*> fn_pcs(fn_ratios.begin(), fn_ratios.begin()+npc);
160  TH1D* fn_pc_rms = make_pc_univ_rms(fn_pcs, fn_cv);
161  double coverage = get_metric(fn_pc_rms, hppfx_rms, near_cv);
162  metric[npc-1] = coverage;
163  std::cout << coverage << std::endl;
164  std::vector<TH1D*> fn_pc_rms_comps = SplitHistograms(fn_pc_rms, labels.size());
165 
166  TCanvas *c1 = new TCanvas("c1", "c1", 200, 10, 1200, 600);
167  c1->Divide(labels.size()/2, 2);
168  for(int pad = 1; pad <= (int)labels.size(); pad++){
169  c1->cd(pad);
170  for(int ppfx_univ = 0; ppfx_univ < (int)ppfx_univ_shifts.size(); ppfx_univ++){
171  ppfx_univ_shifts[ppfx_univ][pad-1]->Draw("hist same");
172  ppfx_univ_shifts[ppfx_univ][pad-1]->SetTitle(labels[pad-1].c_str());
173  ppfx_univ_shifts[ppfx_univ][pad-1]->GetYaxis()->SetRangeUser(0.9, 1.1);
174  ppfx_univ_shifts[ppfx_univ][pad-1]->SetLineColor(kGray);
175  }
176  hppfx_rms_comps[pad-1]->Draw("hist same");
177  fn_pc_rms_comps[pad-1]->Draw("hist same");
178  fn_pc_rms_comps[pad-1]->SetLineColor(kRed);
179  hppfx_rms_comps[pad-1]->SetLineColor(kGreen+2);
180 
181  if(pad == 1){
182  TLegend* leg = new TLegend(0.5, 0.2, 0.9, 0.4);
183  leg->AddEntry(hppfx_rms_comps[pad-1], "Flux Univs RMS", "le");
184  leg->AddEntry(fn_pc_rms_comps[pad-1], "PC Univs RMS", "le");
185  leg->Draw("same");
186  }
187 
188  }
189  c1->Update();
190  std::string filename = "fnplots/fnrms_beam_quadraturesum_compare_"+std::to_string(npc)+"pc.pdf";
191  c1->Print(filename.c_str());
192  }
193  TH1D* hcover = new TH1D(metric);
194  TCanvas* c = new TCanvas("c", "c", 200, 10, 800, 600);
195  c->cd();
196  hcover->Draw("hist");
197  hcover->SetLineColor(kBlack);
198  draw_percentiles(hcover, 1);
199  hcover->GetYaxis()->SetTitle("Fraction of Bins w/ > 99\% Coverage");
200  hcover->GetXaxis()->SetTitle("N_{pc}");
201  draw_h_line(0, npcs, 1, kBlack);
202  c->Update();
203  std::string filename2 = "fnplots/pca_beam_coverage_metric.pdf";
204  c->Print(filename2.c_str());
205  // saveHists->cd();
206  // saveHists->Close();
207 }
208 
210 {
211 
212  TFile* ppfxFile = new TFile("check_cover_pc_fn_hadp_beam.root", "READ");
213 
214  TH1D* hnear = (TH1D*)ppfxFile->Get("cv/near_cv");
215  std::vector<TH1D*> fn_univs;
216  ppfxFile->cd("fn_ratios_ppfx");
217  TIter next_fn_univ(gDirectory->GetListOfKeys());
218  while(TKey* key = (TKey*)next_fn_univ()) {
219  TH1D* h_fn_univ = (TH1D*)key->ReadObj();
220  fn_univs.push_back(h_fn_univ);
221  }
222  std::cout << fn_univs.size() << std::endl;
223 
224  TH1D* fn_cv = (TH1D*)ppfxFile->Get("cv/fn_cv");
225  std::vector<TH1D*> fn_ratios;
226  ppfxFile->cd("fn_ratios");
227  TIter next_fn_ratio(gDirectory->GetListOfKeys());
228  while(TKey* key = (TKey*)next_fn_ratio()) {
229  TH1D* h_fn_ratio = (TH1D*)key->ReadObj();
230  fn_ratios.push_back(h_fn_ratio);
231  }
232 
233  draw_fn_coverage(fn_univs, fn_ratios, fn_cv, 50, hnear);
234 }
TH1D * fn_cv
std::vector< TH1D * > SplitHistograms(TH1D *hist, int nhists)
TH1D * make_pc_univ_rms(std::vector< TH1D * > fn_ratios, TH1D *fn_cv, int nunivs=5000)
const Var weight
void draw_percentiles(TH1 *hist, double maxy=1)
double maxy
void draw_fn_coverage(std::vector< TH1D * > fn_univs, std::vector< TH1D * > fn_ratios, TH1D *fn_cv, int npcs, TH1D *near_cv)
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
string filename
Definition: shutoffs.py:106
TH1D * get_rms(std::vector< TH1D * > fn_univs, TH1D *fn_cv)
void draw_v_line(double x, double y0, double y1, Color_t color)
void draw_h_line(double x0, double x1, double y, Color_t color)
const double j
Definition: BetheBloch.cxx:29
OStream cout
Definition: OStream.cxx:6
double get_metric(TH1D *pc, TH1D *univ, TH1D *weight)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void ppfx_check_var()
c1
Definition: demo5.py:24
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
int nunivs