xsec_uncertainty_per_bin.C
Go to the documentation of this file.
1 #include "Utilities/rootlogon.C"
3 #include "CAFAna/Core/Spectrum.h"
6 #include <iostream>
7 #include <fstream>
8 #include <vector>
9 #include <cmath>
10 #include <stdint.h>
11 #include <string>
12 
13 #include "TFile.h"
14 #include "TTree.h"
15 #include "TCanvas.h"
16 #include "TF1.h"
17 #include "TH1F.h"
18 #include "TH2F.h"
19 #include "TStyle.h"
20 #include "TLegend.h"
21 #include "TGaxis.h"
22 
23 #include <vector>
25 #include "CAFAna/Core/Spectrum.h"
26 #include "CAFAna/Core/SystShifts.h"
27 #include "CAFAna/Analysis/Plots.h"
28 #include "CAFAna/Core/HistAxis.h"
29 #include "CAFAna/Cuts/SpillCuts.h"
30 #include "CAFAna/Cuts/TruthCuts.h"
31 //#include "3FlavorAna/Cuts/NueCutsSecondAna.h"
32 
33 //#include "3FlavorAna/Vars/NueVars.h"
34 //#include "3FlavorAna/Vars/NueVarsExtra.h"
37 
39 #include "CAFAna/Fit/Fit.h"
40 
41 //#include "NDAna/nuecc_inc/NueCCIncVars.h"
42 //#include "NDAna/nuecc_inc/NueCCIncCuts.h"
43 //#include "NDAna/my_nuecc_inc/NueCCIncBins.h"
44 //#include "CAFAna/XSec/FluxMultiverseSyst.h"
46 
47 #include "TH1.h"
48 #include "TH2.h"
49 #include "THStack.h"
50 #include "TCanvas.h"
51 #include "TLegend.h"
52 #include "TROOT.h"
53 #include <iostream>
54 #include <vector>
55 #include <string>
56 
57 #include "CAFAna/Cuts/SpillCuts.h"
58 #include "CAFAna/Core/Utilities.h"
60 #include "CAFAna/Core/Spectrum.h"
61 #include "CAFAna/Cuts/TimingCuts.h"
62 #include "CAFAna/Analysis/Plots.h"
63 #include "CAFAna/Analysis/Style.h"
65 #include "CAFAna/Core/Ratio.h"
66 
67 #include "TH1D.h"
68 #include "TCanvas.h"
69 #include "TGraph.h"
70 #include "TMultiGraph.h"
71 #include "TLine.h"
72 #include "TLegend.h"
73 
74 using namespace ana;
75 
76 const double pot = 8.09e20; // Data POT for ND after 3A
77 const double mc_pot = 3.54206e+21; //MC POT for ND after 3A
78 //const double mc_syst_pot = mc_pot
79 
80 const bool up_down_combined = true; // If true, do up-down for Genie, Calib and Light syst.
81  // Otherwised, treat all up and down shifts individually
82 
83 //-----------------------------
84 TGraphAsymmErrors* MyPlotWithSystErrorBand(TH1F*& nom,
85  std::map<std::string,TH1F*>& ups,
86  std::map<std::string,TH1F*>& dns,
87  int col, int errCol,
88  float headroom, bool newaxis)
89 {
90  if(col == -1){
91  col = kTotalMCColor;
92  errCol = kTotalMCErrorBandColor;
93  }
94  else if(errCol == -1) errCol = col-7; // hopefully a lighter version
95 
96  TCanvas *c = new TCanvas("c","c");
97 
98  nom->SetLineColor(col);
99  nom->GetXaxis()->CenterTitle();
100  nom->GetYaxis()->CenterTitle();
101  if(newaxis) nom->Draw("hist ]["); // Set the axes up
102 
103  double yMax = nom->GetBinContent(nom->GetMaximumBin());
104 
105  TGraphAsymmErrors* g = new TGraphAsymmErrors;
106  return g; // TODO: handle shape and cherenkov, update maps and vectors
107 }
108 
109 //--------------------------------------
110 
111 TH1F* GetSelectedStatisticalUncertainty(Spectrum spectrum_nominal, double pot_function)
112 {
113  TH1F* nominal = (TH1F*)spectrum_nominal.ToTH1(pot_function);
114 
115  //Now get the statistical error
116  TH1F* stat_err_cut_value = (TH1F*)nominal->Clone("error_cut_value");
117 
118  for(int ibin = 1; ibin <= nominal->GetXaxis()->GetNbins(); ibin++){
119  float fSelectedErr = nominal->GetBinContent(ibin);
120  stat_err_cut_value->SetBinContent(ibin,sqrt(fSelectedErr));
121  }
122 
123  TCanvas *c2 = new TCanvas("c2","c2");
124  stat_err_cut_value->Draw("hist");
125  std::string outname = "UncertaintyPerBinPlots/debug/statistical_error_selected_vs_lower_cut.pdf";
126  c2->SaveAs(outname.c_str());
127  std::string outname_root = "UncertaintyPerBinPlots/root/statistical_error_selected_vs_lower_cut.root";
128  c2->SaveAs(outname_root.c_str());
129  c2->Close();
130  delete c2;
131 
132  TH1F* fractional_uncertainty = (TH1F*)stat_err_cut_value->Clone("h1");
133  fractional_uncertainty->Divide(nominal);
134  TCanvas *c3 = new TCanvas("c3","c3");
135  fractional_uncertainty->Draw("hist");
136  outname = "UncertaintyPerBinPlots/debug/relative_statistical_error_selected_vs_lower_cut.pdf";
137  c3->SaveAs(outname.c_str());
138  outname_root = "UncertaintyPerBinPlots/root/relative_statistical_error_selected_vs_lower_cut.root";
139  c3->SaveAs(outname_root.c_str());
140  c3->Close();
141  delete c3;
142 
143  return stat_err_cut_value;
144 
145 }
146 
147 TH1F* GetDenominator( Spectrum spectrum_selected,
148  Spectrum spectrum_background,
149  double pot_function)
150 {
151  TH1F* selected = (TH1F*)spectrum_selected.ToTH1(pot_function);
152  TH1F* background = (TH1F*)spectrum_background.ToTH1(pot_function);
153 
154  TH1F* denominator_cut_value = new TH1F("denom_cut_val",
155  selected->GetXaxis()->GetTitle(),
156  selected->GetXaxis()->GetNbins(),
157  selected->GetXaxis()->GetXmin(),
158  selected->GetXaxis()->GetXmax());
159 
160 
161  //Integrate To Get Total Number of Events selected for a particular cut value
162  for(int ibin = 1; ibin <= selected->GetXaxis()->GetNbins(); ibin++){
163  float fNumSelected = selected->GetBinContent(ibin);
164  float fNumBackground = background->GetBinContent(ibin);
165  denominator_cut_value->SetBinContent(ibin,fNumSelected - fNumBackground);
166  }
167 
168  TCanvas *c1 = new TCanvas("c1","c1");
169  denominator_cut_value->GetYaxis()->SetTitle("N_{sel}-N_{bkg}");
170  denominator_cut_value->GetXaxis()->SetTitle(selected->GetXaxis()->GetTitle());
171  denominator_cut_value->Draw("hist");
172  std::string outname = "UncertaintyPerBinPlots/debug/denominator_vs_lower_cut.pdf";
173  c1->SaveAs(outname.c_str());
174  std::string outname_root = "UncertaintyPerBinPlots/root/denominator_vs_lower_cut.root";
175  c1->SaveAs(outname_root.c_str());
176  c1->Close();
177  delete c1;
178 
179  return denominator_cut_value;
180 }
181 
182 TH1F* GetBackgroundStatisticalUncertainty(Spectrum spectrum_nominal, double pot_function)
183 {
184  TH1F* nominal = (TH1F*)spectrum_nominal.ToTH1(pot_function);
185 
186  //Now get the statistical error
187  TH1F* stat_err_cut_value = (TH1F*)nominal->Clone("error_cut_value_bkg");
188 
189  for(int ibin = 1; ibin <= nominal->GetXaxis()->GetNbins(); ibin++){
190  float fSelectedErr = nominal->GetBinContent(ibin);
191  stat_err_cut_value->SetBinContent(ibin,sqrt(fSelectedErr*pot/mc_pot));
192  }
193 
194  TCanvas *c2 = new TCanvas("c2","c2");
195  stat_err_cut_value->GetYaxis()->SetTitle("#delta N_{bkg}^{stat}");
196  stat_err_cut_value->Draw("hist");
197  std::string outname = "UncertaintyPerBinPlots/debug/statistical_error_background_vs_lower_cut.pdf";
198  c2->SaveAs(outname.c_str());
199  std::string outname_root = "UncertaintyPerBinPlots/root/statistical_error_background_vs_lower_cut.root";
200  c2->SaveAs(outname_root.c_str());
201  c2->Close();
202  delete c2;
203 
204  TH1F* fractional_uncertainty = (TH1F*)stat_err_cut_value->Clone("h1_bkg");
205  fractional_uncertainty->Divide(nominal);
206  TCanvas *c3 = new TCanvas("c3","c3");
207  fractional_uncertainty->GetYaxis()->SetTitle("#frac{#delta N_{bkg}^{stat}}{N_{bkg}}");
208  fractional_uncertainty->Draw("hist");
209  outname = "UncertaintyPerBinPlots/relative_statistical_error_background_vs_lower_cut.pdf";
210  c3->SaveAs(outname.c_str());
211  outname_root = "UncertaintyPerBinPlots/root/relative_statistical_error_background_vs_lower_cut.root";
212  c3->SaveAs(outname_root.c_str());
213  c3->Close();
214  delete c3;
215 
216  return stat_err_cut_value;
217 }
218 
220  TH1F* denominator,
221  std::map<std::string,Spectrum> ShiftUp,
222  std::map<std::string,Spectrum> ShiftDown,
223  double pot_function)
224 {
225  TH1F* nominal = (TH1F*)spectrum_nominal.ToTH1(pot_function);
226  std::map<std::string,TH1F*> hist_shift_up;
227  std::map<std::string,TH1F*> hist_shift_down;
228 
229  for(auto it = ShiftUp.begin(); it != ShiftUp.end(); ++it)
230  hist_shift_up.insert(std::make_pair(it->first,(TH1F*)ShiftUp.at(it->first).ToTH1(pot_function)));
231 
232  for(auto it = ShiftDown.begin(); it != ShiftDown.end(); ++it)
233  hist_shift_down.insert(std::make_pair(it->first,(TH1F*)ShiftDown.at(it->first).ToTH1(pot_function)));
234 
235  //Need histograms to hold all intermediate values
236  //currently there are 5 systematics shifts in both up and down vectors
237  TH1F* nominal_cut_value = (TH1F*)nominal->Clone("nominal_cut_value_syst");
238  TH1F* genie_up_cut_value = (TH1F*)nominal->Clone("genie_up_cut_value");
239  TH1F* genie_down_cut_value = (TH1F*)nominal->Clone("genie_down_cut_value");
240  TH1F* cal_up_cut_value = (TH1F*)nominal->Clone("cal_up_cut_value");
241  TH1F* cal_down_cut_value = (TH1F*)nominal->Clone("cal_down_cut_value");
242  TH1F* light_up_cut_value = (TH1F*)nominal->Clone("light_up_cut_value");
243  TH1F* light_down_cut_value = (TH1F*)nominal->Clone("light_down_cut_value");
244  TH1F* shape_cut_value = (TH1F*)nominal->Clone("shape_cut_value");
245  TH1F* cheren_cut_value = (TH1F*)nominal->Clone("cheren_cut_value");
246 
247  for(int ibin = 1; ibin <= nominal->GetXaxis()->GetNbins(); ibin++){
248  float fNumSelected = nominal->GetBinContent(ibin);
249  nominal_cut_value->SetBinContent(ibin,fNumSelected);
250 
251  float fGenieUp = hist_shift_up.at("genie")->GetBinContent(ibin);
252  float fGenieDown = hist_shift_down.at("genie")->GetBinContent(ibin);
253  genie_up_cut_value->SetBinContent(ibin, fGenieUp);
254  genie_down_cut_value->SetBinContent(ibin,fGenieDown);
255 
256  float fCalUp = hist_shift_up.at("calib")->GetBinContent(ibin);
257  float fCalDown = hist_shift_down.at("calib")->GetBinContent(ibin);
258  cal_up_cut_value->SetBinContent(ibin, fCalUp);
259  cal_down_cut_value->SetBinContent(ibin,fCalDown);
260 
261  float fLightUp = hist_shift_up.at("light")->GetBinContent(ibin);
262  float fLightDown = hist_shift_down.at("light")->GetBinContent(ibin);
263  light_up_cut_value->SetBinContent(ibin, fLightUp);
264  light_down_cut_value->SetBinContent(ibin,fLightDown);
265 
266  float fShape = hist_shift_up.at("shape")->GetBinContent(ibin);
267  shape_cut_value->SetBinContent(ibin, fShape);
268 
269  float fCheren = hist_shift_up.at("cherenkov")->GetBinContent(ibin);
270  cheren_cut_value->SetBinContent(ibin, fCheren);
271  }
272 
273  TCanvas *c_syst_cut = new TCanvas("c_syst_cut","c_syst_cut");
274  nominal_cut_value->SetLineColor(kBlack);
275 
276  genie_up_cut_value->SetLineColor(kRed);
277  cal_up_cut_value->SetLineColor(kGreen);
278  light_up_cut_value->SetLineColor(kOrange);
279  shape_cut_value->SetLineColor(kMagenta);
280  cheren_cut_value->SetLineColor(kBlue);
281 
282  genie_down_cut_value->SetLineColor(kRed+1);
283  cal_down_cut_value->SetLineColor(kGreen+1);
284  light_down_cut_value->SetLineColor(kOrange+1);
285 
286  TLegend *leg_cut = new TLegend(0.65,0.7,0.9,0.9);
287  leg_cut->AddEntry(nominal_cut_value,"Nominal");
288  leg_cut->AddEntry(genie_up_cut_value,"Genie up");
289  leg_cut->AddEntry(genie_down_cut_value,"Genie down");
290  leg_cut->AddEntry(cal_up_cut_value,"Cal up");
291  leg_cut->AddEntry(cal_down_cut_value,"Cal down");
292  leg_cut->AddEntry(light_up_cut_value,"Light up");
293  leg_cut->AddEntry(light_down_cut_value,"Light down");
294  leg_cut->AddEntry(shape_cut_value,"Shape");
295  leg_cut->AddEntry(cheren_cut_value,"Cherenkov");
296 
297  genie_up_cut_value->Draw("");
298  genie_down_cut_value->Draw("same");
299  cal_up_cut_value->Draw("same");
300  cal_down_cut_value->Draw("same");
301  light_up_cut_value->Draw("same");
302  light_down_cut_value->Draw("same");
303  shape_cut_value->Draw("same");
304  cheren_cut_value->Draw("same");
305  nominal_cut_value->Draw("same");
306  leg_cut->Draw("same");
307  c_syst_cut->SaveAs("UncertaintyPerBinPlots/bkg_syst_cut_comparison.pdf");
308  c_syst_cut->SaveAs("UncertaintyPerBinPlots/root/bkg_syst_cut_comparison.root");
309  c_syst_cut->Close();
310  delete c_syst_cut;
311 
312  std::map<std::string,TH1F*> shifts_up;
313  shifts_up.insert(std::make_pair("genie",genie_up_cut_value));
314  shifts_up.insert(std::make_pair("calib",cal_up_cut_value));
315  shifts_up.insert(std::make_pair("light",light_up_cut_value));
316  shifts_up.insert(std::make_pair("shape",shape_cut_value));
317  shifts_up.insert(std::make_pair("cherenkov",cheren_cut_value));
318 
319  std::map<std::string,TH1F*> shifts_down;
320  shifts_down.insert(std::make_pair("genie",genie_down_cut_value));
321  shifts_down.insert(std::make_pair("calib",cal_down_cut_value));
322  shifts_down.insert(std::make_pair("light",light_down_cut_value));
323 
324  new TCanvas();
325  MyPlotWithSystErrorBand(nominal_cut_value, shifts_up, shifts_down, -1,
326  -1, 1.3, true);
327  std::string outname = "UncertaintyPerBinPlots/debug/background_with_syst_errors_vs_lower_cut.pdf";
328  gPad->SaveAs(outname.c_str());
329  std::string outname_root = "UncertaintyPerBinPlots/root/background_with_syst_errors_vs_lower_cut.root";
330  gPad->SaveAs(outname_root.c_str());
331 
332  TH1F* syst_err_cut_value_tot = (TH1F*)nominal_cut_value->Clone("syst_err_cut_value_tot");
333  TH1F* syst_err_cut_value_genie = (TH1F*)nominal_cut_value->Clone("syst_err_cut_value_genie");
334  TH1F* syst_err_cut_value_calib = (TH1F*)nominal_cut_value->Clone("syst_err_cut_value_calib");
335  TH1F* syst_err_cut_value_light = (TH1F*)nominal_cut_value->Clone("syst_err_cut_value_light");
336  TH1F* syst_err_cut_value_shape = (TH1F*)nominal_cut_value->Clone("syst_err_cut_value_shape");
337  TH1F* syst_err_cut_value_cheren = (TH1F*)nominal_cut_value->Clone("syst_err_cut_value_cheren");
338  TH1F* syst_err_cut_value_genie_up = (TH1F*)nominal_cut_value->Clone("syst_err_cut_value_genie_up");
339  TH1F* syst_err_cut_value_calib_up = (TH1F*)nominal_cut_value->Clone("syst_err_cut_value_calib_up");
340  TH1F* syst_err_cut_value_light_up = (TH1F*)nominal_cut_value->Clone("syst_err_cut_value_light_up");
341  TH1F* syst_err_cut_value_genie_down = (TH1F*)nominal_cut_value->Clone("syst_err_cut_value_genie_down");
342  TH1F* syst_err_cut_value_calib_down = (TH1F*)nominal_cut_value->Clone("syst_err_cut_value_calib_down");
343  TH1F* syst_err_cut_value_light_down = (TH1F*)nominal_cut_value->Clone("syst_err_cut_value_light_down");
344 
345  if(up_down_combined) {
346  for(int ibin = 1; ibin <= nominal->GetXaxis()->GetNbins(); ibin++){
347 
348  const double i_nominal = nominal_cut_value->GetBinContent(ibin);
349 
350  double i_bkg_syst_uncert = 0;
351  for(auto it_syst = shifts_up.begin(); it_syst != shifts_up.end(); ++it_syst){
352 
353  //if(it_syst->first == "genie")continue;
354  double isyst_bkg_syst_uncert=0;
355  if(it_syst->first == "shape" || it_syst->first == "cherenkov") // 3==shape or 4==cherenkov
356  isyst_bkg_syst_uncert = shifts_up.at(it_syst->first)->GetBinContent(ibin) - i_nominal;
357  else
358  isyst_bkg_syst_uncert = (shifts_up.at(it_syst->first)->GetBinContent(ibin) -
359  shifts_down.at(it_syst->first)->GetBinContent(ibin)) / 2;
360 
361  i_bkg_syst_uncert += isyst_bkg_syst_uncert*isyst_bkg_syst_uncert;
362  //Debug syst contributions
363  if(it_syst->first == "genie")
364  syst_err_cut_value_genie->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
365  if(it_syst->first == "calib")
366  syst_err_cut_value_calib->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
367  if(it_syst->first == "light")
368  syst_err_cut_value_light->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
369  if(it_syst->first == "shape")
370  syst_err_cut_value_shape->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
371  if(it_syst->first == "cherenkov")
372  syst_err_cut_value_cheren->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
373  } // end for systIdx
374 
375  syst_err_cut_value_tot->SetBinContent(ibin, sqrt(std::abs(i_bkg_syst_uncert)));
376  }
377 
378  TCanvas *c2 = new TCanvas("c2","c2");
379  syst_err_cut_value_tot->Draw("hist");
380  outname = "UncertaintyPerBinPlots/systematic_error_background_vs_lower_cut.pdf";
381  c2->SaveAs(outname.c_str());
382  outname_root = "UncertaintyPerBinPlots/root/systematic_error_background_vs_lower_cut.root";
383  c2->SaveAs(outname_root.c_str());
384  c2->Close();
385  delete c2;
386 
387  TCanvas *c_bkg_syst_decomp = new TCanvas("c_bkg_syst_decomp","c_bkg_syst_decomp");
388  syst_err_cut_value_tot->GetYaxis()->SetTitle("#delta N_{bkg}");
389  syst_err_cut_value_tot->SetLineColor(kBlack);
390  syst_err_cut_value_genie->SetLineColor(kRed);
391  syst_err_cut_value_calib->SetLineColor(kGreen);
392  syst_err_cut_value_light->SetLineColor(kOrange);
393  syst_err_cut_value_shape->SetLineColor(kMagenta);
394  syst_err_cut_value_cheren->SetLineColor(kBlue);
395 
396  TLegend *leg_eff_decomp = new TLegend(0.65,0.7,0.9,0.9);
397  leg_eff_decomp->AddEntry(syst_err_cut_value_tot,"Total");
398  leg_eff_decomp->AddEntry(syst_err_cut_value_genie,"Genie");
399  leg_eff_decomp->AddEntry(syst_err_cut_value_calib,"Calib");
400  leg_eff_decomp->AddEntry(syst_err_cut_value_light,"Light");
401  leg_eff_decomp->AddEntry(syst_err_cut_value_shape,"Shape");
402  leg_eff_decomp->AddEntry(syst_err_cut_value_cheren,"Cherenkov");
403 
404  syst_err_cut_value_tot->Draw("hist");
405  syst_err_cut_value_genie->Draw("histsame");
406  syst_err_cut_value_calib->Draw("histsame");
407  syst_err_cut_value_light->Draw("histsame");
408  syst_err_cut_value_shape->Draw("histsame");
409  syst_err_cut_value_cheren->Draw("histsame");
410  leg_eff_decomp->Draw("same");
411  c_bkg_syst_decomp->SaveAs("UncertaintyPerBinPlots/delta_bkg_syst_decomp.pdf");
412  c_bkg_syst_decomp->SaveAs("UncertaintyPerBinPlots/root/delta_bkg_syst_decomp.root");
413 
414  TH1* copy_syst_err_cut_value_tot = (TH1F*)syst_err_cut_value_tot->Clone("copy_syst_err_cut_value_tot");
415  c_bkg_syst_decomp->SetLeftMargin(0.15);
416  copy_syst_err_cut_value_tot->GetYaxis()->SetTitle("#frac{#delta N_{bkg}}{N_{sel}-N_{bkg}}");
417  copy_syst_err_cut_value_tot->GetYaxis()->SetRangeUser(0,1);
418  copy_syst_err_cut_value_tot->Divide(denominator);
419  syst_err_cut_value_genie->Divide(denominator);
420  syst_err_cut_value_calib->Divide(denominator);
421  syst_err_cut_value_light->Divide(denominator);
422  syst_err_cut_value_shape->Divide(denominator);
423  syst_err_cut_value_cheren->Divide(denominator);
424 
425  copy_syst_err_cut_value_tot->Draw("hist");
426  syst_err_cut_value_genie->Draw("histsame");
427  syst_err_cut_value_calib->Draw("histsame");
428  syst_err_cut_value_light->Draw("histsame");
429  syst_err_cut_value_shape->Draw("histsame");
430  syst_err_cut_value_cheren->Draw("histsame");
431  leg_eff_decomp->Draw("same");
432  c_bkg_syst_decomp->SaveAs("UncertaintyPerBinPlots/bkg_syst_decomp.pdf");
433  c_bkg_syst_decomp->SaveAs("UncertaintyPerBinPlots/root/bkg_syst_decomp.root");
434 
435  c_bkg_syst_decomp->Close();
436  delete c_bkg_syst_decomp;
437  }
438  else {
439  for(int ibin = 1; ibin <= nominal->GetXaxis()->GetNbins(); ibin++){
440 
441  const double i_nominal = nominal_cut_value->GetBinContent(ibin);
442 
443  double i_bkg_syst_uncert = 0;
444 
445  for(auto it_syst = shifts_up.begin(); it_syst != shifts_up.end(); ++it_syst){
446 
447  double isyst_bkg_syst_uncert=0;
448 
449  isyst_bkg_syst_uncert = shifts_up.at(it_syst->first)->GetBinContent(ibin) - i_nominal;
450 
451  i_bkg_syst_uncert += isyst_bkg_syst_uncert*isyst_bkg_syst_uncert;
452 
453  //Debug syst contributions
454  if(it_syst->first == "genie")
455  syst_err_cut_value_genie_up->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
456  if(it_syst->first == "calib")
457  syst_err_cut_value_calib_up->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
458  if(it_syst->first == "light")
459  syst_err_cut_value_light_up->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
460  if(it_syst->first == "shape")
461  syst_err_cut_value_shape->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
462  if(it_syst->first == "cherenkov")
463  syst_err_cut_value_cheren->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
464  } // end for systIdx
465 
466  for(auto it_syst = shifts_down.begin(); it_syst != shifts_down.end(); ++it_syst){
467 
468  double isyst_bkg_syst_uncert=0;
469 
470  isyst_bkg_syst_uncert = shifts_down.at(it_syst->first)->GetBinContent(ibin) - i_nominal;
471 
472  i_bkg_syst_uncert += isyst_bkg_syst_uncert*isyst_bkg_syst_uncert;
473 
474  //Debug syst contributions
475  if(it_syst->first == "genie")
476  syst_err_cut_value_genie_down->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
477  if(it_syst->first == "calib")
478  syst_err_cut_value_calib_down->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
479  if(it_syst->first == "light")
480  syst_err_cut_value_light_down->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
481  } // end for systIdx
482 
483  syst_err_cut_value_tot->SetBinContent(ibin, sqrt(std::abs(i_bkg_syst_uncert)));
484  }
485 
486  TCanvas *c2 = new TCanvas("c2","c2");
487  syst_err_cut_value_tot->Draw("hist");
488  outname = "UncertaintyPerBinPlots/systematic_error_background_vs_lower_cut.pdf";
489  c2->SaveAs(outname.c_str());
490  outname_root = "UncertaintyPerBinPlots/root/systematic_error_background_vs_lower_cut.root";
491  c2->SaveAs(outname_root.c_str());
492  c2->Close();
493  delete c2;
494 
495  TCanvas *c_bkg_syst_decomp = new TCanvas("c_bkg_syst_decomp","c_bkg_syst_decomp");
496  syst_err_cut_value_tot->GetYaxis()->SetTitle("#delta N_{bkg}");
497  syst_err_cut_value_tot->SetLineColor(kBlack);
498  syst_err_cut_value_genie_up->SetLineColor(kRed);
499  syst_err_cut_value_calib_up->SetLineColor(kGreen);
500  syst_err_cut_value_light_up->SetLineColor(kOrange);
501  syst_err_cut_value_genie_down->SetLineColor(kRed+1);
502  syst_err_cut_value_calib_down->SetLineColor(kGreen+1);
503  syst_err_cut_value_light_down->SetLineColor(kOrange+1);
504  syst_err_cut_value_shape->SetLineColor(kMagenta);
505  syst_err_cut_value_cheren->SetLineColor(kBlue);
506 
507  TLegend *leg_eff_decomp = new TLegend(0.65,0.7,0.9,0.9);
508  leg_eff_decomp->AddEntry(syst_err_cut_value_tot,"Total");
509  leg_eff_decomp->AddEntry(syst_err_cut_value_genie_up,"Genie up");
510  leg_eff_decomp->AddEntry(syst_err_cut_value_calib_up,"Cal up");
511  leg_eff_decomp->AddEntry(syst_err_cut_value_light_up,"Light up");
512  leg_eff_decomp->AddEntry(syst_err_cut_value_genie_down,"Genie down");
513  leg_eff_decomp->AddEntry(syst_err_cut_value_calib_down,"Cal down");
514  leg_eff_decomp->AddEntry(syst_err_cut_value_light_down,"Light down");
515  leg_eff_decomp->AddEntry(syst_err_cut_value_shape,"Shape");
516  leg_eff_decomp->AddEntry(syst_err_cut_value_cheren,"Cherenkov");
517 
518  syst_err_cut_value_tot->Draw("hist");
519  syst_err_cut_value_genie_up->Draw("histsame");
520  syst_err_cut_value_calib_up->Draw("histsame");
521  syst_err_cut_value_light_up->Draw("histsame");
522  syst_err_cut_value_genie_down->Draw("histsame");
523  syst_err_cut_value_calib_down->Draw("histsame");
524  syst_err_cut_value_light_down->Draw("histsame");
525  syst_err_cut_value_shape->Draw("histsame");
526  syst_err_cut_value_cheren->Draw("histsame");
527  leg_eff_decomp->Draw("same");
528  c_bkg_syst_decomp->SaveAs("UncertaintyPerBinPlots/delta_bkg_syst_decomp.pdf");
529  c_bkg_syst_decomp->SaveAs("UncertaintyPerBinPlots/root/delta_bkg_syst_decomp.root");
530 
531  TH1* copy_syst_err_cut_value_tot = (TH1F*)syst_err_cut_value_tot->Clone("copy_syst_err_cut_value_tot");
532  c_bkg_syst_decomp->SetLeftMargin(0.15);
533  copy_syst_err_cut_value_tot->GetYaxis()->SetTitle("#frac{#delta N_{bkg}}{N_{sel}-N_{bkg}}");
534  copy_syst_err_cut_value_tot->GetYaxis()->SetRangeUser(0,1);
535  copy_syst_err_cut_value_tot->Divide(denominator);
536  syst_err_cut_value_genie_up->Divide(denominator);
537  syst_err_cut_value_calib_up->Divide(denominator);
538  syst_err_cut_value_light_up->Divide(denominator);
539  syst_err_cut_value_genie_down->Divide(denominator);
540  syst_err_cut_value_calib_down->Divide(denominator);
541  syst_err_cut_value_light_down->Divide(denominator);
542  syst_err_cut_value_shape->Divide(denominator);
543  syst_err_cut_value_cheren->Divide(denominator);
544 
545  copy_syst_err_cut_value_tot->Draw("hist");
546  syst_err_cut_value_genie_up->Draw("histsame");
547  syst_err_cut_value_calib_up->Draw("histsame");
548  syst_err_cut_value_light_up->Draw("histsame");
549  syst_err_cut_value_genie_down->Draw("histsame");
550  syst_err_cut_value_calib_down->Draw("histsame");
551  syst_err_cut_value_light_down->Draw("histsame");
552  syst_err_cut_value_shape->Draw("histsame");
553  syst_err_cut_value_cheren->Draw("histsame");
554  leg_eff_decomp->Draw("same");
555  c_bkg_syst_decomp->SaveAs("UncertaintyPerBinPlots/bkg_syst_decomp.pdf");
556  c_bkg_syst_decomp->SaveAs("UncertaintyPerBinPlots/root/bkg_syst_decomp.root");
557 
558  c_bkg_syst_decomp->Close();
559  delete c_bkg_syst_decomp;
560 
561  }
562  return syst_err_cut_value_tot;
563 }
564 
566  Spectrum spec_truth_nominal,
567  std::map<std::string,Spectrum> ShiftUp,
568  std::map<std::string,Spectrum> ShiftDown,
569  std::map<std::string,Spectrum> TruthUp,
570  std::map<std::string,Spectrum> TruthDown,
571  double pot_function)
572 {
573  TH1F* signal_nominal = (TH1F*)spec_signal_nominal.ToTH1(pot_function);
574  TH1F* truth_nominal = (TH1F*)spec_truth_nominal.ToTH1(pot_function);
575 
576  std::map<std::string,TH1F*> hist_shift_up;
577  std::map<std::string,TH1F*> hist_shift_down;
578  std::map<std::string,TH1F*> truth_shift_up;
579  std::map<std::string,TH1F*> truth_shift_down;
580 
581  for(auto it = ShiftUp.begin(); it != ShiftUp.end(); ++it){
582  hist_shift_up.insert(std::make_pair(it->first,(TH1F*)ShiftUp.at(it->first).ToTH1(pot_function)));
583  truth_shift_up.insert(std::make_pair(it->first,(TH1F*)TruthUp.at(it->first).ToTH1(pot_function)));
584  }
585  for(auto it = ShiftDown.begin(); it != ShiftDown.end(); ++it) {
586  hist_shift_down.insert(std::make_pair(it->first,(TH1F*)ShiftDown.at(it->first).ToTH1(pot_function)));
587  truth_shift_down.insert(std::make_pair(it->first,(TH1F*)TruthDown.at(it->first).ToTH1(pot_function)));
588  }
589 
590  //Need histograms to hold all intermediate values
591  //currently there are 5 systematics shifts in both up and down vectors
592  //TH1* copy_syst_err_cut_value_eff_tot = (TH1F*)syst_err_cut_value_eff_tot->Clone("copy_syst_err_cut_value_eff_tot");
593 
594  TH1F* nominal_cut_value_eff = (TH1F*)signal_nominal->Clone("nominal_cut_value_eff");
595  TH1F* genie_up_cut_value_eff = (TH1F*)hist_shift_up.at("genie")->Clone("genie_up_cut_value_eff");
596  TH1F* genie_down_cut_value_eff = (TH1F*)hist_shift_down.at("genie")->Clone("genie_down_cut_value_eff");
597  TH1F* cal_up_cut_value_eff = (TH1F*)hist_shift_up.at("calib")->Clone("cal_up_cut_value_eff");
598  TH1F* cal_down_cut_value_eff = (TH1F*)hist_shift_down.at("calib")->Clone("cal_down_cut_value_eff");
599  TH1F* light_up_cut_value_eff = (TH1F*)hist_shift_up.at("light")->Clone("light_up_cut_value_eff");
600  TH1F* light_down_cut_value_eff = (TH1F*)hist_shift_down.at("light")->Clone("light_down_cut_value_eff");
601  TH1F* shape_cut_value_eff = (TH1F*)hist_shift_up.at("shape")->Clone("shape_cut_value_eff");
602  TH1F* cheren_cut_value_eff = (TH1F*)hist_shift_up.at("cherenkov")->Clone("cheren_cut_value_eff");
603 
604  TCanvas *c = new TCanvas("c","c");
605  TH1F* truth_genie_up = (TH1F*)truth_shift_up.at("genie")->Clone("truth_genie_up");
606  TH1F* truth_genie_down = (TH1F*)truth_shift_down.at("genie")->Clone("truth_genie_down");
607  TH1F* truth_cal_up = (TH1F*)truth_shift_up.at("calib")->Clone("truth_cal_up");
608  TH1F* truth_cal_down = (TH1F*)truth_shift_down.at("calib")->Clone("truth_cal_down");
609  TH1F* truth_light_up = (TH1F*)truth_shift_up.at("light")->Clone("truth_light_up");
610  TH1F* truth_light_down = (TH1F*)truth_shift_down.at("light")->Clone("truth_light_down");
611  TH1F* truth_shape = (TH1F*)truth_shift_up.at("shape")->Clone("truth_shape");
612  TH1F* truth_cheren = (TH1F*)truth_shift_up.at("cherenkov")->Clone("truth_cheren");
613  truth_nominal->SetLineColor(kBlack);
614  truth_genie_up->SetLineColor(kRed);
615  truth_cal_up->SetLineColor(kGreen);
616  truth_light_up->SetLineColor(kOrange);
617  truth_shape->SetLineColor(kMagenta);
618  truth_cheren->SetLineColor(kBlue);
619  truth_genie_down->SetLineColor(kRed+1);
620  truth_cal_down->SetLineColor(kGreen+1);
621  truth_light_down->SetLineColor(kOrange+1);
622  truth_genie_up->SetTitle("True Signal Events");
623  truth_genie_up->Draw("hist");
624  truth_genie_down->Draw("hist same");
625  truth_cal_up->Draw("hist same");
626  truth_cal_down->Draw("hist same");
627  truth_light_up->Draw("hist same");
628  truth_light_down->Draw("hist same");
629  truth_shape->Draw("hist same");
630  truth_cheren->Draw("hist same");
631  truth_nominal->Draw("hist same");
632  TLegend *leg1 = new TLegend(0.65,0.7,0.9,0.9);
633  leg1->AddEntry(truth_nominal,"Nominal");
634  leg1->AddEntry(truth_genie_up,"Genie up");
635  leg1->AddEntry(truth_genie_down,"Genie down");
636  leg1->AddEntry(truth_cal_up,"Cal up");
637  leg1->AddEntry(truth_cal_down,"Cal down");
638  leg1->AddEntry(truth_light_up,"Light up");
639  leg1->AddEntry(truth_light_down,"Light down");
640  leg1->AddEntry(truth_shape,"Shape up");
641  leg1->AddEntry(truth_cheren,"Cherenkov up");
642  leg1->Draw();
643  c->SaveAs("UncertaintyPerBinPlots/true_syst_comparison.pdf");
644  c->SaveAs("UncertaintyPerBinPlots/root/true_syst_comparison.root");
645  c->Close();
646  delete c;
647 
648 
649 
650  nominal_cut_value_eff->Divide(truth_nominal);
651  genie_up_cut_value_eff->Divide(truth_shift_up.at("genie"));
652  genie_down_cut_value_eff->Divide(truth_shift_down.at("genie"));
653  cal_up_cut_value_eff->Divide(truth_shift_up.at("calib"));
654  cal_down_cut_value_eff->Divide(truth_shift_down.at("calib"));
655  light_up_cut_value_eff->Divide(truth_shift_up.at("light"));
656  light_down_cut_value_eff->Divide(truth_shift_down.at("light"));
657  shape_cut_value_eff->Divide(truth_shift_up.at("shape"));
658  cheren_cut_value_eff->Divide(truth_shift_up.at("cherenkov"));
659 
660  TCanvas *c_eff = new TCanvas("c_eff","c_eff");
661 
662  nominal_cut_value_eff->SetLineColor(kBlack);
663 
664  genie_up_cut_value_eff->SetLineColor(kRed);
665  cal_up_cut_value_eff->SetLineColor(kGreen);
666  light_up_cut_value_eff->SetLineColor(kOrange);
667  shape_cut_value_eff->SetLineColor(kMagenta);
668  cheren_cut_value_eff->SetLineColor(kBlue);
669 
670  genie_down_cut_value_eff->SetLineColor(kRed+1);
671  cal_down_cut_value_eff->SetLineColor(kGreen+1);
672  light_down_cut_value_eff->SetLineColor(kOrange+1);
673 
674  TLegend *leg = new TLegend(0.65,0.7,0.9,0.9);
675  leg->AddEntry(nominal_cut_value_eff,"Nominal");
676  leg->AddEntry(genie_up_cut_value_eff,"Genie up");
677  leg->AddEntry(genie_down_cut_value_eff,"Genie down");
678  leg->AddEntry(cal_up_cut_value_eff,"Cal up");
679  leg->AddEntry(cal_down_cut_value_eff,"Cal down");
680  leg->AddEntry(light_up_cut_value_eff,"Light up");
681  leg->AddEntry(light_down_cut_value_eff,"Light down");
682  leg->AddEntry(shape_cut_value_eff,"Shape up");
683  leg->AddEntry(cheren_cut_value_eff,"Cherenkov up");
684 
685  //nominal_cut_value_eff->GetYaxis()->SetRangeUser(0,0.08);
686  cheren_cut_value_eff->GetYaxis()->SetTitle("#epsilon");
687  cheren_cut_value_eff->Draw("hist");
688  genie_down_cut_value_eff->Draw("hist same");
689  cal_up_cut_value_eff->Draw("hist same");
690  cal_down_cut_value_eff->Draw("hist same");
691  light_up_cut_value_eff->Draw("hist same");
692  light_down_cut_value_eff->Draw("hist same");
693  shape_cut_value_eff->Draw("hist same");
694  genie_up_cut_value_eff->Draw("hist same");
695  nominal_cut_value_eff->Draw("hist same");
696  leg->Draw("same");
697  c_eff->SaveAs("UncertaintyPerBinPlots/eff_syst_comparison.pdf");
698  c_eff->SaveAs("UncertaintyPerBinPlots/root/eff_syst_comparison.root");
699  c_eff->Close();
700  delete c_eff;
701 
702  std::map<std::string,TH1F*> shifts_up;
703  shifts_up.insert(std::make_pair("genie",genie_up_cut_value_eff));
704  shifts_up.insert(std::make_pair("calib",cal_up_cut_value_eff));
705  shifts_up.insert(std::make_pair("light",light_up_cut_value_eff));
706  shifts_up.insert(std::make_pair("shape",shape_cut_value_eff));
707  shifts_up.insert(std::make_pair("cherenkov",cheren_cut_value_eff));
708 
709  std::map<std::string,TH1F*> shifts_down;
710  shifts_down.insert(std::make_pair("genie",genie_down_cut_value_eff));
711  shifts_down.insert(std::make_pair("calib",cal_down_cut_value_eff));
712  shifts_down.insert(std::make_pair("light",light_down_cut_value_eff));
713 
714  new TCanvas();
715  MyPlotWithSystErrorBand(nominal_cut_value_eff, shifts_up,
716  shifts_down, -1,
717  -1, 1.3, true);
718  std::string outname = "UncertaintyPerBinPlots/debug/efficiency_with_syst_errors_vs_lower_cut.pdf";
719  gPad->SaveAs(outname.c_str());
720  std::string outname_root = "UncertaintyPerBinPlots/root/efficiency_with_syst_errors_vs_lower_cut.root";
721  gPad->SaveAs(outname_root.c_str());
722 
723  TH1F* syst_err_cut_value_eff_tot = (TH1F*)nominal_cut_value_eff->Clone("syst_err_cut_value_eff_tot");
724  TH1F* syst_err_cut_value_eff_truth = (TH1F*)nominal_cut_value_eff->Clone("syst_err_cut_value_eff_truth");
725  TH1F* syst_err_cut_value_eff_genie = (TH1F*)nominal_cut_value_eff->Clone("syst_err_cut_value_eff_genie");
726  TH1F* syst_err_cut_value_eff_calib = (TH1F*)nominal_cut_value_eff->Clone("syst_err_cut_value_eff_calib");
727  TH1F* syst_err_cut_value_eff_light = (TH1F*)nominal_cut_value_eff->Clone("syst_err_cut_value_eff_light");
728  TH1F* syst_err_cut_value_eff_shape = (TH1F*)nominal_cut_value_eff->Clone("syst_err_cut_value_eff_shape");
729  TH1F* syst_err_cut_value_eff_cheren = (TH1F*)nominal_cut_value_eff->Clone("syst_err_cut_value_eff_cheren");
730  TH1F* syst_err_cut_value_eff_genie_up = (TH1F*)nominal_cut_value_eff->Clone("syst_err_cut_value_eff_genie_up");
731  TH1F* syst_err_cut_value_eff_calib_up = (TH1F*)nominal_cut_value_eff->Clone("syst_err_cut_value_eff_calib_up");
732  TH1F* syst_err_cut_value_eff_light_up = (TH1F*)nominal_cut_value_eff->Clone("syst_err_cut_value_eff_light_up");
733  TH1F* syst_err_cut_value_eff_genie_down = (TH1F*)nominal_cut_value_eff->Clone("syst_err_cut_value_eff_genie_down");
734  TH1F* syst_err_cut_value_eff_calib_down = (TH1F*)nominal_cut_value_eff->Clone("syst_err_cut_value_eff_calib_down");
735  TH1F* syst_err_cut_value_eff_light_down = (TH1F*)nominal_cut_value_eff->Clone("syst_err_cut_value_eff_light_down");
736 
737 
738  if(up_down_combined) {
739  for(int ibin = 1; ibin <= signal_nominal->GetXaxis()->GetNbins(); ibin++){
740 
741  const double i_nominal = nominal_cut_value_eff->GetBinContent(ibin);
742 
743  double i_bkg_syst_uncert = 0;
744 
745  for(auto it_syst = shifts_up.begin(); it_syst != shifts_up.end(); ++it_syst){
746 
747  double isyst_bkg_syst_uncert=0;
748  if(it_syst->first == "genie")
749  syst_err_cut_value_eff_genie_up->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
750  if(it_syst->first == "calib")
751  syst_err_cut_value_eff_calib_up->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
752  if(it_syst->first == "light")
753  syst_err_cut_value_eff_light_up->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
754  if(it_syst->first == "shape")
755  syst_err_cut_value_eff_shape->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
756  if(it_syst->first == "cherenkov")
757  syst_err_cut_value_eff_cheren->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
758  } // end for syst loop
759 
760  for(auto it_syst = shifts_down.begin(); it_syst != shifts_down.end(); ++it_syst){
761 
762  double isyst_bkg_syst_uncert=0;
763 
764  isyst_bkg_syst_uncert = shifts_down.at(it_syst->first)->GetBinContent(ibin)-i_nominal;
765  //std::cout << " down " << it_syst->first << " " << isyst_bkg_syst_uncert << std::endl;
766  i_bkg_syst_uncert += isyst_bkg_syst_uncert*isyst_bkg_syst_uncert;
767 
768  //Debug syst contributions
769  if(it_syst->first == "genie")
770  syst_err_cut_value_eff_genie_down->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
771  if(it_syst->first == "calib")
772  syst_err_cut_value_eff_calib_down->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
773  if(it_syst->first == "light")
774  syst_err_cut_value_eff_light_down->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
775  } // end for syst loop
776 
777  syst_err_cut_value_eff_tot->SetBinContent(ibin, sqrt(std::abs(i_bkg_syst_uncert)));
778  }
779  TCanvas *c2 = new TCanvas("c2","c2");
780  syst_err_cut_value_eff_tot->Draw("hist");
781  outname = "UncertaintyPerBinPlots/debug/systematic_error_efficiency_vs_lower_cut.pdf";
782  c2->SaveAs(outname.c_str());
783  outname_root = "UncertaintyPerBinPlots/root/systematic_error_efficiency_vs_lower_cut.root";
784  c2->SaveAs(outname_root.c_str());
785  c2->Close();
786  delete c2;
787 
788  TCanvas *c_eff_decomp = new TCanvas("c_eff_decomp","c_eff_decomp");
789  syst_err_cut_value_eff_tot->GetYaxis()->SetTitle("#delta#epsilon");
790  syst_err_cut_value_eff_tot->SetLineColor(kBlack);
791  syst_err_cut_value_eff_genie_up->SetLineColor(kRed);
792  syst_err_cut_value_eff_calib_up->SetLineColor(kGreen);
793  syst_err_cut_value_eff_light_up->SetLineColor(kOrange);
794  syst_err_cut_value_eff_genie_down->SetLineColor(kRed+1);
795  syst_err_cut_value_eff_calib_down->SetLineColor(kGreen+1);
796  syst_err_cut_value_eff_light_down->SetLineColor(kOrange+1);
797  syst_err_cut_value_eff_shape->SetLineColor(kMagenta);
798  syst_err_cut_value_eff_cheren->SetLineColor(kBlue);
799 
800  TLegend *leg_eff_decomp = new TLegend(0.65,0.7,0.9,0.9);
801  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_tot,"Total");
802  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_genie_up,"Genie up");
803  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_calib_up,"Calib up");
804  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_light_up,"Light up");
805  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_genie_down,"Genie down");
806  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_calib_down,"Calib down");
807  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_light_down,"Light down");
808  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_shape,"Shape");
809  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_cheren,"Cherenkov");
810 
811  syst_err_cut_value_eff_tot->Draw("hist");
812  syst_err_cut_value_eff_genie_up->Draw("histsame");
813  syst_err_cut_value_eff_calib_up->Draw("histsame");
814  syst_err_cut_value_eff_light_up->Draw("histsame");
815  syst_err_cut_value_eff_genie_down->Draw("histsame");
816  syst_err_cut_value_eff_calib_down->Draw("histsame");
817  syst_err_cut_value_eff_light_down->Draw("histsame");
818  syst_err_cut_value_eff_shape->Draw("histsame");
819  syst_err_cut_value_eff_cheren->Draw("histsame");
820  leg_eff_decomp->Draw("same");
821  c_eff_decomp->SaveAs("UncertaintyPerBinPlots/efficiency_syst_error_numerator.pdf");
822  c_eff_decomp->SaveAs("UncertaintyPerBinPlots/root/efficiency_syst_error_numerator.root");
823 
824  TH1* copy_syst_err_cut_value_eff_tot = (TH1F*)syst_err_cut_value_eff_tot->Clone("copy_syst_err_cut_value_eff_tot");
825  copy_syst_err_cut_value_eff_tot->GetYaxis()->SetTitle("#frac{#delta#epsilon}{#epsilon}");
826  //copy_syst_err_cut_value_eff_tot->GetYaxis()->SetRangeUser(0,1);
827  copy_syst_err_cut_value_eff_tot->Divide(nominal_cut_value_eff);
828  syst_err_cut_value_eff_genie_up->Divide(nominal_cut_value_eff);
829  syst_err_cut_value_eff_calib_up->Divide(nominal_cut_value_eff);
830  syst_err_cut_value_eff_light_up->Divide(nominal_cut_value_eff);
831  syst_err_cut_value_eff_genie_down->Divide(nominal_cut_value_eff);
832  syst_err_cut_value_eff_calib_down->Divide(nominal_cut_value_eff);
833  syst_err_cut_value_eff_light_down->Divide(nominal_cut_value_eff);
834  syst_err_cut_value_eff_shape->Divide(nominal_cut_value_eff);
835  syst_err_cut_value_eff_cheren->Divide(nominal_cut_value_eff);
836 
837  copy_syst_err_cut_value_eff_tot->SetMaximum(0.4);
838  copy_syst_err_cut_value_eff_tot->Draw("hist");
839  syst_err_cut_value_eff_genie_up->Draw("histsame");
840  syst_err_cut_value_eff_calib_up->Draw("histsame");
841  syst_err_cut_value_eff_light_up->Draw("histsame");
842  syst_err_cut_value_eff_genie_down->Draw("histsame");
843  syst_err_cut_value_eff_calib_down->Draw("histsame");
844  syst_err_cut_value_eff_light_down->Draw("histsame");
845  syst_err_cut_value_eff_shape->Draw("histsame");
846  syst_err_cut_value_eff_cheren->Draw("histsame");
847  leg_eff_decomp->Draw("same");
848  c_eff_decomp->SaveAs("UncertaintyPerBinPlots/efficiency_syst_error_decomp.pdf");
849  c_eff_decomp->SaveAs("UncertaintyPerBinPlots/root/efficiency_syst_error_decomp.root");
850 
851  c_eff_decomp->Close();
852  delete c_eff_decomp;
853  }
854  return syst_err_cut_value_eff_tot;
855 }
856 
857 /////////////////////////////////////////////////////////
858 /////////////////////////////////////////////////////////
859 /////////////////////////////////////////////////////////
860 /////////////////////////////////////////////////////////
861 /////////////////////////////////////////////////////////
862 /////////////////////////////////////////////////////////
863 
864 
866 {
867  // std::string mkdir_name = "mkdir -p UncertaintyPerBinPlots_" + param + "/root";
868  // gSystem->Exec(mkdir_name.c_str());
869 
870  gSystem->Exec("rm -rf UncertaintyPerBinPlots");
871  gSystem->Exec("mkdir -p UncertaintyPerBinPlots/root");
872  gSystem->Exec("mkdir -p UncertaintyPerBinPlots/debug");
873 
874  std::string param_sel;
875  std::string param_obs_sig;
876  std::string param_obs_bkg;
877  std::string param_truth;
878 
879  bool left2right = false;
880 
881  if(param=="energy1") {
882  param_sel="RecoPionEnergySelAnalysisBinning";
883  param_obs_sig="RecoPionEnergySigAnalysisBinning";
884  param_obs_bkg="RecoPionEnergyBkgAnalysisBinning";
885  param_truth="RecoPionEnergyTruAnalysisBinning";
886  }
887  else if(param=="energy2") {
888  param_sel="RecoPionEnergySelAnalysisBinning1";
889  param_obs_sig="RecoPionEnergySigAnalysisBinning1";
890  param_obs_bkg="RecoPionEnergyBkgAnalysisBinning1";
891  param_truth="RecoPionEnergyTruAnalysisBinning1";
892  }
893  else if(param=="energy3") {
894  param_sel="RecoPionEnergySelAnalysisBinning2";
895  param_obs_sig="RecoPionEnergySigAnalysisBinning2";
896  param_obs_bkg="RecoPionEnergyBkgAnalysisBinning2";
897  param_truth="RecoPionEnergyTruAnalysisBinning2";
898  }
899  else{
900  cout<<"why"<<endl;
901  std::cout << "ERROR: Parameter to optimize against is not specified or not supported yet!" << std::endl;
902  return;
903  }
904 
905 
906  TFile *file1 = new TFile("/nova/app/users/projas/dev_19-01-22/ccpiinc_mc_studies-full-calibdown-10k-newKEbinning1.root", "read");
907  TFile *file2 = new TFile("/nova/app/users/projas/dev_19-01-22/ccpiinc_mc_studies-full-calibup-10k-newKEbinning1.root", "read");
908  TFile *file3 = new TFile("/nova/app/users/projas/dev_19-01-22/ccpiinc_mc_studies-full-ckvproton-10k-newKEbinning1.root", "read");
909  TFile *file4 = new TFile("/nova/app/users/projas/dev_19-01-22/ccpiinc_mc_studies-full-negoffset-10k-newKEbinning1.root", "read");
910  TFile *file5 = new TFile("/nova/app/users/projas/dev_19-01-22/ccpiinc_mc_studies-full-posoffset-10k-newKEbinning1.root", "read");
911  TFile *file6 = new TFile("/nova/app/users/projas/dev_19-01-22/ccpiinc_mc_studies-full-calibshape-10k-newKEbinning1.root");
912  TFile *filerw = new TFile("/nova/app/users/projas/dev_19-01-22/ccpiinc_mc_studies-full-geniemv-newKEbinning1.root", "read");
913  //*/
914 
915 
916  //O == Preselected 1 == Presel && Signal 2 == Presel & Bkgrd
917 
918  std::vector<GenieMultiverseSpectra> multiversespec;
919  std::vector<Spectrum> calibupspect;
920  std::vector<Spectrum> calibdownspect;
921  std::vector<Spectrum> lightupspect;
922  std::vector<Spectrum> lightdownspect;
923  std::vector<Spectrum> calibshapespect;
924  std::vector<Spectrum> cherenspect;
925 
926  std::vector<std::unique_ptr<ana::Spectrum>>genie_spects_sel;
927  std::vector<std::unique_ptr<ana::Spectrum>>genie_spects_obbkg;
928  std::vector<std::unique_ptr<ana::Spectrum>>genie_spects_obsig;
929  std::vector<std::unique_ptr<ana::Spectrum>>genie_spects_tru;
930  std::vector<std::unique_ptr<ana::Spectrum>>genie_spects_eff;
931 
932 
933  if(param=="energy1") {
934  for(int k = 0; k < 100; k++){
935  char name[50];
936  sprintf(name, "mc_geniemv_recoenergy1_%i_sel",k);
937  genie_spects_sel.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
938  sprintf(name, "mc_geniemv_recoenergy1_%i_sig",k);
939  genie_spects_obsig.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
940  auto tmpspecnum=Spectrum::LoadFrom(filerw->GetDirectory(name));
941  sprintf(name, "mc_geniemv_recoenergy1_%i_bkg",k);
942  genie_spects_obbkg.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
943  sprintf(name, "mc_geniemv_recoenergy1_%i_tru",k);
944  genie_spects_tru.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
945  auto tmpspecdenom=Spectrum::LoadFrom(filerw->GetDirectory(name));
946  TH1* tmpspecnumhist = tmpspecnum->ToTH1(mc_pot);
947  TH1* tmpspecdenomhist = tmpspecdenom->ToTH1(mc_pot);
948  TH1D* tmpeffhist = new TH1D("tmpeffhist",
949  TString(";") + tmpspecdenomhist->GetXaxis()->GetTitle(),
950  tmpspecdenomhist->GetXaxis()->GetNbins(),
951  tmpspecdenomhist->GetXaxis()->GetXmin(),
952  tmpspecdenomhist->GetXaxis()->GetXmax());
953  for(int ibin=1;ibin<tmpspecnumhist->GetXaxis()->GetNbins();ibin++)tmpeffhist->SetBinContent(ibin,tmpspecnumhist->Integral(ibin,tmpspecnumhist->GetXaxis()->GetNbins())/tmpspecdenomhist->Integral(ibin,tmpspecdenomhist->GetXaxis()->GetNbins()));
954  Spectrum* tmpeff = new Spectrum(tmpeffhist,tmpspecnum->POT(),tmpspecnum->Livetime());
955  //std::unique_ptr<Spectrum>tmpeff1 = std::make_unique<Spectrum>tmpeff;
956  genie_spects_eff.push_back((std::unique_ptr<Spectrum>)tmpeff);
957  delete tmpeffhist;
958  }
959  }
960 
961  if(param=="energy2") {
962  for(int k = 0; k < 100; k++){
963  char name[50];
964  sprintf(name, "mc_geniemv_recoenergy2_%i_sel",k);
965  genie_spects_sel.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
966  sprintf(name, "mc_geniemv_recoenergy2_%i_sig",k);
967  genie_spects_obsig.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
968  auto tmpspecnum=Spectrum::LoadFrom(filerw->GetDirectory(name));
969  sprintf(name, "mc_geniemv_recoenergy2_%i_bkg",k);
970  genie_spects_obbkg.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
971  sprintf(name, "mc_geniemv_recoenergy2_%i_tru",k);
972  genie_spects_tru.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
973  auto tmpspecdenom=Spectrum::LoadFrom(filerw->GetDirectory(name));
974  TH1* tmpspecnumhist = tmpspecnum->ToTH1(mc_pot);
975  TH1* tmpspecdenomhist = tmpspecdenom->ToTH1(mc_pot);
976  TH1D* tmpeffhist = new TH1D("tmpeffhist",
977  TString(";") + tmpspecdenomhist->GetXaxis()->GetTitle(),
978  tmpspecdenomhist->GetXaxis()->GetNbins(),
979  tmpspecdenomhist->GetXaxis()->GetXmin(),
980  tmpspecdenomhist->GetXaxis()->GetXmax());
981  for(int ibin=1;ibin<tmpspecnumhist->GetXaxis()->GetNbins();ibin++)tmpeffhist->SetBinContent(ibin,tmpspecnumhist->Integral(ibin,tmpspecnumhist->GetXaxis()->GetNbins())/tmpspecdenomhist->Integral(ibin,tmpspecdenomhist->GetXaxis()->GetNbins()));
982  Spectrum* tmpeff = new Spectrum(tmpeffhist,tmpspecnum->POT(),tmpspecnum->Livetime());
983  //std::unique_ptr<Spectrum>tmpeff1 = std::make_unique<Spectrum>tmpeff;
984  genie_spects_eff.push_back((std::unique_ptr<Spectrum>)tmpeff);
985  delete tmpeffhist;
986  }
987  }
988 
989  if(param=="energy3") {
990  for(int k = 0; k < 100; k++){
991  char name[50];
992  sprintf(name, "mc_geniemv_recoenergy3_%i_sel",k);
993  genie_spects_sel.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
994  sprintf(name, "mc_geniemv_recoenergy3_%i_sig",k);
995  genie_spects_obsig.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
996  auto tmpspecnum=Spectrum::LoadFrom(filerw->GetDirectory(name));
997  sprintf(name, "mc_geniemv_recoenergy3_%i_bkg",k);
998  genie_spects_obbkg.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
999  sprintf(name, "mc_geniemv_recoenergy3_%i_tru",k);
1000  genie_spects_tru.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
1001  auto tmpspecdenom=Spectrum::LoadFrom(filerw->GetDirectory(name));
1002  TH1* tmpspecnumhist = tmpspecnum->ToTH1(mc_pot);
1003  TH1* tmpspecdenomhist = tmpspecdenom->ToTH1(mc_pot);
1004  TH1D* tmpeffhist = new TH1D("tmpeffhist",
1005  TString(";") + tmpspecdenomhist->GetXaxis()->GetTitle(),
1006  tmpspecdenomhist->GetXaxis()->GetNbins(),
1007  tmpspecdenomhist->GetXaxis()->GetXmin(),
1008  tmpspecdenomhist->GetXaxis()->GetXmax());
1009  for(int ibin=1;ibin<tmpspecnumhist->GetXaxis()->GetNbins();ibin++)tmpeffhist->SetBinContent(ibin,tmpspecnumhist->Integral(ibin,tmpspecnumhist->GetXaxis()->GetNbins())/tmpspecdenomhist->Integral(ibin,tmpspecdenomhist->GetXaxis()->GetNbins()));
1010  Spectrum* tmpeff = new Spectrum(tmpeffhist,tmpspecnum->POT(),tmpspecnum->Livetime());
1011  //std::unique_ptr<Spectrum>tmpeff1 = std::make_unique<Spectrum>tmpeff;
1012  genie_spects_eff.push_back((std::unique_ptr<Spectrum>)tmpeff);
1013  delete tmpeffhist;
1014  }
1015  }
1016 
1017  GenieMultiverseSpectra sMulti_sel = GenieMultiverseSpectra(100,genie_spects_sel,true);
1018  GenieMultiverseSpectra sMulti_sig = GenieMultiverseSpectra(100,genie_spects_obsig,true);
1019  GenieMultiverseSpectra sMulti_bkg = GenieMultiverseSpectra(100,genie_spects_obbkg,true);
1020  GenieMultiverseSpectra sMulti_truth = GenieMultiverseSpectra(100,genie_spects_tru,true);
1021  GenieMultiverseSpectra sMulti_eff = GenieMultiverseSpectra(100,genie_spects_eff,true);
1022 
1023  multiversespec.push_back(sMulti_sel);
1024  multiversespec.push_back(sMulti_sig);
1025  multiversespec.push_back(sMulti_bkg);
1026  multiversespec.push_back(sMulti_eff);
1027 
1028  calibupspect.push_back(*Spectrum::LoadFrom(file4->GetDirectory(param_sel.c_str())));
1029  calibdownspect.push_back(*Spectrum::LoadFrom(file5->GetDirectory(param_sel.c_str())));
1030  lightupspect.push_back(*Spectrum::LoadFrom(file1->GetDirectory(param_sel.c_str())));
1031  lightdownspect.push_back(*Spectrum::LoadFrom(file2->GetDirectory(param_sel.c_str())));
1032  calibshapespect.push_back(*Spectrum::LoadFrom(file6->GetDirectory(param_sel.c_str())));
1033  cherenspect.push_back(*Spectrum::LoadFrom(file3->GetDirectory(param_sel.c_str())));
1034 
1035  calibupspect.push_back(*Spectrum::LoadFrom(file4->GetDirectory(param_obs_sig.c_str())));
1036  calibdownspect.push_back(*Spectrum::LoadFrom(file5->GetDirectory(param_obs_sig.c_str())));
1037  lightupspect.push_back(*Spectrum::LoadFrom(file1->GetDirectory(param_obs_sig.c_str())));
1038  lightdownspect.push_back(*Spectrum::LoadFrom(file2->GetDirectory(param_obs_sig.c_str())));
1039  calibshapespect.push_back(*Spectrum::LoadFrom(file6->GetDirectory(param_obs_sig.c_str())));
1040  cherenspect.push_back(*Spectrum::LoadFrom(file3->GetDirectory(param_obs_sig.c_str())));
1041 
1042  calibupspect.push_back(*Spectrum::LoadFrom(file4->GetDirectory(param_obs_bkg.c_str())));
1043  calibdownspect.push_back(*Spectrum::LoadFrom(file5->GetDirectory(param_obs_bkg.c_str())));
1044  lightupspect.push_back(*Spectrum::LoadFrom(file1->GetDirectory(param_obs_bkg.c_str())));
1045  lightdownspect.push_back(*Spectrum::LoadFrom(file2->GetDirectory(param_obs_bkg.c_str())));
1046  calibshapespect.push_back(*Spectrum::LoadFrom(file6->GetDirectory(param_obs_bkg.c_str())));
1047  cherenspect.push_back(*Spectrum::LoadFrom(file3->GetDirectory(param_obs_bkg.c_str())));
1048 
1049  auto truth_nominal = *sMulti_truth.Nominal();
1050  auto truth_genieup = *sMulti_truth.UpperSigma();
1051  auto truth_geniedwn = *sMulti_truth.LowerSigma();
1052  auto truth_calup = *Spectrum::LoadFrom(file4->GetDirectory(param_truth.c_str()));
1053  auto truth_caldwn = *Spectrum::LoadFrom(file5->GetDirectory(param_truth.c_str()));
1054  auto truth_lightup = *Spectrum::LoadFrom(file1->GetDirectory(param_truth.c_str()));
1055  auto truth_lightdwn = *Spectrum::LoadFrom(file2->GetDirectory(param_truth.c_str()));
1056  auto truth_calibshape = *Spectrum::LoadFrom(file6->GetDirectory(param_truth.c_str()));
1057  auto truth_cherenkov = *Spectrum::LoadFrom(file3->GetDirectory(param_truth.c_str()));
1058 
1059  TCanvas * c100 = new TCanvas("c100","c100");
1060  auto signal_multiverse = multiversespec[1].Nominal()->ToTH1(pot);
1061  auto background_multiverse = multiversespec[2].Nominal()->ToTH1(pot);
1062  background_multiverse->Draw("hist");
1063  signal_multiverse->Draw("hist same");
1064  c100->SaveAs("UncertaintyPerBinPlots/debug/nominal_distributions.pdf");
1065  c100->SaveAs("UncertaintyPerBinPlots/root/nominal_distributions.root");
1066  c100->Close();
1067  delete c100;
1068 
1069  std::cout << truth_nominal.POT() << std::endl;
1070 
1071  std::map<std::string,Spectrum> systs_up_bkg;
1072  std::map<std::string,Spectrum> systs_dwn_bkg;
1073  std::map<std::string,Spectrum> systs_up_signal;
1074  std::map<std::string,Spectrum> systs_dwn_signal;
1075  std::map<std::string,Spectrum> systs_up_truth;
1076  std::map<std::string,Spectrum> systs_dwn_truth;
1077 
1078  //Let's plot the error bands first
1079  int makePlotNumber = 0;
1080  std::string makePlotString = "presel_";
1081  //bool shouldPlot = false;
1082  bool shouldPlot = true;
1083  std::string outname;
1084  std::string outname_root;
1085 
1086  //First Up GENIE
1087  systs_up_truth.insert(std::make_pair("genie",truth_genieup));
1088  systs_dwn_truth.insert(std::make_pair("genie",truth_geniedwn));
1089  systs_up_signal.insert(std::make_pair("genie",*multiversespec[1].UpperSigma()));
1090  systs_dwn_signal.insert(std::make_pair("genie",*multiversespec[1].LowerSigma()));
1091  systs_up_bkg.insert(std::make_pair("genie",*multiversespec[2].UpperSigma()));
1092  systs_dwn_bkg.insert(std::make_pair("genie",*multiversespec[2].LowerSigma()));
1093 
1094  if(shouldPlot){
1095  for(uint i=0; i<3; i++){
1096  if (i==0) makePlotString = "presel_";
1097  if (i==1) makePlotString = "preselsig_";
1098  if (i==2) makePlotString = "preselbkg_";
1099  std::map<std::string,Spectrum> genie_up;
1100  genie_up.insert(std::make_pair("genie",*multiversespec[i].UpperSigma()));
1101  std::map<std::string,Spectrum> genie_dwn;
1102  genie_dwn.insert(std::make_pair("genie",*multiversespec[i].LowerSigma()));
1103 
1104  new TCanvas();
1105  //PlotWithSystErrorBand(*(multiversespec[i].Nominal()),
1106  //genie_up, genie_dwn,pot);
1107  outname = "UncertaintyPerBinPlots/debug/" + makePlotString + "genie_errorband.pdf";
1108  gPad->SaveAs(outname.c_str());
1109  outname_root = "UncertaintyPerBinPlots/root/" + makePlotString + "genie_errorband.root";
1110  gPad->SaveAs(outname_root.c_str());
1111  }
1112  }
1113 
1114  //Next Calibration Normalization
1115  systs_up_truth.insert(std::make_pair("calib",truth_calup));
1116  systs_dwn_truth.insert(std::make_pair("calib",truth_caldwn));
1117  systs_up_signal.insert(std::make_pair("calib",calibupspect[1]));
1118  systs_dwn_signal.insert(std::make_pair("calib",calibdownspect[1]));
1119  systs_dwn_bkg.insert(std::make_pair("calib",calibdownspect[2]));
1120  systs_up_bkg.insert(std::make_pair("calib",calibupspect[2]));
1121 
1122  if(shouldPlot){
1123  for(uint i=0; i<3; i++){
1124  if (i==0) makePlotString = "presel_";
1125  if (i==1) makePlotString = "preselsig_";
1126  if (i==2) makePlotString = "preselbkg_";
1127 
1128  std::map<std::string,Spectrum> cal_up;
1129  cal_up.insert(std::make_pair("calib", calibupspect[i]));
1130  std::map<std::string,Spectrum> cal_dwn;
1131  cal_dwn.insert(std::make_pair("calib", calibdownspect[i]));
1132 
1133  new TCanvas();
1134  //PlotWithSystErrorBand(*(multiversespec[i].Nominal()),
1135  // cal_up, cal_dwn,
1136  // pot);
1137  //PlotWithSystErrorBand(*(multiversespec[i].Nominal()),
1138  //calibupspect[i],calibdownspect[i], pot);
1139  outname = "UncertaintyPerBinPlots/debug/" + makePlotString + "calibration_errorband.pdf";
1140  gPad->SaveAs(outname.c_str());
1141  outname_root = "UncertaintyPerBinPlots/root/" + makePlotString + "calibration_errorband.root";
1142  gPad->SaveAs(outname_root.c_str());
1143  new TCanvas();
1144  }
1145  }
1146 
1147  //Next Light Levels
1148 
1149  systs_up_truth.insert(std::make_pair("light",truth_lightup));
1150  systs_dwn_truth.insert(std::make_pair("light",truth_lightdwn));
1151  systs_up_signal.insert(std::make_pair("light",lightupspect[1]));
1152  systs_dwn_signal.insert(std::make_pair("light",lightdownspect[1]));
1153  systs_up_bkg.insert(std::make_pair("light",lightupspect[2]));
1154  systs_dwn_bkg.insert(std::make_pair("light",lightdownspect[2]));
1155 
1156  if(shouldPlot){
1157  for(uint i=0; i<3; i++){
1158  if (i==0) makePlotString = "presel_";
1159  if (i==1) makePlotString = "preselsig_";
1160  if (i==2) makePlotString = "preselbkg_";
1161 
1162  std::map<std::string,Spectrum> light_up;
1163  light_up.insert(std::make_pair("light", lightupspect[i]));
1164  std::map<std::string,Spectrum> light_dwn;
1165  light_dwn.insert(std::make_pair("light", lightdownspect[i]));
1166 
1167  new TCanvas();
1168  // PlotWithSystErrorBand(*(multiversespec[i].Nominal()),
1169  // light_up, light_dwn, pot);
1170  outname = "UncertaintyPerBinPlots/debug/" + makePlotString + "light_level_errorband.pdf";
1171  gPad->SaveAs(outname.c_str());
1172  outname_root = "UncertaintyPerBinPlots/root/" + makePlotString + "light_level_errorband.root";
1173  gPad->SaveAs(outname_root.c_str());
1174  }
1175  }
1176 
1177  //Next CalibrationShape
1178  systs_up_truth.insert(std::make_pair("shape",truth_calibshape));
1179  systs_up_signal.insert(std::make_pair("shape",calibshapespect[1]));
1180  systs_up_bkg.insert(std::make_pair("shape",calibshapespect[2]));
1181 
1182  if(shouldPlot){
1183  for(uint i=0; i<3; i++){
1184  if (i==0) makePlotString = "presel_";
1185  if (i==1) makePlotString = "preselsig_";
1186  if (i==2) makePlotString = "preselbkg_";
1187 
1188  std::map<std::string,Spectrum> shape;
1189  shape.insert(std::make_pair("shape", calibshapespect[i]));
1190 
1191  new TCanvas();
1192  // PlotWithSystErrorBand(*(multiversespec[i].Nominal()),
1193  // shape, shape, pot);
1194  outname = "UncertaintyPerBinPlots/debug/" + makePlotString + "calibration_shape_errorband.pdf";
1195  gPad->SaveAs(outname.c_str());
1196  outname_root = "UncertaintyPerBinPlots/root/" + makePlotString + "calibration_shape_errorband.root";
1197  gPad->SaveAs(outname_root.c_str());
1198  }
1199  }
1200 
1201  //Next Cherenekov
1202  systs_up_truth.insert(std::make_pair("cherenkov",truth_cherenkov));
1203  systs_up_signal.insert(std::make_pair("cherenkov",cherenspect[1]));
1204  systs_up_bkg.insert(std::make_pair("cherenkov",cherenspect[2]));
1205 
1206  if(shouldPlot){
1207  for(uint i=0; i<3; i++){
1208  if (i==0) makePlotString = "presel_";
1209  if (i==1) makePlotString = "preselsig_";
1210  if (i==2) makePlotString = "preselbkg_";
1211  std::map<std::string,Spectrum> cheren;
1212  cheren.insert(std::make_pair("cherenkov", calibshapespect[i]));
1213 
1214  new TCanvas();
1215  // PlotWithSystErrorBand(*(smultiversespec[i].Nominal()),
1216  // cheren, cheren, pot);
1217  outname = "UncertaintyPerBinPlots/debug/" + makePlotString + "cherenekov_errorband.pdf";
1218  gPad->SaveAs(outname.c_str());
1219  outname_root = "UncertaintyPerBinPlots/root/" + makePlotString + "cherenekov_errorband.root";
1220  gPad->SaveAs(outname_root.c_str());
1221  }
1222  }
1223 
1224  TCanvas *c_sel = new TCanvas("c_sel","c_sel");
1225  TH1F* sel_nominal = (TH1F*)multiversespec[0].Nominal()->ToTH1(pot);
1226  TH1F* sel_genie_up = (TH1F*)multiversespec[0].UpperSigma()->ToTH1(pot);
1227  TH1F* sel_genie_down = (TH1F*)multiversespec[0].LowerSigma()->ToTH1(pot);
1228  TH1F* sel_cal_up = (TH1F*)calibupspect[0].ToTH1(pot);
1229  TH1F* sel_cal_down = (TH1F*)calibdownspect[0].ToTH1(pot);
1230  TH1F* sel_light_up = (TH1F*)lightupspect[0].ToTH1(pot);
1231  TH1F* sel_light_down = (TH1F*)lightdownspect[0].ToTH1(pot);
1232  TH1F* sel_shape = (TH1F*)calibshapespect[0].ToTH1(pot);
1233  TH1F* sel_cherenkov = (TH1F*)cherenspect[0].ToTH1(pot);
1234 
1235  sel_nominal->SetLineColor(kBlack);
1236  sel_genie_up->SetLineColor(kRed);
1237  sel_genie_down->SetLineColor(kRed+1);
1238  sel_cal_up->SetLineColor(kGreen);
1239  sel_cal_down->SetLineColor(kGreen+1);
1240  sel_light_up->SetLineColor(kOrange);
1241  sel_light_down->SetLineColor(kOrange+1);
1242  sel_shape->SetLineColor(kMagenta);
1243  sel_cherenkov->SetLineColor(kBlue);
1244 
1245  TLegend *leg = new TLegend(0.5,0.15,0.8,0.35);
1246  leg->AddEntry(sel_nominal,"Nominal");
1247  leg->AddEntry(sel_genie_up,"Genie up");
1248  leg->AddEntry(sel_genie_down,"Genie down");
1249  leg->AddEntry(sel_cal_up,"Cal up");
1250  leg->AddEntry(sel_cal_down,"Cal down");
1251  leg->AddEntry(sel_light_up,"Light up");
1252  leg->AddEntry(sel_light_down,"Light down");
1253  leg->AddEntry(sel_shape,"Shape");
1254  leg->AddEntry(sel_cherenkov,"Cherenkov");
1255 
1256  sel_genie_up->SetTitle("Reco Selected Events");
1257  sel_genie_up->Draw("hist");
1258  sel_genie_down->Draw("histsame");
1259  sel_cal_up->Draw("histsame");
1260  sel_cal_down->Draw("histsame");
1261  sel_light_up->Draw("histsame");
1262  sel_light_down->Draw("histsame");
1263  sel_shape->Draw("histsame");
1264  sel_cherenkov->Draw("histsame");
1265  sel_nominal->Draw("histsame");
1266 
1267  leg->Draw("same");
1268  c_sel->SaveAs("UncertaintyPerBinPlots/selected_comparison.pdf");
1269  c_sel->SaveAs("UncertaintyPerBinPlots/root/selected_comparison.root");
1270  c_sel->Close();
1271  delete c_sel;
1272 
1273  TCanvas *c_sig = new TCanvas("c_sig","c_sig");
1274  TH1F* sig_nominal = (TH1F*)multiversespec[1].Nominal()->ToTH1(pot);
1275  TH1F* sig_genie_up = (TH1F*)multiversespec[1].UpperSigma()->ToTH1(pot);
1276  TH1F* sig_genie_down = (TH1F*)multiversespec[1].LowerSigma()->ToTH1(pot);
1277  TH1F* sig_cal_up = (TH1F*)calibupspect[1].ToTH1(pot);
1278  TH1F* sig_cal_down = (TH1F*)calibdownspect[1].ToTH1(pot);
1279  TH1F* sig_light_up = (TH1F*)lightupspect[1].ToTH1(pot);
1280  TH1F* sig_light_down = (TH1F*)lightdownspect[1].ToTH1(pot);
1281  TH1F* sig_shape = (TH1F*)calibshapespect[1].ToTH1(pot);
1282  TH1F* sig_cherenkov = (TH1F*)cherenspect[1].ToTH1(pot);
1283 
1284  sig_nominal->SetLineColor(kBlack);
1285  sig_genie_up->SetLineColor(kRed);
1286  sig_genie_down->SetLineColor(kRed+1);
1287  sig_cal_up->SetLineColor(kGreen);
1288  sig_cal_down->SetLineColor(kGreen+1);
1289  sig_light_up->SetLineColor(kOrange);
1290  sig_light_down->SetLineColor(kOrange+1);
1291  sig_shape->SetLineColor(kMagenta);
1292  sig_cherenkov->SetLineColor(kBlue);
1293 
1294  sig_genie_up->SetTitle("Reco Signal Events");
1295  sig_genie_up->Draw("hist");
1296  sig_genie_down->Draw("histsame");
1297  sig_cal_up->Draw("histsame");
1298  sig_cal_down->Draw("histsame");
1299  sig_light_up->Draw("histsame");
1300  sig_light_down->Draw("histsame");
1301  sig_shape->Draw("histsame");
1302  sig_cherenkov->Draw("histsame");
1303  sig_nominal->Draw("histsame");
1304 
1305  leg->Draw("same");
1306  c_sig->SaveAs("UncertaintyPerBinPlots/signal_comparison.pdf");
1307  c_sig->SaveAs("UncertaintyPerBinPlots/root/signal_comparison.root");
1308  c_sig->Close();
1309  delete c_sig;
1310 
1311  TCanvas *c_bkg = new TCanvas("c_bkg","c_bkg");
1312  TH1F* bkg_nominal = (TH1F*)multiversespec[2].Nominal()->ToTH1(pot);
1313  TH1F* bkg_genie_up = (TH1F*)multiversespec[2].UpperSigma()->ToTH1(pot);
1314  TH1F* bkg_genie_down = (TH1F*)multiversespec[2].LowerSigma()->ToTH1(pot);
1315  TH1F* bkg_cal_up = (TH1F*)calibupspect[2].ToTH1(pot);
1316  TH1F* bkg_cal_down = (TH1F*)calibdownspect[2].ToTH1(pot);
1317  TH1F* bkg_light_up = (TH1F*)lightupspect[2].ToTH1(pot);
1318  TH1F* bkg_light_down = (TH1F*)lightdownspect[2].ToTH1(pot);
1319  TH1F* bkg_shape = (TH1F*)calibshapespect[2].ToTH1(pot);
1320  TH1F* bkg_cherenkov = (TH1F*)cherenspect[2].ToTH1(pot);
1321 
1322  bkg_nominal->SetLineColor(kBlack);
1323  bkg_genie_up->SetLineColor(kRed);
1324  bkg_genie_down->SetLineColor(kRed+1);
1325  bkg_cal_up->SetLineColor(kGreen);
1326  bkg_cal_down->SetLineColor(kGreen+1);
1327  bkg_light_up->SetLineColor(kOrange);
1328  bkg_light_down->SetLineColor(kOrange+1);
1329  bkg_shape->SetLineColor(kMagenta);
1330  bkg_cherenkov->SetLineColor(kBlue);
1331 
1332  bkg_genie_up->SetTitle("Reco Background Events");
1333  bkg_genie_up->Draw("hist");
1334  bkg_genie_down->Draw("hist same");
1335  bkg_cal_up->Draw("hist same");
1336  bkg_cal_down->Draw("hist same");
1337  bkg_light_up->Draw("hist same");
1338  bkg_light_down->Draw("hist same");
1339  bkg_shape->Draw("hist same");
1340  bkg_cherenkov->Draw("hist same");
1341  bkg_nominal->Draw("hist same");
1342 
1343  leg->Draw("same");
1344  c_bkg->SaveAs("UncertaintyPerBinPlots/background_comparison.pdf");
1345  c_bkg->SaveAs("UncertaintyPerBinPlots/root/background_comparison.root");
1346  c_bkg->Close();
1347  delete c_bkg;
1348 
1349  TH1F* selected_statistical = GetSelectedStatisticalUncertainty(*(multiversespec[0].Nominal()), pot);
1350  TH1F* background_statistical = GetBackgroundStatisticalUncertainty(*(multiversespec[2].Nominal()), mc_pot);
1351  TH1F* denominator = GetDenominator(*(multiversespec[0].Nominal()),
1352  *(multiversespec[2].Nominal()),
1353  pot);
1354  TH1F* background_systematic = GetBackgroundSystematicUncertainty(*(multiversespec[2].Nominal()),denominator, systs_up_bkg, systs_dwn_bkg, pot);
1355  TH1F* efficiency_systematic = GetEfficiencySystematicUncertainty(*(multiversespec[1].Nominal()), truth_nominal, systs_up_signal, systs_dwn_signal, systs_up_truth, systs_dwn_truth, pot);
1356  //TH1F* efficiency_denominator = GetEfficiencyDenominator(*(multiversespec[1].Nominal()), systs_up_truth.at("genie"), pot);
1357  //TH1F* efficiency_denominator = GetEfficiencyDenominator(*(multiversespec[1].Nominal()), truth_nominal, pot);
1358 
1359  TCanvas *c0 = new TCanvas("c0","c0");
1360  c0->SetLeftMargin(0.15);
1361  TH1D* cut_val_eff = multiversespec[1].Nominal()->ToTH1(pot);
1362  cut_val_eff->GetYaxis()->SetTitle("#epsilon and Purity");
1363  cut_val_eff->GetXaxis()->SetTitle(multiversespec[1].Nominal()->ToTH1(pot)->GetXaxis()->GetTitle());
1364  cut_val_eff->SetMaximum(1);
1365  TH1D* nominal_selected = multiversespec[0].Nominal()->ToTH1(pot);
1366  TH1D* nominal_signal = multiversespec[1].Nominal()->ToTH1(pot);
1367  TH1D* cut_val_purity = multiversespec[1].Nominal()->ToTH1(pot);
1368  TH1D* cut_val_signal_stats = multiversespec[1].Nominal()->ToTH1(pot);
1369  for(int ibin = 1; ibin <= nominal_selected->GetXaxis()->GetNbins(); ibin++){
1370  double fNumSelected = nominal_selected->Integral(ibin,nominal_selected->GetXaxis()->GetNbins());
1371  double fNumSelected1 = nominal_signal->Integral(ibin,nominal_selected->GetXaxis()->GetNbins());
1372  if(fNumSelected==0){
1373  fNumSelected=1;
1374  fNumSelected1=0;
1375  }
1376  cut_val_purity->SetBinContent(ibin,fNumSelected1/fNumSelected);
1377  cut_val_eff->SetBinContent(ibin,fNumSelected1/nominal_signal->Integral(1,nominal_selected->GetXaxis()->GetNbins()));
1378  if(fNumSelected1==0)cut_val_signal_stats->SetBinContent(ibin,0);
1379  else cut_val_signal_stats->SetBinContent(ibin,1/sqrt(fNumSelected1));
1380  }
1381  cut_val_eff->Draw("hist");
1382  cut_val_purity->Draw("hist same");
1383  cut_val_signal_stats->Draw("hist same");
1384  outname = "UncertaintyPerBinPlots/eff_and_pur_curves.pdf";
1385  Simulation();
1386  c0->SaveAs(outname.c_str());
1387  outname_root = "UncertaintyPerBinPlots/root/eff_and_pur_curves.root";
1388  c0->SaveAs(outname_root.c_str());
1389  c0->Close();
1390  delete c0;
1391 
1392  TCanvas *c1 = new TCanvas("c1","c1");
1393  c1->SetLeftMargin(0.15);
1394  c1->SetLogy();
1395  selected_statistical->GetYaxis()->SetTitle("#frac{#sqrt{N_{sel}}}{(N_{sel} - N_{bkg})}");
1396  selected_statistical->GetXaxis()->SetTitle(multiversespec[1].Nominal()->ToTH1(pot)->GetXaxis()->GetTitle());
1397  selected_statistical->Divide(denominator);
1398  selected_statistical->Draw("hist");
1399  outname = "UncertaintyPerBinPlots/stat_uncert_selected_contribution.pdf";
1400  Simulation();
1401  c1->SaveAs(outname.c_str());
1402  outname_root = "UncertaintyPerBinPlots/root/stat_uncert_selected_contribution.root";
1403  c1->SaveAs(outname_root.c_str());
1404  c1->Close();
1405  delete c1;
1406 
1407  TCanvas *c2 = new TCanvas("c2","c2");
1408  c2->SetLeftMargin(0.15);
1409  c2->SetLogy();
1410  background_statistical->GetYaxis()->SetTitle("#frac{#sqrt{N_{bkg}}}{(N_{sel} - N_{bkg})}");
1411  background_statistical->GetXaxis()->SetTitle(multiversespec[1].Nominal()->ToTH1(pot)->GetXaxis()->GetTitle());
1412  background_statistical->Divide(denominator);
1413  background_statistical->Draw("hist");
1414  outname = "UncertaintyPerBinPlots/stat_uncert_background_contribution.pdf";
1415  Simulation();
1416  c2->SaveAs(outname.c_str());
1417  outname_root = "UncertaintyPerBinPlots/root/stat_uncert_background_contribution.root";
1418  c2->SaveAs(outname_root.c_str());
1419  c2->Close();
1420  delete c2;
1421 
1422  TCanvas *c3 = new TCanvas("c3","c3");
1423  c3->SetLeftMargin(0.15);
1424  background_systematic->GetYaxis()->SetTitle("#frac{#delta N_{bkg}}{(N_{sel} - N_{bkg})}");
1425  background_systematic->GetXaxis()->SetTitle(multiversespec[1].Nominal()->ToTH1(pot)->GetXaxis()->GetTitle());
1426  background_systematic->Divide(denominator);
1427  background_systematic->SetMaximum(1);
1428  background_systematic->Draw("hist");
1429  Simulation();
1430  outname = "UncertaintyPerBinPlots/syst_uncert_background_contribution.pdf";
1431  c3->SaveAs(outname.c_str());
1432  outname_root = "UncertaintyPerBinPlots/root/syst_uncert_background_contribution.root";
1433  c3->SaveAs(outname_root.c_str());
1434  c3->Close();
1435  delete c3;
1436 
1437  TCanvas *c4 = new TCanvas("c4","c4");
1438  c4->SetLeftMargin(0.15);
1439  efficiency_systematic->GetYaxis()->SetTitle("#frac{#delta #epsilon}{#epsilon}");
1440  efficiency_systematic->GetXaxis()->SetTitle(multiversespec[1].Nominal()->ToTH1(pot)->GetXaxis()->GetTitle());
1441  //efficiency_systematic->Divide(efficiency_denominator);
1442  efficiency_systematic->Divide(truth_nominal.ToTH1(pot));
1443  efficiency_systematic->SetMaximum(1);
1444  efficiency_systematic->Draw("hist");
1445  Simulation();
1446  outname = "UncertaintyPerBinPlots/syst_uncert_efficiency_contribution.pdf";
1447  c4->SaveAs(outname.c_str());
1448  outname_root = "UncertaintyPerBinPlots/root/syst_uncert_efficiency_contribution.root";
1449  c4->SaveAs(outname_root.c_str());
1450  c4->Close();
1451  delete c4;
1452 
1453  TH1D* dsigma = (TH1D*) nominal_signal->Clone("dsigma");
1454  dsigma->SetYTitle("#frac{#delta#sigma}{#sigma}");
1455 
1456  for( int ibin = 1; ibin <= efficiency_systematic->GetXaxis()->GetNbins(); ibin++){
1457  double first = selected_statistical->GetBinContent(ibin);
1458  double second = background_statistical->GetBinContent(ibin);
1459  double third = background_systematic->GetBinContent(ibin);
1460  double fourth = efficiency_systematic->GetBinContent(ibin);
1461 
1462  double final = first*first + second*second + third*third + fourth*fourth;
1463 
1464  final = sqrt(final);
1465 
1466  dsigma->SetBinContent(ibin,final);
1467  }
1468 
1469  TCanvas *c5 = new TCanvas("c5","c5");
1470  c5->SetLeftMargin(0.15);
1471  dsigma->SetMaximum(1);
1472  dsigma->Draw("hist");
1473  Simulation();
1474  outname = "UncertaintyPerBinPlots/xsec_total_uncertainty.pdf";
1475  c5->SaveAs(outname.c_str());
1476  outname_root = "UncertaintyPerBinPlots/root/xsec_total_uncertainty.root";
1477  c5->SaveAs(outname_root.c_str());
1478  c5->Close();
1479  delete c5;
1480 
1481  TCanvas *c6 = new TCanvas("c6","c6");
1482  c6->SetLeftMargin(0.15);
1483  dsigma->SetMaximum(1);
1484  cut_val_eff->GetYaxis()->SetTitle("#epsilon and Purity and #frac{#delta#sigma}{#sigma}");
1485  cut_val_eff->SetLineColor(kRed);
1486  cut_val_eff->Draw("hist");
1487  cut_val_purity->SetLineColor(kCyan);
1488  cut_val_purity->Draw("hist same");
1489  cut_val_signal_stats->SetLineColor(kGreen);
1490  cut_val_signal_stats->Draw("hist same");
1491  dsigma->Draw("hist same");
1492  leg = new TLegend(.18,.45,.35,.75);
1493  leg->AddEntry(cut_val_eff,"#epsilon","l");
1494  leg->AddEntry(cut_val_purity,"Purity","l");
1495  leg->AddEntry(cut_val_signal_stats,"#deltaN_{Sig}^{stat}/N_{Sig}","l");
1496  leg->AddEntry(dsigma,"#delta#sigma/#sigma","l");
1497  leg->SetBorderSize(0); //no border for legend
1498  leg->SetFillColor(0); //fill colour is white
1499  leg->Draw("hist");
1500  Simulation();
1501  outname = "UncertaintyPerBinPlots/eff_pur_xsec_curves.pdf";
1502  c6->SaveAs(outname.c_str());
1503  outname_root = "UncertaintyPerBinPlots/root/eff_pur_xsec_curves.root";
1504  c6->SaveAs(outname_root.c_str());
1505  c6->Close();
1506  delete c6;
1507 
1508  float tmpminval=5;
1509  float tmpminbin=0;
1510  for(int i=1;i<dsigma->GetXaxis()->GetNbins();i++){
1511  if(dsigma->GetBinContent(i)==0)continue;
1512  if(tmpminval < dsigma->GetBinContent(i))continue;
1513  tmpminval=dsigma->GetBinContent(i);
1514  tmpminbin=i;
1515  }
1516  cout<<"The Bin# containing the minimum value is: "<<tmpminbin<<" with a relative uncertainty of "<<dsigma->GetBinContent(tmpminbin)<<endl;
1517 
1518  file1->Close();
1519  file2->Close();
1520  file3->Close();
1521  file4->Close();
1522  file5->Close();
1523  file6->Close();
1524  filerw->Close();
1525 
1526 }
void Simulation()
Definition: tools.h:16
const XML_Char * name
Definition: expat.h:151
enum BeamMode kOrange
TH1F * GetBackgroundSystematicUncertainty(Spectrum spectrum_nominal, TH1F *denominator, std::map< std::string, Spectrum > ShiftUp, std::map< std::string, Spectrum > ShiftDown, double pot_function)
enum BeamMode kRed
constexpr T fourth(T x)
Definition: pow.h:33
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
set< int >::iterator it
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
T sqrt(T number)
Definition: d0nt_math.hpp:156
TLegend * leg1
Definition: plot_hist.C:105
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
float abs(float number)
Definition: d0nt_math.hpp:39
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
c2
Definition: demo5.py:33
TH1F * GetDenominator(Spectrum spectrum_selected, Spectrum spectrum_background, double pot_function)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
Int_t col[ntarg]
Definition: Style.C:29
correl_xv GetXaxis() -> SetDecimals()
std::vector< float > Spectrum
Definition: Constants.h:610
const double mc_pot
TGraphAsymmErrors * MyPlotWithSystErrorBand(TH1F *&nom, std::map< std::string, TH1F * > &ups, std::map< std::string, TH1F * > &dns, int col, int errCol, float headroom, bool newaxis)
OStream cout
Definition: OStream.cxx:6
TH1F * GetEfficiencySystematicUncertainty(Spectrum spec_signal_nominal, Spectrum spec_truth_nominal, std::map< std::string, Spectrum > ShiftUp, std::map< std::string, Spectrum > ShiftDown, std::map< std::string, Spectrum > TruthUp, std::map< std::string, Spectrum > TruthDown, double pot_function)
const bool up_down_combined
const Spectrum * UpperSigma(BandOptions opt=kBandFromNominal) const
c1
Definition: demo5.py:24
const Spectrum * LowerSigma(BandOptions opt=kBandFromNominal) const
const Color_t kTotalMCColor
Definition: Style.h:16
void xsec_uncertainty_per_bin(std::string param="energy")
TH1F * GetSelectedStatisticalUncertainty(Spectrum spectrum_nominal, double pot_function)
enum BeamMode kGreen
enum BeamMode kBlue
TH1F * GetBackgroundStatisticalUncertainty(Spectrum spectrum_nominal, double pot_function)
const double pot
const Spectrum * Nominal() const
return the nominal universe
enum BeamMode string
unsigned int uint